Mixtures in Madrid (3)

For my second lecture today, I need to plot a likelihood surface for a basic two-component mixture with only the means unknown: here is the R code to speed up things

llsurf=function(trumyn=2.,wayt=.3,var2=1.,ssiz=500){
# draws the log-likelihood surface and a random sample
sd2=sqrt(var2)
parti=(runif(ssiz)>wayt)
sampl=(1-parti)*rnorm(ssiz)+parti*(trumyn+sd2*rnorm(ssiz))
mu2=mu1=seq(min(sampl),max(sampl),.1)
mo1=mu1%*%t(rep(1,length(mu2)))
mo2=(rep(1,length(mu2)))%*%t(mu2)
ca1=-0.5*mo1*mo1
ca2=-0.5*mo2*mo2
like=.1*(ca1+ca2) # log prior N(0,10)
for (i in 1:ssiz)
like=like+log(wayt*dnorm(sampl[i],mo1,sd2)+(1-wayt)*dnorm(sampl[i],mo2,1))
par(mar=c(4,4,1,1))
image(mu1,mu2,like,xlab=expression(mu[1]),ylab=expression(mu[2]),col=heat.colors(250))
contour(mu1,mu2,like,levels=seq(min(like),max(like),(max(like)-min(like))/50),add=T)
sampl
}

resulting in an outcome (very) similar to the above. The EM steps are then an easy

mu=cmu=rnorm(2,sd=5)
for (t in 1:30){
 z=0.3*dnorm(sampl,mu[1])/(0.3*dnorm(sampl,mu[1])+0.7*dnorm(sampl,mu[2]))
 cmu=rbind(cmu,mu)
 mu[1]=sum(z*sampl)/sum(z)
 mu[2]=sum((1-z)*sampl)/sum(1-z)
 }
lines(cmu[,1],cmu[,2])

while the Gibbs sampler is

mu=cmu=rnorm(2,sd=5)
for (t in 1:10^3){
 z=0.3*dnorm(sampl,mu[1])/(0.3*dnorm(sampl,mu[1])+0.7*dnorm(sampl,mu[2]))
 u=runif(length(z))
 mu[1]=sum((u<z)*sampl)/(.01+sum(u<z))+rnorm(1)/sqrt(.01+sum(u<z))
 mu[2]=sum((u>z)*sampl)/(.01+sum(u>z))+rnorm(1)/sqrt(.01+sum(u>z))
 cmu=rbind(cmu,mu)
}
lines(cmu[,1],cmu[,2])

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.