Algorithme EM implémenté manuellement

20

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

logLc(Ψ)=i=1gj=1nzij{logπi+logfi(yi;θ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.zij1i0fiππ1π2

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

Q(Ψ;Ψ(0))=EΨ(0){logLc(|Ψ)|y}.
ce qui conduit, après quelques dérivations au résultat (page 49):

τi(yj;Ψ(k))=πi(k)fi(yj;θi(k)f(yj;Ψ(k)=πi(k)fi(yj;θi(k)h=1gπh(k)fh(yj;θh(k))
dans le cas de deux Gaussiens (page 82):

τi(yj;Ψ)=πiϕ(yj;μi,Σi)h=1gπhϕ(yj;μh,Σh)
L' étape M est maintenant la maximisation de Q (page 49):

Q(Ψ;Ψ(k))=je=1gj=1nτje(yj;Ψ(k)){Journalπje+JournalFje(yj;θje)}.
Cela conduit à (dans le cas de deux Gaussiens) (page 82):

μje(k+1)=j=1nτjej(k)yjj=1nτjej(k)Σje(k+1)=j=1nτjej(k)(yj-μje(k+1))(yj-μje(k+1))Tj=1nτjej(k)
et nous le savons (p. 50)

πje(k+1)=j=1nτje(yj;Ψ(k))n(je=1,,g).
Nous répétons les étapes E, M jusqu'à ce queL(Ψ(k+1))-L(Ψ(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?

Stat Tistician
la source
Le problème n'est pas d'ordre statistique, mais plutôt numérique. Vous devez ajouter des contingences pour les vraisemblances inférieures à la précision de la machine dans votre code.
JohnRos
pourquoi n'essayez-vous pas de tester la fonction mixtools avec un exemple très simple qui peut être vérifié à la main, par exemple cinq ou dix valeurs et deux séries de temps, tout d'abord. puis, si vous trouvez que cela fonctionne, généralisez votre code et vérifiez à chaque étape.

Réponses:

17

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).σΣσ

  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.τ 2τ1τ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 Histogramme

zhanxw
la source
@zahnxw merci pour votre réponse, cela signifie-t-il donc que mon code est incorrect? L'idée basi ne fonctionne donc pas?
Stat Tistician
"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." Eh bien, c'est mon code? le loglik [] est défini comme je l'ai déclaré dans le code que j'ai posté?
Stat Tistician
1
@StatTistician l'idée est correcte, mais l'implémentation a des défauts. Par exemple, vous n'avez pas envisagé de sous-débit. De plus, la variable k en boucle est source de confusion, vous définissez d'abord loglik [1] et loglik [2], après avoir entré la boucle while, vous définissez à nouveau loglik [1]. Ce n'est pas la façon naturelle de procéder. Ma suggestion concernant l'initialisation de loglik [] signifie code:, loklik <- rep(NA, 100)qui pré-allouera loglik [1], loglik [2] ... loglik [100]. Je pose cette question car dans votre code d'origine, je n'ai pas trouvé la délocalisation de loglik, peut-être que le code est tronqué lors du collage?
zhanxw
Comme je l'ai posté ci-dessous: Merci pour votre aide, mais je laisse tomber ce sujet, car il est trop avancé pour moi.
Stat Tistician
Existe-t-il maintenant un moyen de déterminer quelle partie des données appartient à quel mélange?
Cardinal
2

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 def(y;θ)exp(0.5(yμ)2/σ2)μyτ 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τ

τJournal(F(y|θ))

évaluer

Journal(F(y|θ)τ) .

Mathématiquement la même chose, mais pensez à ce qui se passe lorsque et sont . Actuellement, vous obtenez:F(y|θ)τ0

  • 0Journal(0)=0(-jenF)=NuneN

mais avec tau déplacé vous obtenez

  • Journal(00)=Journal(1)=0

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

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

τJournal(F(y|θ))

=τJournal(exp(-0,5(y-μ)2/σ2)/2πσ2)

=-0,5τJournal(2πσ2)-0,5τ(y-μ)2σ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(y-μ)2σ2=-0,5402=-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 .Journal(exp(-800))=Journal(0)=-jenF

Tapoter
la source
mh, pour être honnête: je ne suis pas assez bon pour faire fonctionner cette chose. Ce qui m'intéressait, c'est: puis-je obtenir le même résultat avec mon algorithme que la version implémentée du package mixtools. Mais de mon point de vue, cela semble demander la lune. Mais je pense que vous avez fait des efforts dans votre réponse, donc je l'accepterai! Merci!
Stat Tistician