Warum muss der EM-Algorithmus iterativ sein?

9

Angenommen, Sie haben eine Population mit Einheiten mit jeweils einer Zufallsvariablen . Sie beobachten Werte für jede Einheit, für die . Wir wollen eine Schätzung von .NXiPoisson(λ)n=Nn0Xi>0λ

Es gibt Methoden für Momente und bedingte Möglichkeiten mit maximaler Wahrscheinlichkeit, um die Antwort zu erhalten, aber ich wollte den EM-Algorithmus ausprobieren. Ich erhalte den EM-Algorithmus als wobei der Index den Wert aus der vorherigen Iteration des Algorithmus angibt und in Bezug auf konstant ist Die Parameter. (Ich denke tatsächlich, dass das in dem Bruch in Klammern , aber das scheint nicht genau zu sein; eine Frage für ein anderes Mal).

Q(λ1,λ)=λ(n+nexp(λ1)1)+log(λ)i=1nxi+K,
1Knn+1

Um dies konkret zu machen, nehmen wir an, dass , . Natürlich werden und nicht beobachtet und ist zu schätzen.n=10xi=20Nn0λ

Wenn ich die folgende Funktion iteriere und den Maximalwert der vorherigen Iteration einfüge, erreiche ich die richtige Antwort (überprüft durch CML, MOM und eine einfache Simulation):

EmFunc <- function(lambda, lambda0){
  -lambda * (10 + 10 / (exp(lambda0) - 1)) + 20 * log(lambda)
}

lambda0 <- 2
lambda  <- 1

while(abs(lambda - lambda0) > 0.0001){
  lambda0 <- lambda
  iter    <- optimize(EmFunc, lambda0 = lambda0, c(0,4), maximum = TRUE)
  lambda  <- iter$maximum
}

> iter
$maximum
[1] 1.593573

$objective
[1] -10.68045

Dies ist jedoch ein einfaches Problem. Lassen Sie uns einfach maximieren, ohne zu iterieren:

MaxFunc <- function(lambda){
  -lambda * (10 + 10 / (exp(lambda) - 1)) + 20 * log(lambda)
}

optimize(MaxFunc, c(0,4), maximum = TRUE)
$maximum
[1] 2.393027

$objective
[1] -8.884968

Der Wert der Funktion ist höher als bei der nicht iterativen Prozedur und das Ergebnis stimmt nicht mit den anderen Methoden überein. Warum gibt das zweite Verfahren eine andere und (ich nehme an) falsche Antwort?

Charlie
quelle

Antworten:

6

Wenn Sie Ihre Zielfunktion für den EM-Algorithmus gefunden haben, haben Sie vermutlich die Anzahl der Einheiten mit , die ich als latenten Parameter . In diesem Fall gehe ich (wieder) davon aus, dass eine reduzierte Form des erwarteten Wertes über der mit gegebenen Wahrscheinlichkeit darstellt . Dies ist nicht das gleiche wie die volle Wahrscheinlichkeit, weil das treadted wird als gegeben.xi=0yQy λ1λ1

Daher können Sie für die volle Wahrscheinlichkeit verwenden, da es keine Informationen darüber enthält, wie das Ändern von die Verteilung von ändert (und Sie möchten auch die wahrscheinlichsten Werte von auswählen, wenn Sie die volle Wahrscheinlichkeit maximieren). Aus diesem Grund unterscheidet sich die volle maximale Wahrscheinlichkeit für das auf Null abgeschnittene Poisson von Ihrer Funktion, und Sie erhalten eine andere (und falsche) Antwort, wenn Sie maximieren .QλyyQf(λ)=Q(λ,λ)

Numerisch gesehen führt das Maximieren von zwangsläufig zu einer Zielfunktion, die mindestens so groß ist wie Ihr EM-Ergebnis, und wahrscheinlich größer, da es keine Garantie dafür gibt, dass der EM-Algorithmus zu einem Maximum von konvergiert - es soll nur zu konvergieren ein Maximum der Wahrscheinlichkeitsfunktion !f(λ)f

Jayk
quelle