EM algorithm for Gaussian mix Can anyone help me with the R code to implement EM algorithm.  I got different value if I chose different starting value; clearly this is not good. And the value of $\mu$, $\sigma$ goes to NA after few iterations. Here is my code
da1=read.table("anythin.Rdata", header=TRUE)
y=as.vector(da1[,2])
n=length(y)
mu=matrix(NA,1000,2)
sigma=matrix(NA,1000,2)
w=matrix(NA,1000,2)
mu[1,]=c(2,4)
sigma[1,]=c(0.5,0.1)
w[1,]=c(0.5,0.5)
xi1=0
xi2=0
for (i in 2:1000){

    # E step
    xi1=w[i-1,1]*dnorm(y,mean=mu[i-1,1],sd=sigma[i-1,1])/(w[i-1,1]*dnorm(y,mean=mu[i-1,1],sd=sigma[i-1,1])+w[i-1,2]*dnorm(y,mean=mu[i-1,2],sd=sigma[i-1,2]))
    xi2=w[i-1,2]*dnorm(y,mean=mu[i-1,2],sd=sigma[i-1,2])/(w[i-1,1]*dnorm(y,mean=mu[i-1,1],sd=sigma[i-1,1])+w[i-1,2]*dnorm(y,mean=mu[i-1,2],sd=sigma[i-1,2]))

    # M step
    w[i,1]=sum(xi1)
    w[i,2]=sum(xi2)

    mu[i,1]=sum(xi1*y)/sum(xi1)
    mu[i,2]=sum(xi2*y)/sum(xi2)

    sigma[i,1]=sum(xi1*(mu[i,1]-y)^2)/sum(xi1)
    sigma[i,2]=sum(xi2*(mu[i,2]-y)^2)/sum(xi2)
}

 A: "Can anyone help me with my code" is a StackOverflow question.
@StasK says "this may be data specific".
It IS data specific. EM is a gradient descent. It isn't "global optimal" but "local optimal". Just like Newtons method it can be launched into the boonies, especially with home-brew formulations not built to resist it.
The simplest EM for a Gaussian Mixture is actually k-means.  It assumes constant variance, so there isn't a variance parameter update.
K-means works like this:


*

*Specify number of clusters

*Randomly assign points to a cluster

*Using cluster membership compute cluster centers, assign to clusters

*Using point position and cluster center, assign points to nearest cluster

*Repeat 3 and 4 until "convergence"


This is the "happy path", the ideal case.  It is incompatible with reality.  In reality there are 24 unhappy paths for each happy path.
Here are two pretty common "unhappy paths" for the GMM


*

*you can get oscillations.  It can get down as "tight" as data allows and then oscillate between values

*you can get a single point in a cluster and send the variance toward zero. 

