Los estimadores de máxima verosimilitud tienen excelentes propiedades, pero a veces involucran un problema de optimización difícil de resolver. El algoritmo EM (Expectation-Maximization) es un método para calcular el estimador de máxima verosimilitud indicado en casos en los que falta cierta información tal que el problema sería mucho más sencillo si se dispusiera de ella. Esto ocurre, por ejemplo, cuando hay datos faltantes, cuando los datos están censurados o cuando corresponden a una mezcla de distribuciones, sin que se sepa a qué distribución corresponde cada observación.
El algoritmo fue propuesto por Dempster, Laird y Rubin (1977) en un famoso artículo que acumula más de 70.000 citas en Google Scholar. En 2021, Nan Laird recibió el Premio Internacional de Estadística. En este enlace se puede encontrar un artículo de Altea Lorenzo muy recomendable sobre su vida y contribuciones.
Explico a continuación brevemente la idea del algoritmo y lo ilustro con el ejemplo más sencillo que se me ha ocurrido.
Disponemos de una muestra \(x=(x_1,\ldots,x_n)\) procedente de una distribución que depende de un parámetro \(\theta\) y queremos calcular el estimador de máxima verosimilitud de \(\theta\).
Llamamos \(y=(y_1,\ldots,y_m)\) a otros datos correspondientes a información adicional que convertiría el problema en otro mucho más fácil. Normalmente estos datos \(y\) son en realidad inaccesibles, pero eso no importa en la aplicación del método.
En esta situación podemos plantearnos dos problemas:
Aunque \(y\) sea en realidad inobservable, si maximizar \(\ell_c\) es más fácil que maximizar \(\ell\) la siguiente idea puede ser útil:
Idea: Sustituir la verosimilitud completa por una aproximación a partir de la información que realmente tenemos. Para ello, se aproxima \(\ell_c(\theta;X,Y)\) calculando su esperanza condicionada dado \(X=x\). Como hay que elegir un valor del parámetro, se usa el más reciente obtenido en la aplicación del algoritmo.
La aplicación de la idea anterior lleva al algoritmo EM:
Iniciamos el método con algún valor sensato \(\theta_0\). Sea \(\theta_k\) el valor del parámetro tras \(k\) iteraciones. En cada iteración hay que dar los dos pasos siguientes:
\[Q(\theta;\theta_k,x):= \mbox{E}_{Y|x;\theta_k}[\ell_c(\theta;x,Y)].\]
Se reiteran los pasos E y M hasta que se cumple algún criterio de convergencia.
Bajo condiciones de regularidad (hay muchos artículos que analizan el comportamiento de este método) puede demostrarse que en cada paso se mejora el objetivo del problema que realmente queremos resolver, es decir, \(\ell(\theta_{k+1}; x)\geq \ell(\theta_k; x)\) y que la sucesión de valores generada converge a \(\hat\theta\), el estimador de máxima verosimilitud correspondiente a la muestra \(x\) que realmente tenemos.
Este es un ejemplo sencillo que sirve para ilustrar las ideas anteriores. Para resolver el mismo problema podría usarse cualquier otro método de optimización razonable.
Supongamos que las observaciones \(x_1,\ldots,x_n\) son independientes y vienen con probabilidad \(1-p\) de una distribución con densidad \(f_0\) y con probabilidad \(p\) de otra distribución \(f_1\). Tanto \(f_0\) como \(f_1\) son conocidas y el objetivo es estimar la proporción de observaciones \(p\) que vienen de \(f_1\). La función de densidad de los datos es \(f(x) = (1-p)f_0(x) + pf_1(x)\) (para comprobarlo, se usa la fórmula de la probabilidad total). Por lo tanto, la log-verosimilitud que tenemos que maximizar es \[\ell(p; x) = \sum_{i=1}^n \log [(1-p)f_0(x_i) + p f_1(x_i)].\] No es posible encontrar explícitamente el valor de \(p\) que maximiza esta función por lo que hay que recurrir a algún método numérico.
Si alguien nos dijera a qué distribución pertenece cada observación el problema sería inmediato de resolver. Escribimos \(y_i=1\) si la observación \(x_i\) viene de \(f_1\) e \(y_i=0\) si viene de \(f_0\). La verosimilitud completa correspondiente a \((x,y)\) es \[ L_c(p;x,y) := \prod_{i=1}^n [(1-p)f_0(x_i)]^{1-y_i}[p f_1(x_i)]^{y_i}, \] y tomando logaritmos, \[ \ell_c(p; x,y) = \sum_{i=1}^n [(1-y_i)\log(1-p)+y_i\log p ] + C, \] donde \(C\) no depende del parámetro \(p\), por lo que es irrelevante para nuestros propósitos. Derivando respecto a \(p\) es fácil ver que el óptimo se alcanza en el valor \(n^{-1}\sum_{i=1}^n y_i\). Como no podía ser de otra manera, si conocemos de dónde viene cada observación estimamos \(p\) como la proporción de valores muestrales que vienen de \(f_1\).
Nosotros no sabemos de qué distribución viene cada \(x_i\) pero si aplicamos la fórmula de Bayes y usamos el valor \(p_k\) del parámetro,
\[ y_{i,k} :=\mbox{E}(Y_i|x_i;p_k) = \frac{p_kf_1(x_i)}{p_kf_1(x_i) + (1-p_k)f_0(x_i)}. \tag{1}\]
Como consecuencia, los dos pasos del algoritmo se reducen a:
PASO E:
Salvo sumandos que no dependen de \(p\),
\[Q(p;p_k) = \sum_{i=1}^n [(1-y_{i,k})\log(1-p) + y_{i,k}\log p],\] donde \(y_{i,k}\) se calcula mediante la ecuación (1).
Paso M:
Maximizamos \(Q(p;p_k)\). Para ello derivamos e igualamos a cero y obtenemos \[ p_{k+1} = \frac{1}{n}\sum_{i=1}^n y_{i,k}. \] Un ejemplo numérico con datos simulados se muestra a continuación. Primero simulamos \(n=100\) datos que pueden venir de una normal estándar (con probabilidad \(p=0.8\)), o de una normal de media 4 y varianza 1 (con probabilidad \(1-p=0.2\)):
# Generamos datos de una N(0,1) o de una N(4,1)
set.seed(156)
n <- 100
p <- 0.8 # objetivo a estimar
media <- 4
y <- rbinom(n, 1, p)
x <- rnorm(n)
n0 <- sum(y==0)
x[y==0] <- rnorm(n0, mean = media)
# Representa histograma de los datos
library(ggplot2)
ggplot(data.frame(x)) +
geom_histogram(aes(x = x), col = 'black', fill = 'olivedrab4', bins = 20)
En el histograma se observa claramente la bimodalidad correspondiente a las dos posibles distribuciones de las que pueden venir los datos. Ahora fingimos que no conocemos que \(p=0.8\) y aplicamos el algoritmo EM a la muestra x (no usamos para nada y, que es la variable que nos informa de qué distribución viene cada observación). Usamos el valor inicial \(p_0=1/2\).
# Iniciamos el algoritmo
p_k <- 0.5
tol <- 0.1
# Pasos E-M
while(tol > 0.001){
p_0 <- p_k
y_k <- p_0*dnorm(x)/(p_k*dnorm(x) + (1-p_0)*dnorm(x, mean = media)) # Paso E
p_k <- mean(y_k) # Paso M
tol <- abs(p_k-p_0)
}
p_k
#> [1] 0.8179144Se obtiene el valor 0.818, cercano al valor verdadero. Hemos supuesto las dos distribuciones conocidas. Una situación de más interés y más complicada sería no suponer las medias conocidas y tratar de estimar por máxima verosimilitud tanto \(p\) como las dos medias.