Algoritmo EM para la mezcla gaussiana

Question

Puede alguien ayudarme con el código R para implementar el algoritmo EM. Obtuve un valor diferente si elegí un valor inicial diferente; claramente esto no es bueno. Y el valor de $\mu$ , $\sigma$ pasa a NA después de algunas iteraciones. Este es mi código

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)
}

Answer 1

"¿Puede alguien ayudarme con mi código?" es una pregunta de StackOverflow.

@StasK dice que "esto puede ser un dato específico".
Es ES datos específicos. EM es un descenso de gradiente. No es "óptimo global" sino "óptimo local". Al igual que el método de Newton, puede lanzarse al vacío, especialmente con formulaciones caseras que no están hechas para resistirlo.

El EM más sencillo para una mezcla gaussiana es en realidad k-means. Asume una varianza constante, por lo que no hay una actualización del parámetro de varianza.

K-means funciona así:

1. Especifique el número de clusters
2. Asignación aleatoria de puntos a una agrupación
3. Utilizando la membresía de los clusters, calcular los centros de los clusters, asignarlos a los clusters
4. Utilizando la posición de los puntos y el centro del cluster, asigna los puntos al cluster más cercano
5. Repetir 3 y 4 hasta la "convergencia"

Este es el "camino feliz", el caso ideal. Es incompatible con la realidad. En la realidad hay 24 caminos infelices por cada camino feliz.

Aquí hay dos "caminos infelices" bastante comunes para el MMG

• puede obtener oscilaciones. Puede bajar tan "apretado" como los datos lo permitan y luego oscilar entre los valores
• puede obtener un solo punto en un grupo y enviar la varianza hacia cero.