Algorithme EM implémenté manuellement

Je veux implémenter l'algorithme EM manuellement, puis le comparer aux résultats normalmixEMdu mixtoolspackage. Bien sûr, je serais heureux si les deux aboutissaient aux mêmes résultats. La référence principale est Geoffrey McLachlan (2000), Finite Mixture Models .

J'ai une densité de mélange de deux gaussiens, sous forme générale, la log-vraisemblance est donnée par (McLachlan page 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)\}.$$ Les sont , si l'observation est de la ^ème densité de composantes, sinon . Le est la densité de la distribution normale. Le est la proportion du mélange, donc est la probabilité, qu'une observation soit de la première distribution gaussienne et est la probabilité, qu'une observation soit de la deuxième distribution gaussienne.$z_{ij}$$1$$i$$0$$f_i$$\pi$$\pi_1$$\pi_2$

L' étape E est maintenant le calcul de l'espérance conditionnelle:

$$Q(\Psi;\Psi^{(0)}) = E_{\Psi(0)}\{\log L_c(|\Psi)|y\}.$$ ce qui conduit, après quelques dérivations au résultat (page 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}$$ dans le cas de deux Gaussiens (page 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)}$$ L' étape M est maintenant la maximisation de Q (page 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)\}.$$ Cela conduit à (dans le cas de deux Gaussiens) (page 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}$$ et nous le savons (p. 50)

$$\pi_i^{(k+1)} = \frac{\sum_{j=1}^n \tau_i(y_j;\Psi^{(k)})}{n}\qquad (i = 1, \ldots, g).$$ Nous répétons les étapes E, M jusqu'à ce que$L(\Psi^{(k+1)})-L(\Psi^{(k)})$ soit petit.

J'ai essayé d'écrire un code R (les données peuvent être trouvées ici ).

# 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

L'algorithme ne fonctionne pas, car certaines observations ont la probabilité de zéro et le logarithme de cela est -Inf. Où est mon erreur?

Vous avez plusieurs problèmes dans le code source:

1. Comme l'a souligné @Pat, vous ne devez pas utiliser log (dnorm ()) car cette valeur peut facilement aller à l'infini. Vous devez utiliser logmvdnorm

2. Lorsque vous utilisez sum , n'oubliez pas de supprimer les valeurs infinies ou manquantes

3. La variable k en boucle est incorrecte, vous devez mettre à jour loglik [k + 1] mais vous mettez à jour loglik [k]

4. Les valeurs initiales de votre méthode et de vos mixtools sont différentes. Vous utilisez dans votre méthode, mais utilisez pour mixtools (c.-à-d. Écart type, du manuel mixtools).$\Sigma$$\sigma$

5. Vos données ne ressemblent pas à un mélange de normal (vérifiez l'histogramme que j'ai tracé à la fin). Et un composant du mélange a un très petit sd, j'ai donc arbitrairement ajouté une ligne pour définir et pour qu'ils soient égaux pour certains échantillons extrêmes. Je les ajoute juste pour m'assurer que le code peut fonctionner.τ $\tau_1$$\tau_2$

Je vous suggère également de mettre des codes complets (par exemple, comment vous initialisez loglik []) dans votre code source et de mettre le code en retrait pour en faciliter la lecture.

Après tout, merci d'avoir introduit le package mixtools , et je prévois de les utiliser dans mes futures recherches.

J'ai également mis mon code de travail pour référence:

# 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

Je reçois toujours une erreur lorsque j'essaie d'ouvrir votre fichier .rar, mais cela peut juste être moi qui fais quelque chose de stupide.

Je ne vois pas d'erreurs évidentes dans votre code. Une raison possible pour laquelle vous obtenez des zéros est due à la précision en virgule flottante. N'oubliez pas que lorsque vous calculez , vous évaluez . Cela ne prend pas une très grande différence entre et pour que cela soit arrondi à 0 lorsque vous le faites sur un ordinateur. Ceci est doublement perceptible dans les modèles de mélange, car certaines de vos données ne seront pas "affectées" à chaque composant du mélange et peuvent donc se retrouver très loin de lui. En théorie, ces points devraient également se retrouver avec une faible valeur de$f(y;\theta)$$\exp(-0.5(y-\mu)^2/\sigma^2)$$\mu$$y$$\tau$ lorsque vous évaluez la probabilité du journal, neutralisant le problème - mais grâce à l'erreur en virgule flottante, la quantité a déjà été évaluée comme -Inf à ce stade, donc tout se rompt :).

Si tel est le problème, il existe des solutions possibles:

L'une consiste à déplacer votre à l'intérieur du logarithme. Donc, au lieu d'évaluer$\tau$

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

évaluer

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

Mathématiquement la même chose, mais pensez à ce qui se passe lorsque et sont . Actuellement, vous obtenez:$f(y|\theta)$$\tau$$\approx 0$

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

mais avec tau déplacé vous obtenez

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

en supposant que R évalue (je ne sais pas si c'est le cas ou j'ai tendance à utiliser matlab)$0^0 = 1$

Une autre solution consiste à développer le contenu à l'intérieur du logarithme. En supposant que vous utilisez des logarithmes naturels:

$\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}$ .

Mathématiquement identique, mais devrait être plus résistant aux erreurs en virgule flottante car vous avez évité de calculer une grande puissance négative. Cela signifie que vous ne pouvez plus utiliser la fonction d'évaluation de norme intégrée, mais si ce n'est pas un problème, c'est probablement la meilleure réponse. Par exemple, disons que nous avons la situation où

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

Évaluez cela comme je l'ai suggéré, et vous obtenez -800. Cependant, dans matlab, si nous exp le prenons le journal, nous obtenons .$\log(\exp(-800)) = \log(0) = -Inf$