Algoritmo EM implementado manualmente

Quiero implementar el algoritmo EM manualmente y luego compararlo con los resultados normalmixEMdel mixtoolspaquete. Por supuesto, sería feliz si ambos conducen a los mismos resultados. La referencia principal es Geoffrey McLachlan (2000), Modelos de mezclas finitas .

Tengo una densidad de mezcla de dos gaussianos, en forma general, la probabilidad logarítmica está dada por (McLachlan página 48):

\log L_{c} (Ψ) = \sum_{i = 1}^{g} \sum_{j = 1}^{n} z_{i j} {\log π_{i} + \log f_{i} (y_{i}; θ_{i})} .

$\log L_c(\Psi) = \sum_{i=1}^g \sum_{j=1}^n z_{ij}\{\log \pi_i + \log f_i(y_i;\theta_i)\}.$ Los son , si la observación fue de la densidad de componentes^ésima , de lo contrario

. El

es la densidad de la distribución normal. La

es la proporción de la mezcla, por lo que

es la probabilidad de que una observación sea de la primera distribución gaussiana y

es la probabilidad de que una observación sea de la segunda distribución gaussiana.

z_{i j}

$z_{ij}$

1

$1$

i

$i$

0

$0$

f_{i}

$f_i$

π

$\pi$

π_{1}

$\pi_1$

π_{2}

$\pi_2$

El paso E ahora es el cálculo de la expectativa condicional:

Q (Ψ; Ψ^{(0)}) = E_{Ψ (0)} {\log L_{c} (| Ψ) | y} .

$Q(\Psi;\Psi^{(0)}) = E_{\Psi(0)}\{\log L_c(|\Psi)|y\}.$ que conduce, después de algunas derivaciones al resultado (página 49):

\begin{aligned} τ_{i} (y_{j}; Ψ^{(k)}) & = \frac{π_{i}^{(k)} f_{i} (y_{j}; θ_{i}^{(k)}}{f (y_{j}; Ψ^{(k)}} \\ = \frac{π_{i}^{(k)} f_{i} (y_{j}; θ_{i}^{(k)}}{\sum_{h = 1}^{g} π_{h}^{(k)} f_{h} (y_{j}; θ_{h}^{(k)})} \end{aligned}

$\begin{align} \tau_i(y_j;\Psi^{(k)}) &= \frac{\pi_i^{(k)}f_i(y_j;\theta_i^{(k)}}{f(y_j;\Psi^{(k)}} \\[8pt] &= \frac{\pi_i^{(k)}f_i(y_j;\theta_i^{(k)}}{\sum_{h=1}^g \pi_h^{(k)}f_h(y_j;\theta_h^{(k)})} \end{align}$ en el caso de dos gaussianos (página 82):

τ_{i} (y_{j}; Ψ) = \frac{π_{i} ϕ (y_{j}; μ_{i}, Σ_{i})}{\sum_{h = 1}^{g} π_{h} ϕ (y_{j}; μ_{h}, Σ_{h})}

$\tau_i(y_j;\Psi) = \frac{\pi_i \phi(y_j;\mu_i,\Sigma_i)}{\sum_{h=1}^g \pi_h\phi(y_j; \mu_h,\Sigma_h)}$ Elpaso M ahora es la maximización de Q (página 49):

Q (Ψ; Ψ^{(k)}) = \sum_{i = 1}^{g} \sum_{j = 1}^{n} τ_{i} (y_{j}; Ψ^{(k)}) {\log π_{i} + \log f_{i} (y_{j}; θ_{i})} .

$Q(\Psi;\Psi^{(k)}) = \sum_{i=1}^g\sum_{j=1}^n\tau_i(y_j;\Psi^{(k)})\{\log \pi_i + \log f_i(y_j;\theta_i)\}.$ Esto lleva a (en el caso de dos gaussianos) (página 82):

\begin{aligned} μ_{i}^{(k + 1)} & = \frac{\sum_{j = 1}^{n} τ_{i j}^{(k)} y_{j}}{\sum_{j = 1}^{n} τ_{i j}^{(k)}} \\ Σ_{i}^{(k + 1)} & = \frac{\sum_{j = 1}^{n} τ_{i j}^{(k)} (y_{j} - μ_{i}^{(k + 1)}) (y_{j} - μ_{i}^{(k + 1)})^{T}}{\sum_{j = 1}^{n} τ_{i j}^{(k)}} \end{aligned}

$\begin{align} \mu_i^{(k+1)} &= \frac{\sum_{j=1}^n \tau_{ij}^{(k)}y_j}{\sum_{j=1}^n \tau_{ij}^{(k)}} \\[8pt] \Sigma_i^{(k+1)} &= \frac{\sum_{j=1}^n \tau_{ij}^{(k)}(y_j - \mu_i^{(k+1)})(y_j - \mu_i^{(k+1)})^T}{\sum_{j=1}^n \tau_{ij}^{(k)}} \end{align}$ y sabemos que (p. 50)

π_{i}^{(k + 1)} = \frac{\sum_{j = 1}^{n} τ_{i} (y_{j}; Ψ^{(k)})}{n} (i = 1, \dots, g) .

$\pi_i^{(k+1)} = \frac{\sum_{j=1}^n \tau_i(y_j;\Psi^{(k)})}{n}\qquad (i = 1, \ldots, g).$ Repetimos los pasos E, M hasta que es pequeño.

L (Ψ^{(k + 1)}) - L (Ψ^{(k)})

$L(\Psi^{(k+1)})-L(\Psi^{(k)})$

Traté de escribir un código R (los datos se pueden encontrar aquí ).

# 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

El algoritmo no funciona, ya que algunas observaciones tienen la probabilidad de cero y el registro de esto sí -Inf. ¿Dónde está mi error?

r expectation-maximization gaussian-mixture Stat Tistician
fuente

El problema no es estadístico, sino numérico. Debe agregar contingencias para probabilidades más pequeñas que la precisión de la máquina en su código.

JohnRos

¿Por qué no intentas probar la función mixtools con un ejemplo muy simple que se puede verificar a mano, digamos solo cinco o diez valores y dos series de tiempo, primero. luego, si encuentra que funciona allí, generalice su código y verifique en cada paso.

Respuestas:

Tiene varios problemas en el código fuente:

Como señaló @Pat, no debe usar log (dnorm ()) ya que este valor puede llegar fácilmente al infinito. Deberías usar logmvdnorm
Cuando use la suma , tenga en cuenta que elimina los valores infinitos o faltantes
Su variable de bucle k está mal, debe actualizar loglik [k + 1] pero actualiza loglik [k]
Los valores iniciales para su método y mixtools son diferentes. Está usando en su método, pero está usando para mixtools (es decir, desviación estándar, del manual de mixtools). $\Sigma$ $\sigma$
Sus datos no se ven como una mezcla de lo normal (verifique el histograma que tracé al final). Y un componente de la mezcla tiene un SD muy pequeño, así que agregué arbitrariamente una línea para configurar y para que sean iguales para algunas muestras extremas. Los agrego solo para asegurarme de que el código pueda funcionar. $\tau_1$ $\tau_2$

También le sugiero que ponga códigos completos (p. Ej., Cómo inicializa loglik []) en su código fuente y sangría el código para que sea fácil de leer.

Después de todo, gracias por presentar el paquete mixtools , y planeo usarlos en mi investigación futura.

También pongo mi código de trabajo para su referencia:

# EM algorithm manually
# dat is the data
setwd("~/Downloads/")
load("datem.Rdata")
x <- dat

# initial values
pi1<-0.5
pi2<-0.5
mu1<--0.01
mu2<-0.01
sigma1<-sqrt(0.01)
sigma2<-sqrt(0.02)
loglik<- rep(NA, 1000)
loglik[1]<-0
loglik[2]<-mysum(pi1*(log(pi1)+log(dnorm(dat,mu1,sigma1))))+mysum(pi2*(log(pi2)+log(dnorm(dat,mu2,sigma2))))

mysum <- function(x) {
  sum(x[is.finite(x)])
}
logdnorm <- function(x, mu, sigma) {
  mysum(sapply(x, function(x) {logdmvnorm(x, mu, sigma)}))  
}
tau1<-0
tau2<-0
#k<-1
k<-2

# loop
while(abs(loglik[k]-loglik[k-1]) >= 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))
  tau1[is.na(tau1)] <- 0.5
  tau2[is.na(tau2)] <- 0.5

  # M step
  pi1<-mysum(tau1)/length(dat)
  pi2<-mysum(tau2)/length(dat)

  mu1<-mysum(tau1*x)/mysum(tau1)
  mu2<-mysum(tau2*x)/mysum(tau2)

  sigma1<-mysum(tau1*(x-mu1)^2)/mysum(tau1)
  sigma2<-mysum(tau2*(x-mu2)^2)/mysum(tau2)

  #  loglik[k]<-sum(tau1*(log(pi1)+log(dnorm(x,mu1,sigma1))))+sum(tau2*(log(pi2)+log(dnorm(x,mu2,sigma2))))
  loglik[k+1]<-mysum(tau1*(log(pi1)+logdnorm(x,mu1,sigma1)))+mysum(tau2*(log(pi2)+logdnorm(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

Historgram Histograma

zhanxw
fuente

@zahnxw gracias por tu respuesta, ¿eso significa que mi código está equivocado? ¿Entonces la idea básica no está funcionando?

Stat Tistician

"También le sugiero que ponga códigos completos (p. Ej., Cómo inicializa loglik []) en su código fuente y sangra el código para que sea fácil de leer". Bueno, este es mi código? el loglik [] se define como lo declare en el código que publiqué?

Stat Tistician

@StatTistician, la idea es correcta, pero la implementación tiene fallas. Por ejemplo, no consideró el flujo insuficiente. Además, su variable de bucle k es confusa, primero configura loglik [1] y loglik [2], después de ingresar al bucle while, configura nuevamente loglik [1]. Esta no es la forma natural de hacerlo. Mi sugerencia sobre la inicialización de loglik [] significa código:, loklik <- rep(NA, 100)que asignará previamente loglik [1], loglik [2] ... loglik [100]. Planteo esa pregunta porque en su código original, no encontré la delcaración de loglik, ¿tal vez el código se trunca durante el pegado?

zhanxw

Como publiqué a continuación: Gracias por su ayuda, pero estoy abandonando este tema, ya que es demasiado avanzado para mí.

Stat Tistician

¿Existe ahora una manera de determinar qué parte de los datos pertenece a qué mezcla?

Cardenal

Sigo recibiendo un error al intentar abrir su archivo .rar, pero eso puede ser solo que estoy haciendo algo tonto.

$f(y;\theta)$ $\exp(-0.5(y-\mu)^2/\sigma^2)$ $\mu$ $y$ $\tau$

Si ese es el problema, hay algunas soluciones posibles:

$\tau$

$\tau \log(f(y|\theta))$

evaluar

$\log \left( f(y|\theta)^\tau \right)$

$f(y|\theta)$ $\tau$ $\approx 0$

$0 \log (0) = 0 (-Inf) = NaN$

pero con tau movido obtienes

$\log \left( 0^0\right) = \log(1) = 0$

$0^0 = 1$

Otra solución es expandir las cosas dentro del logaritmo. Asumiendo que estás usando logaritmos naturales:

$\tau \log(f(y|\theta))$

$= \tau \log(\exp(-0.5(y-\mu)^2/\sigma^2)/\sqrt{2\pi\sigma^2})$

$= -0.5\tau \log(2 \pi\sigma^2) - 0.5 \tau \frac{(y-\mu)^2}{\sigma^2}$

Matemáticamente es lo mismo, pero debería ser más resistente a los errores de coma flotante, ya que ha evitado calcular una gran potencia negativa. Esto significa que ya no puede usar la función de evaluación de normas incorporada, pero si eso no es un problema, esta es probablemente la mejor respuesta. Por ejemplo, digamos que tenemos la situación donde

$-0.5\frac{(y-\mu)^2}{\sigma^2} = -0.5*40^2 = -800$

$\log(\exp(-800)) = \log(0) = -Inf$

Palmadita
fuente

mh, para ser sincero: no soy lo suficientemente bueno para hacer que esto funcione. Lo que me interesaba es: ¿Puedo obtener el mismo resultado con mi algoritmo que la versión implementada del paquete mixtools. Pero desde mi punto de vista, esto parece estar pidiendo la luna. Pero creo que pones esfuerzo en tu respuesta, ¡así que lo aceptaré! ¡Gracias!

Stat Tistician