¿Por qué el algoritmo EM tiene que ser iterativo?

9

Suponga que tiene una población con unidades, cada una con una variable aleatoria . Usted observa valores de para cualquier unidad para la cual . Queremos una estimación de .X iPoisson ( λ ) n = N - n 0 X i > 0 λNXiPoisson(λ)n=Nn0Xi>0λ

Hay métodos de momentos y formas condicionales de máxima probabilidad de obtener la respuesta, pero quería probar el algoritmo EM. Obtengo que el algoritmo EM es donde el subíndice indica el valor de la iteración anterior del algoritmo y es constante con respecto a Los parametros. (De hecho, creo que la en la fracción entre paréntesis debería ser-1Knn+1

Q(λ1,λ)=λ(n+nexp(λ1)1)+log(λ)i=1nxi+K,
1Knn+1 , pero eso no parece exacto; una pregunta para otro momento).

Para hacer esto concreto, suponga que , . Por supuesto, y no son observados y se debe estimar .x i = 20 N n 0 λn=10xi=20Nn0λ

Cuando itero la siguiente función, conectando el valor máximo de la iteración anterior, llego a la respuesta correcta (verificada por CML, MOM y una simulación 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

Pero este es un problema simple; vamos a maximizar sin iterar:

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

El valor de la función es más alto que en el procedimiento no iterativo y el resultado es inconsistente con las otras metodologías. ¿Por qué el segundo procedimiento da una respuesta diferente y (supongo) incorrecta?

Charlie
fuente

Respuestas:

6

Cuando haya encontrado su función objetivo para el algoritmo EM, supongo que trató el número de unidades con , que llamaré , como su parámetro latente. En este caso, estoy (nuevamente) asumiendo que representa una forma reducida del valor esperado sobre de la probabilidad dada . Esto no es lo mismo que la probabilidad total, porque esoy Q y λ - 1 λ - 1xi=0yQy λ1λ1 se pisa según lo dado.

Por lo tanto, no puede usar para la probabilidad completa, ya que no contiene información sobre cómo cambiar cambia la distribución de (y desea seleccionar también los valores más probables de cuando maximiza la probabilidad completa). Esta es la razón por la cual la probabilidad máxima completa para el Poisson truncado cero difiere de su función , y por qué obtiene una respuesta diferente (e incorrecta) cuando maximiza .λ y y Q f ( λ ) = Q ( λ , λ )QλyyQf(λ)=Q(λ,λ)

Numéricamente, maximizar necesariamente dará como resultado una función objetivo al menos tan grande como su resultado EM, y probablemente más grande ya que no hay garantía de que el algoritmo EM converja a un máximo de , solo se supone que converge a ¡un máximo de la función de probabilidad !ff(λ)f

Jayk
fuente