Ich möchte den EM-Algorithmus manuell implementieren und ihn dann mit den Ergebnissen des normalmixEM
of- mixtools
Pakets 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):
Der E- Schritt ist nun die Berechnung der bedingten Erwartung:
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?
quelle
Antworten:
You have several problems in the source code:
As @Pat pointed out, you should not use log(dnorm()) as this value can easily go to infinity. You should use logmvdnorm
When you use sum, be aware to remove infinite or missing values
You looping variable k is wrong, you should update loglik[k+1] but you update loglik[k]
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).
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:
Historgramm
quelle
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.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 y Wenn 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
bewerten
Mathematisch dasselbe, aber denken Sie darüber nach, was wann passiertf( y| θ) und τ sind ≈ 0 . Derzeit erhalten Sie:
but with tau moved you get
assuming R evaluates00=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:
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
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 ) = - In f .
quelle