6.1 Algoritmo de Esperanza-Maximización (EM)

Recordemos que el problema de maximización de verosimilitud se vuelve más difícil en la presencia de datos faltantes.

Ejemplo. Escogemos al azar una moneda y luego tiramos un volado con esa moneda (modelo Moneda -> Sol):

Si parametrizamos con \(\theta=P(moneda=A)\), \(\theta_A=P(sol|moneda=A)\) y \(\theta_B=P(sol|moneda=B)\), la log verosimilitud para los datos completos es \[\log (\theta\theta_A) + \log(\theta(1-\theta_A)) + 2\log((1-\theta)(1-\theta_B))+ \log((1-\theta)\theta_B)\] \[= 2\log\theta+3\log(1-\theta)+\log\theta_A+\log(1-\theta_A)+2\log(1-\theta_B)+\log\theta_B\] En este caso es fácil encontrar los estimadores de máxima verosimilitud. Ahora, para incluir la contribución del caso con faltante suponemos MAR y promediamos los valores faltantes (como indica la fórmula en (5), cambiando la integral por suma):

\[p_{\theta}(x^{6}_{sol}=0 )=p_{\theta}(x^{6}_{sol}=0 |x^{6}_{moneda}=A) p_{\theta}(x^{6}_{moneda}=A) + p_{\theta}(x^{6}_{sol}=0 |x^{6}_{moneda}=B) p_{\theta}(x^{6}_{moneda}=B),\] y ahora buscamos maximizar \[p_{\theta}(y_{obs})=2\log\theta+3\log(1-\theta)+\log\theta_A+\log(1-\theta_A)+2\log(1-\theta_B)+\log\theta_B + \log((1-\theta_A)\theta + (1-\theta_B)(1-\theta)).\]

Notemos que esta expresión es más difícil de maximizar. El algoritmo EM da un algoritmo iterativo para encontrar máximos locales de la verosimilitud \(p_{\theta} (y_{obs})\). Este algoritmo sigue una estrategia de optimización motivada por la noción misma de datos faltantes y considerando la distribución condicional de los faltantes dado lo observado.

El algoritmo de esperanza-maximización (EM) es un proceso iterativo para maximizar verosimilitud en presencia de datos faltantes. Cada iteración consiste de dos pasos:

  • Calcular el valor esperado de la log-verosimilitud promediando sobre datos faltantes con una aproximación de la solución (\(\theta^{(t)}\)).

  • Obtener una nueva aproximación maximizando la cantidad resultante en el paso previo.

Podemos ver que el algorimto EM que el paso de esperanza consiste en un soft assignment de cada valor faltante de acuerdo al modelo de los datos y los datos observados.

Denotemos por \(\theta^{(t)}\) el maximizador estimado en la iteración \(t\) (\(t=0,1,...\)). En lugar de tratar directamente con \(\log p_{\theta}(y_{obs})=\log p(y_{obs}\theta)\) el algoritmo EM trabaja sobre \(Q(\theta \vert \theta^{(t)})\) que definimos como la esperanza de la log-verosimilitud de los datos completos \(y=(y_{obs},y_{falta})\) condicional a los datos observados \(y_{obs}\):

\[Q(\theta \vert \theta^{(t)})=E\big[\log\mathcal{L}(\theta|y) \big|y_{obs}, \theta^{(t)}\big]\]

\[=E \big[\log p(y|\theta) \big |y_{obs}, \theta^{(t)}\big]\]

\[=\int_{y_{falta}} [\log p(y|\theta)] p(y_{falta}|y_{obs}, \theta^{(t)})dy_{falta}.\]

Es así que en el paso esperanza calculamos \[Q(\theta|\theta^{(t)})\] y en el paso maximización: \[\theta^{(t+1)} = argmax_{\theta}Q(\theta|\theta^{(t)})\]

Ejemplo. En nuestro ejemplo anterior, haríamos lo siguiente. Primero definimos la función que calcula el esperado de la log verosimilitud:

Ahora escogemos una solución inicial y calculamos el esperado

Optimizamos este valor esperado:

Y vemos que recuperamos la misma solución que en el caso de maximizar usando la función optim.

Veamos porque funciona trabajar con \(Q(\theta|\theta^{(t)})\). Si escribimos: \[ p(y_{obs}|\theta)=\frac{p(y|\theta)}{p(y_{falta}|y_{obs},\theta)}.\]

Tomando logaritmos, \[\log{p}(y_{obs}|\theta)=\log{p}(y|\theta)-log{p}(y_{falta}|y_{obs}, \theta).\]

Y tomando el valor esperado con respecto la distribuión de \(y_{falta}|y_{obs}\) con parámetro \(\theta^{(t)}\) (que consideramos como una aproximación de los verdaderos parámetros), obtenemos:

\[E\big[\log p(y_{obs}|\theta)\big | y_{obs}, \theta^{(t)}\big]= E\big[\log p(y|\theta)\big | y_{obs}, \theta^{(t)}\big] - E\big[\log p(y_{falta}|\theta)\big | y_{obs}, \theta^{(t)}\big]\]

\[=\int_{y_{falta}}[\log p(y|\theta)] p(y_{falta}|y_{obs}, \theta^{(t)})dy_{falta} - \int_{y_{falta}}[\log{p}(y_{falta}| y_{obs},\theta)]p(y_{falta}|y_{obs},\theta^{(t)})dy_{falta}\] \[= Q(\theta|\theta^{(t)}) - H(\theta|\theta^{(t)})\]

La igualdad anterior es cierta para cualquier valor de \(\theta\), en particular, \(\theta=\theta^{(t)}\) por lo que: \[\log{p}(y_{obs}|\theta) -\log{p}(y_{obs}|\theta^{(t)})= Q(\theta|\theta^{(t)}) - Q(\theta^{(t)}|\theta^{(t)}) - [H(\theta|\theta^{(t)}) - H(\theta^{(t)}|\theta^{(t)})]\] utilizando la desigualdad de Jensen se puede demostrar que \(H(\theta|\theta^{(t)}) \leq H(\theta^{(t)}|\theta^{(t)})\), por lo que (en términos de la verosimilutd) \[\log\mathcal{L}(\theta|y_{obs}) -\log\mathcal{L}(\theta^{(t)}|t_{obs}) \geq Q(\theta|\theta^{(t)}) - Q(\theta^{(t)}|\theta^{(t)})\] y vemos que si incrementamos \(Q(\theta|\theta^{(t)})\) entonces también incrementamos \(\mathcal{L}(\theta|y_{obs})\).

6.1.1 Observaciones del algoritmo EM

  • El algoritmo EM es una generalización natural de la estiamación por máxima verosimilitud cuando hay presencia de datos faltantes.

  • El problema de maximización que ataca EM es más complejo que el problema de máxima verosimilitud con datos completos. En el caso usual \(\log p(x|\theta)\) tiene un único máximo global que muchas veces se puede encontrar de forma cerrada, por su parte la verosimilitud con datos faltantes tiene muchos máximos locales y no tiene solución cerrada.

  • El algoritmo EM reduce el problema de múltiples máximos a una secuencia de subproblemas (\(Q(\theta|\theta^{(t)}\)) cada uno con un único máximo global.
    Estos subproblemas son tales que se garantiza convergencia de las soluciones \(\theta^{(1)}, \theta^{(2)}, ...\) a un máximo, pues a cada paso se incrementa monótonamente el valor de la log-verosimilitud de los datos observados.

  • El algoritmo EM garantiza convergencia únicamente a un máximo local pues la secuencia de \(\theta^{(t)}\) depende del valor con el que se inicializa el proceso \(\theta^{(1)}\). Por tanto, es conveniente repetir el algoritmo con varios valores iniciales.

  • Una manera de entender el algoritmo EM es pensar que en cada paso de esperanza imputas los faltantes, sin embargo, en lugar de una imputación dura se realiza un soft assignment de cada valor faltante. El soft assignment se calcula obteniendo la probabilidad de cada alternativa para el valor faltante (usando el parámetro \(\theta^{(t)})\)) para crear una base de datos ponderada que considera todos los posibles escenarios de los faltantes. Al usar los valores esperados (o datos ponderados) en lugar de realizar una imputación dura el algoritmo toma en cuenta la confianza del modelo cada vez que completa los datos.