Por que o algoritmo EM precisa ser iterativo?

9

Suponha que você tenha uma população com unidades, cada uma com uma variável aleatória . Você observa valores para qualquer unidade para a qual . Queremos uma estimativa de .X i isson Poisson ( λ ) n = N - n 0 X i > 0 λNXiPoisson(λ)n=Nn0Xi>0λ

Existem métodos de momentos e formas condicionais de probabilidade máxima de obter a resposta, mas eu queria experimentar o algoritmo EM. Eu recebo o algoritmo EM para que o índice indica o valor da iteração anterior do algoritmo e é constante em relação a os parametros. (Na verdade, acho que o na fração entre parênteses deve ser , mas isso não parece preciso; uma pergunta para outra hora).-1Knn+1

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

Para tornar isso concreto, suponha que , . Obviamente, e não são observados e deve ser estimado.x i = 20 N n 0 λn=10xi=20Nn0λ

Quando eu itero a função a seguir, inserindo o valor máximo da iteração anterior, alcanço a resposta correta (verificada por CML, MOM e uma simulação simples):

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

Mas este é um problema simples; vamos maximizar sem 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

O valor da função é maior que no procedimento não iterativo e o resultado é inconsistente com as outras metodologias. Por que o segundo procedimento está dando uma resposta diferente e (presumo) incorreta?

Charlie
fonte

Respostas:

6

Quando você encontra sua função objetivo para o algoritmo EM, presumo que você tenha tratado o número de unidades com , que chamarei de , como seu parâmetro latente. Nesse caso, estou (novamente) assumindo que representa uma forma reduzida do valor esperado sobre da probabilidade fornecida . Isso não é o mesmo que a probabilidade total, porque esse é trocado como determinado.y Q y λ - 1 λ - 1xi=0yQy λ1λ1

Portanto, você não pode usar para a probabilidade total, pois isso não contém informações sobre como a alteração altera a distribuição de (e você deseja selecionar os valores mais prováveis ​​de também ao maximizar a probabilidade total). É por isso que a probabilidade máxima máxima para o Poisson truncado zero difere da sua função e por que você recebe uma resposta diferente (e incorreta) quando maximiza .λ y y Q f ( λ ) = Q ( λ , λ )QλyyQf(λ)=Q(λ,λ)

Numericamente, maximizar necessariamente resultará em uma função objetiva pelo menos tão grande quanto o resultado EM, e provavelmente maior, pois não há garantia de que o algoritmo EM convergirá para um máximo de - é suposto convergir apenas para um máximo da função de probabilidade !ff(λ)f

jayk
fonte