# [R] EM algorithm implemented manually

Stat Tistician statisticiangermany at gmail.com
Thu Apr 4 13:44:11 CEST 2013

```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

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