Załóżmy, że masz populację jednostki, każda z losową zmienną . Ty obserwujeszwartości dla dowolnej jednostki, dla której . Chcemy oszacowania .
Istnieją metody chwil i warunkowe maksymalne prawdopodobieństwo uzyskania odpowiedzi, ale chciałem wypróbować algorytm EM. Otrzymuję algorytm EM: gdzie indeks dolny wskazuje wartość z poprzedniej iteracji algorytmu, a jest stały w odniesieniu do Parametry. (Właściwie uważam, że we frakcji w nawiasach powinno wynosić , ale to nie wydaje się dokładne; pytanie na inny raz).
Aby uczynić to konkretnym, załóżmy, że , . Oczywiście, i nie są obserwowane i należy oszacować .
Gdy iteruję następującą funkcję, podając maksymalną wartość z poprzedniej iteracji, docieram do prawidłowej odpowiedzi (zweryfikowanej przez CML, MOM i prostą symulację):
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
Ale to jest prosty problem; zmaksymalizujmy bez iteracji:
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
Wartość funkcji jest wyższa niż w procedurze iteracyjnej, a wynik jest niezgodny z innymi metodologiami. Dlaczego druga procedura daje inną i (jak sądzę) nieprawidłową odpowiedź?