[R] EM algorithm implemented manually

[Stat Tistician]

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>