Pourquoi l'algorithme EM doit-il être itératif?

9

Supposons que vous ayez une population avec unités, chacune avec une variable aléatoire . Vous observez valeurs pour toute unité pour laquelle . Nous voulons une estimation de .Nn = N - n 0 X i > 0 λXiPoisson(λ)n=Nn0Xi>0λ

Il existe une méthode des moments et des moyens conditionnels de probabilité maximale d'obtenir la réponse, mais je voulais essayer l'algorithme EM. J'obtiens l'algorithme EM comme où l' indice indique la valeur de l'itération précédente de l'algorithme et est constant par rapport à Les paramètres. (Je pense en fait que le dans la fraction entre parenthèses devrait être , mais cela ne semble pas exact; une question pour une autre fois).-1Knn+1

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

Pour rendre cela concret, supposons que , . Bien sûr, et sont pas observés et doit être estimé.x i = 20 N n 0n=10xi=20Nn0λ

Lorsque j'itère la fonction suivante, en branchant la valeur maximale de l'itération précédente, j'atteins la bonne réponse (vérifiée par CML, MOM et une simulation simple):

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

Mais c'est un problème simple; maximisons sans itérer:

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

La valeur de la fonction est plus élevée que dans la procédure non itérative et le résultat est incompatible avec les autres méthodologies. Pourquoi la deuxième procédure donne-t-elle une réponse différente et (je présume) incorrecte?

Charlie
la source

Réponses:

6

Lorsque vous avez trouvé votre fonction objective pour l'algorithme EM, je suppose que vous avez traité le nombre d'unités avec , que j'appellerai , comme paramètre latent. Dans ce cas, je suppose (encore) que représente une forme réduite de la valeur attendue sur de la probabilité donnée . Ce n'est pas la même chose que la vraisemblance complète, parce que est traité comme indiqué.y Q y λ - 1 λ - 1xi=0yQy λ1λ1

Par conséquent, vous ne pouvez pas utiliser pour la pleine probabilité, car il ne contient pas d'informations sur la façon dont la modification de modifie la distribution de (et vous souhaitez également sélectionner les valeurs les plus probables de lorsque vous maximisez la pleine probabilité). C'est pourquoi la probabilité maximale totale pour le Poisson tronqué zéro diffère de votre fonction , et pourquoi vous obtenez une réponse différente (et incorrecte) lorsque vous maximisez .λ y y Q f ( λ ) = Q ( λ , λ )QλyyQf(λ)=Q(λ,λ)

Numériquement, maximiser entraînera nécessairement une fonction objective au moins aussi grande que votre résultat EM, et probablement plus grande car il n'y a aucune garantie que l'algorithme EM convergera vers un maximum de - il est seulement supposé converger vers un maximum de la fonction de vraisemblance !ff(λ)f

jayk
la source