EM-Algorithmus manuell implementiert

20

Ich möchte den EM-Algorithmus manuell implementieren und ihn dann mit den Ergebnissen des normalmixEMof- mixtoolsPakets vergleichen. Natürlich würde ich mich freuen, wenn beide zu den gleichen Ergebnissen führen würden. Die Hauptreferenz ist Geoffrey McLachlan (2000), Finite Mixture Models .

Ich habe eine Mischungsdichte von zwei Gaußschen, in allgemeiner Form ist die log-Wahrscheinlichkeit gegeben durch (McLachlan Seite 48):

logLc(Ψ)=i=1gj=1nzij{logπi+logfi(yi;θi)}.
Diezij sind1 , wenn die Beobachtung aus deri tenKomponentendichtestammte, sonst0 . Dasfi ist die Dichte der Normalverteilung. Dasπ ist das Mischungsverhältnis, also istπ1 die Wahrscheinlichkeit, dass eine Beobachtung aus der ersten Gaußschen Verteilung stammt, undπ2 ist die Wahrscheinlichkeit, dass eine Beobachtung aus der zweiten Gaußschen Verteilung stammt.

Der E- Schritt ist nun die Berechnung der bedingten Erwartung:

Q(Ψ;Ψ(0))=EΨ(0){logLc(|Ψ)|y}.
was nach ein paar Ableitungen zum Ergebnis führt (Seite 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))
im Fall von zwei Gaußschen (Seite 82):

τi(yj;Ψ)=πiϕ(yj;μi,Σi)h=1gπhϕ(yj;μh,Σh)
DerM-Schritt ist jetzt die Maximierung von Q (Seite 49):

Q(Ψ;Ψ(k))=i=1gj=1nτi(yj;Ψ(k)){logπi+logfi(yj;θi)}.
Dies führt zu (im Fall von zwei Gaußschen) (Seite 82):

μich(k+1)=j=1nτichj(k)yjj=1nτichj(k)Σich(k+1)=j=1nτichj(k)(yj-μich(k+1))(yj-μich(k+1))Tj=1nτichj(k)
and we know that (p. 50)

πich(k+1)=j=1nτich(yj;Ψ(k))n(ich=1,,G).
We repeat the E, M steps until L(Ψ(k+1))-L(Ψ(k)) is small.

I tried to write a R code (data can be found here).

# 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

The algorithm is not working, since some observations have the likelihood of zero and the log of this is -Inf. Where is my mistake?

Stat Tistician
quelle
The problem is not a statistical one, but rather a numerical one. You should add contingencies for likelihoods smaller than machine precision in your code.
JohnRos
why dont you try veryfying the mixtools function with a very simple example that can be verified by hand , say just five or ten values and two timeseries,first. then, if you find it works there, generalize your code and verify at each step.

Antworten:

17

You have several problems in the source code:

  1. As @Pat pointed out, you should not use log(dnorm()) as this value can easily go to infinity. You should use logmvdnorm

  2. When you use sum, be aware to remove infinite or missing values

  3. You looping variable k is wrong, you should update loglik[k+1] but you update loglik[k]

  4. The initial values for your method and mixtools are different. You are using Σ in your method, but using σ for mixtools(i.e. standard deviation, from mixtools manual).

  5. Your data do not look like a mixture of normal (check histogram I plotted at the end). And one component of the mixture has very small s.d., so I arbitrarily added a line to set τ1 and τ2 to be equal for some extreme samples. I add them just to make sure the code can work.

Ich schlage auch vor, dass Sie vollständige Codes (z. B. wie Sie loglik [] initialisieren) in Ihren Quellcode einfügen und den Code einrücken, um das Lesen zu vereinfachen.

Immerhin vielen Dank, dass Sie das mixtools- Paket eingeführt haben und ich plane, es für meine zukünftige Forschung zu verwenden.

Ich habe auch meinen Arbeitscode als Referenz angegeben:

# 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

Historgramm Histogram

zhanxw
quelle
@zahnxw danke für deine antwort, heißt das also, dass mein code falsch ist? Die Grundidee funktioniert also nicht?
Statistiker
"Ich schlage auch vor, dass Sie vollständige Codes (z. B. wie Sie loglik [] initialisieren) in Ihren Quellcode einfügen und den Code einrücken, um das Lesen zu vereinfachen." Na das ist mein Code? Das Loglik [] ist so definiert, wie ich es in dem von mir geposteten Code deklariert habe.
Statistiker
1
@StatTistician die Idee ist richtig, aber die Umsetzung hat Mängel. Beispielsweise haben Sie einen Unterlauf nicht berücksichtigt. Außerdem ist die Schleifenvariable k verwirrend. Sie setzen zuerst loglik [1] und loglik [2], nachdem Sie die while-Schleife eingegeben haben, setzen Sie loglik [1] erneut. Dies ist nicht der natürliche Weg. Mein Vorschlag zum Initialisieren von loglik [] bedeutet code loklik <- rep(NA, 100):, der loglik [1], loglik [2] ... loglik [100] vorbelegt. Ich stelle diese Frage, weil ich in Ihrem ursprünglichen Code die Delkaration von loglik nicht gefunden habe. Vielleicht wird der Code beim Einfügen abgeschnitten.
Zhanxw
Wie ich unten gepostet habe: Danke für Ihre Hilfe, aber ich lösche dieses Thema, da es für mich zu fortgeschritten ist.
Statistiker
Gibt es jetzt eine Möglichkeit zu bestimmen, welcher Teil der Daten zu welcher Mischung gehört?
Kardinal
2

Beim Versuch, Ihre .rar-Datei zu öffnen, wird immer wieder eine Fehlermeldung angezeigt, aber möglicherweise tue ich nur etwas Dummes.

Ich kann keine offensichtlichen Fehler in Ihrem Code sehen. Ein möglicher Grund, warum Sie Nullen erhalten, liegt in der Gleitkommapräzision. Denken Sie daran, wenn Sie rechnenf(y;θ), Sie bewerten exp(-0,5(y-μ)2/σ2). Es macht keinen großen Unterschied zwischenμ und yWenn Sie dies auf einem Computer tun, wird dies auf 0 abgerundet. Dies macht sich in Mischungsmodellen doppelt bemerkbar, da einige Ihrer Daten nicht jeder Mischungskomponente "zugeordnet" werden und daher sehr weit davon entfernt sein können. Theoretisch sollten diese Punkte auch einen niedrigen Wert von habenτ Wenn Sie die Log-Wahrscheinlichkeit auswerten, um dem Problem entgegenzuwirken, wurde die Menge dank des Gleitkomma-Fehlers zu diesem Zeitpunkt bereits als -Inf ausgewertet, sodass alles kaputt geht :).

Wenn dies das Problem ist, gibt es einige mögliche Lösungen:

Eine ist, deine zu bewegen τinnerhalb des Logarithmus. Also anstatt zu bewerten

τLog(f(y|θ))

bewerten

log(f(y|θ)τ).

Mathematisch dasselbe, aber denken Sie darüber nach, was wann passiert f(y|θ) und τ sind 0. Derzeit erhalten Sie:

  • 0log(0)=0(Inf)=NaN

but with tau moved you get

  • log(00)=log(1)=0

assuming R evaluates 00=1 (I don't know if it does or not as I tend to use matlab)

Another solution is to expand out the stuff inside the logarithm. Assuming you're using natural logarithms:

τlog(f(y|θ))

=τlog(exp(0.5(yμ)2/σ2)/2πσ2)

=0.5τlog(2πσ2)0.5τ(yμ)2σ2.

Mathematisch dasselbe, sollte aber gegenüber Gleitkommafehlern widerstandsfähiger sein, da Sie die Berechnung einer großen negativen Potenz vermieden haben. Dies bedeutet, dass Sie die eingebaute Normauswertungsfunktion nicht mehr verwenden können. Wenn dies jedoch kein Problem darstellt, ist dies wahrscheinlich die bessere Antwort. Nehmen wir zum Beispiel an, wir haben die Situation, in der

-0,5(y-μ)2σ2=-0,5402=-800.

Bewerten Sie das, wie ich es vorgeschlagen habe, und Sie erhalten -800. In Matlab erhalten wir jedoch, wenn wir das Protokoll herausnehmenLog(exp(-800))=Log(0)=-ichnf.

Klopfen
quelle
mh, um ehrlich zu sein: Ich bin nicht gut genug, um dieses Ding zum Laufen zu bringen. Was mich interessiert hat, ist: Kann ich mit meinem Algorithmus das gleiche Ergebnis erzielen wie mit der implementierten Version des mixtools-Pakets? Aber aus meiner Sicht scheint dies nach dem Mond zu fragen. Aber ich denke, Sie geben sich Mühe, und ich werde es akzeptieren! Vielen Dank!
Statistiker