I want to implement the EM algorithm manually, with my own loops and so.
Afterwards, I want to compare it to the normalmixEM output of mixtools
package.
Since the notation is very advanced, I used LaTex and attached the two
screenshots I also attached the data. Alternatively, the data can be found
here:http://uploadeasy.net/upload/py6j4.rar and the screenshots here:
http://uploadeasy.net/upload/9i04.PNG
http://uploadeasy.net/upload/vloku.PNG
I want to fit two gaussian mixtures.
My main question is: Did I implement the EM algorithm correctly? It does
not work, because I get such odd values especially for sigma, that the
likelihood of some observations is zero and therefore the log is -Inf.
Where is my mistake?
My code:
# EM algorithm manually
# dat is the data
# initial values
pi1<-0.5
pi2<-0.5
mu1<--0.01
mu2<-0.01
sigma1<-0.01
sigma2<-0.02
loglik[1]<-0
loglik[2]<-sum(pi1*(log(pi1)+log(dnorm(dat,mu1,sigma1))))+sum(pi2*(log(pi2)+log(dnorm(dat,mu2,sigma2))))
tau1<-0
tau2<-0
k<-1
# loop
while(abs(loglik[k+1]-loglik[k])>= 0.00001) {
# E step
tau1<-pi1*dnorm(dat,mean=mu1,sd=sigma1)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))
tau2<-pi2*dnorm(dat,mean=mu2,sd=sigma2)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))
# M step
pi1<-sum(tau1)/length(dat)
pi2<-sum(tau2)/length(dat)
mu1<-sum(tau1*x)/sum(tau1)
mu2<-sum(tau2*x)/sum(tau2)
sigma1<-sum(tau1*(x-mu1)^2)/sum(tau1)
sigma2<-sum(tau2*(x-mu2)^2)/sum(tau2)
loglik[k]<-sum(tau1*(log(pi1)+log(dnorm(x,mu1,sigma1))))+sum(tau2*(log(pi2)+log(dnorm(x,mu2,sigma2))))
k<-k+1
}
# compare
library(mixtools)
gm<-normalmixEM(x,k=2,lambda=c(0.5,0.5),mu=c(-0.01,0.01),sigma=c(0.01,0.02))
gm$lambda
gm$mu
gm$sigma
gm$loglik
-------------- next part --------------
A non-text attachment was scrubbed...
Name: no1..PNG
Type: image/png
Size: 93293 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20130404/528e57ad/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: no2.PNG
Type: image/png
Size: 46672 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20130404/528e57ad/attachment-0001.png>