Quiero implementar el algoritmo EM manualmente y luego compararlo con los resultados del normalmixEM de mixtools paquete. Por supuesto, me alegraría que ambos condujeran 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 log-verosimilitud viene dada por (McLachlan página 48):
$$ \log L_c(\Psi) = \sum_{i=1}^g \sum_{j=1}^n z_{ij}\{\log \pi_i + \log f_i(y_i;\theta_i)\}. $$ En $z_{ij}$ son $1$ si la observación procedía del $i$ th densidad de componentes, en caso contrario $0$ . En $f_i$ es la densidad de la distribución normal. La dirección $\pi$ es la proporción de la mezcla, por lo que $\pi_1$ es la probabilidad de que una observación proceda de la primera distribución gaussiana y $\pi_2$ es la probabilidad de que una observación proceda de la segunda distribución gaussiana.
En E paso es ahora, el cálculo de la expectativa condicional:
$$ Q(\Psi;\Psi^{(0)}) = E_{\Psi(0)}\{\log L_c(|\Psi)|y\}. $$ que conduce, tras algunas derivaciones al resultado (página 49):
\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 gaussianas (página 82):
$$ \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)} $$ En M es ahora la maximización de Q (página 49):
$$ 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 conduce a (en el caso de dos gaussianas) (página 82):
\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)
$$ \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 $L(\Psi^{(k+1)})-L(\Psi^{(k)})$ es pequeño.
He intentado 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$loglikEl algoritmo no funciona, ya que algunas observaciones tienen probabilidad cero y el logaritmo de ésta es -Inf . ¿Dónde está mi error?
0 votos
El problema no es estadístico, sino numérico. Deberías añadir contingencias para probabilidades menores que la precisión de la máquina en tu código.
0 votos
¿por qué no pruebas a utilizar la función mixtools con un ejemplo muy simple que pueda verificarse a mano, por ejemplo, sólo cinco o diez valores y dos series temporales, primero. luego, si ves que funciona, generaliza tu código y verifícalo en cada paso.