[R] [R-SIG-Finance] EM algorithm with R manually implemented?

[R. Michael Weylandt]

Moved to R-help because there's no obvious financial content.

Michael

On Sat, Apr 6, 2013 at 10:56 AM, Stat Tistician
<statisticiangermany at gmail.com> wrote:
> Hi,
> 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
>
>
> Thanks
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-SIG-Finance at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.