Algoritmo de expectativa-maximización

En estadística, un algoritmo de expectativa-maximización (EM) es un método iterativo para encontrar la máxima verosimilitud (local) o el máximo a posteriori. (MAP) estimaciones de parámetros en modelos estadísticos, donde el modelo depende de variables latentes no observadas. La iteración EM alterna entre realizar un paso de expectativa (E), que crea una función para la expectativa de la probabilidad logarítmica evaluada utilizando la estimación actual de los parámetros, y un paso de maximización (M), que calcula los parámetros que maximizan la probabilidad logarítmica esperada. probabilidad encontrada en el paso E. Estas estimaciones de parámetros se utilizan luego para determinar la distribución de las variables latentes en el siguiente paso E.

Historia

El algoritmo EM fue explicado y recibió su nombre en un artículo clásico de 1977 escrito por Arthur Dempster, Nan Laird y Donald Rubin. Señalaron que el método había sido "propuesto muchas veces en circunstancias especiales" por autores anteriores. Uno de los primeros es el método de recuento de genes para estimar las frecuencias alélicas de Cedric Smith. Otro fue propuesto por H.O. Hartley en 1958, y Hartley y Hocking en 1977, de donde se originaron muchas de las ideas del artículo de Dempster-Laird-Rubin. Otro de S.K Ng, Thriyambakam Krishnan y G.J McLachlan en 1977. Las ideas de Hartley pueden ampliarse a cualquier distribución discreta agrupada. Rolf Sundberg publicó un tratamiento muy detallado del método EM para familias exponenciales en su tesis y en varios artículos, tras su colaboración con Per Martin-Löf y Anders Martin-Löf. El artículo de Dempster-Laird-Rubin de 1977 generalizó el método y esbozó un análisis de convergencia para una clase más amplia de problemas. El artículo de Dempster-Laird-Rubin estableció el método EM como una importante herramienta de análisis estadístico. Véase también Meng y van Dyk (1997).

El análisis de convergencia del algoritmo Dempster-Laird-Rubin era defectuoso y C. F. Jeff Wu publicó un análisis de convergencia correcto en 1983. La prueba de Wu estableció la convergencia del método EM también fuera de la familia exponencial, como afirma Dempster-Laird-Rubin.

Introducción

El algoritmo EM se utiliza para encontrar parámetros de máxima verosimilitud (locales) de un modelo estadístico en los casos en que las ecuaciones no se pueden resolver directamente. Normalmente, estos modelos involucran variables latentes además de parámetros desconocidos y observaciones de datos conocidos. Es decir, existen valores faltantes entre los datos o el modelo puede formularse de manera más simple suponiendo la existencia de más puntos de datos no observados. Por ejemplo, un modelo de mezcla se puede describir de manera más simple suponiendo que cada punto de datos observado tiene un punto de datos no observado correspondiente, o variable latente, que especifica el componente de la mezcla al que pertenece cada punto de datos.

Encontrar una solución de máxima verosimilitud generalmente requiere tomar las derivadas de la función de verosimilitud con respecto a todos los valores desconocidos, los parámetros y las variables latentes, y resolver simultáneamente las ecuaciones resultantes. En modelos estadísticos con variables latentes, esto suele ser imposible. En cambio, el resultado suele ser un conjunto de ecuaciones entrelazadas en las que la solución de los parámetros requiere los valores de las variables latentes y viceversa, pero al sustituir un conjunto de ecuaciones por otro se produce una ecuación sin solución.

El algoritmo EM parte de la observación de que existe una manera de resolver numéricamente estos dos conjuntos de ecuaciones. Uno puede simplemente elegir valores arbitrarios para uno de los dos conjuntos de incógnitas, usarlos para estimar el segundo conjunto, luego usar estos nuevos valores para encontrar una mejor estimación del primer conjunto y luego seguir alternando entre los dos hasta que los valores resultantes sean ambos. convergen a puntos fijos. No es obvio que esto vaya a funcionar, pero se puede demostrar en este contexto. Además, se puede demostrar que la derivada de la probabilidad es (arbitrariamente cercana a) cero en ese punto, lo que a su vez significa que el punto es un máximo local o un punto de silla. En general, pueden ocurrir múltiples máximos, sin garantía de que se encuentre el máximo global. Algunas probabilidades también tienen singularidades, es decir, máximos sin sentido. Por ejemplo, una de las soluciones que puede encontrar EM en un modelo mixto implica establecer que uno de los componentes tenga varianza cero y que el parámetro medio para el mismo componente sea igual a uno de los datos. puntos.

Descripción

Los símbolos

Dado el modelo estadístico que genera un conjunto X{displaystyle mathbf {X} de datos observados, un conjunto de datos latentes no conservados o valores perdidos Z{displaystyle mathbf {Z}, y un vector de parámetros desconocidos Silencio Silencio {displaystyle {boldsymbol {theta }, junto con una función de probabilidad L()Silencio Silencio ;X,Z)=p()X,Z▪ ▪ Silencio Silencio ){displaystyle L({boldsymbol {theta };mathbf {X}mathbf {Z})=p(mathbf {X}mathbf {Z} mid {boldsymbol {theta }}}}}}}, la estimación de probabilidad máxima (MLE) de los parámetros desconocidos se determina maximizando la probabilidad marginal de los datos observados

L()Silencio Silencio ;X)=p()X▪ ▪ Silencio Silencio )=∫ ∫ p()X,Z▪ ▪ Silencio Silencio )dZ=∫ ∫ p()X▪ ▪ Z,Silencio Silencio )p()Z▪ ▪ Silencio Silencio )dZ{fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif}

Sin embargo, esta cantidad es a menudo intratable desde Z{displaystyle mathbf {Z} no se conserva y la distribución de Z{displaystyle mathbf {Z} se desconoce antes de alcanzar Silencio Silencio {displaystyle {boldsymbol {theta }.

El algoritmo EM

El algoritmo EM busca encontrar el MLE de la probabilidad marginal aplicando iterativamente estos dos pasos:

Paso de espera (E paso): Definir Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t)){displaystyle Q({boldsymbol {theta }midboldsymbol {theta }}{(t)}}}} como el valor esperado de la función de probabilidad de registro Silencio Silencio {displaystyle {boldsymbol {theta }, con respecto a la distribución condicional actual Z{displaystyle mathbf {Z} dado X{displaystyle mathbf {X} y las estimaciones actuales de los parámetros Silencio Silencio ()t){displaystyle {boldsymbol {theta} {} {}}}}} {fnMicrosoft Sans Serif}:

Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))=EZ♪ ♪ p()⋅ ⋅ SilencioX,Silencio Silencio ()t))⁡ ⁡ [log⁡ ⁡ p()X,ZSilencioSilencio Silencio )]{displaystyle Q({boldsymbol {theta }midboldsymbol {theta }}{(t)}}=operatorname {E} _{mathbf {Z} sim p(cdot)}cdot [mathbf {X}{boldsymbol {theta }} {cdot)}left[log p(mathbf {X}mathbf {Z} Silencio {boldsymbol {theta})derecha],}

Paso de maximización (M paso): Encontrar los parámetros que maximizan esta cantidad:

Silencio Silencio ()t+1)=argmaxSilencio Silencio Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t)){displaystyle {boldsymbol {theta}{(t+1)}={underset {fnMicrosoft {fnMicrosoft} - ¿Qué? ¿Qué?

De manera más sucinta, podemos escribirlo como una ecuación:

Silencio Silencio ()t+1)=argmaxSilencio Silencio EZ♪ ♪ p()⋅ ⋅ SilencioX,Silencio Silencio ()t))⁡ ⁡ [log⁡ ⁡ p()X,ZSilencioSilencio Silencio )]{displaystyle {boldsymbol {theta}{(t+1)}={underset {fnMicrosoft {fnMicrosoft} ### ## ## ## ## ## ## ##operatorname {arg,max} ## {E} _{mathbf {Z} sim p(cdot)}cdot [mathbf {X}{boldsymbol {theta }} {cdot)}left[log p(mathbf {X}mathbf {Z} Silencio {boldsymbol {theta})derecha],}

Interpretación de las variables

Los modelos típicos a los que se aplica EM Z{displaystyle mathbf {Z} como variable latente que indica la pertenencia a uno de los grupos:

1. Los puntos de datos observados X{displaystyle mathbf {X} puede ser discreto (tomar valores en un conjunto finito o contablemente infinito) o continuo (tomar valores en un conjunto incontablemente infinito). Asociado con cada punto de datos puede ser un vector de observaciones.
2. Los valores perdidos (variables latentes) Z{displaystyle mathbf {Z} son discretas, extraídas de un número fijo de valores, y con una variable latente por unidad observada.
3. Los parámetros son continuos, y son de dos tipos: Parámetros asociados a todos los puntos de datos, y aquellos asociados con un valor específico de una variable latente (es decir, asociados a todos los puntos de datos cuya variable latente correspondiente tiene ese valor).

Sin embargo, es posible aplicar EM a otros tipos de modelos.

La motivación es la siguiente. Si el valor de los parámetros Silencio Silencio {displaystyle {boldsymbol {theta } es conocido, generalmente el valor de las variables latente Z{displaystyle mathbf {Z} se puede encontrar maximizando la probabilidad de registro sobre todos los valores posibles de Z{displaystyle mathbf {Z}, o simplemente por iterating encima Z{displaystyle mathbf {Z} o a través de un algoritmo como el algoritmo Viterbi para los modelos de Markov ocultos. Por el contrario, si conocemos el valor de las variables latentes Z{displaystyle mathbf {Z}, podemos encontrar una estimación de los parámetros Silencio Silencio {displaystyle {boldsymbol {theta } bastante fácilmente, normalmente agrupando los puntos de datos observados según el valor de la variable latente asociada y promediando los valores, o alguna función de los valores, de los puntos en cada grupo. Esto sugiere un algoritmo iterativo, en el caso en que ambos Silencio Silencio {displaystyle {boldsymbol {theta } y Z{displaystyle mathbf {Z} are unknown:

1. Primero, inicialice los parámetros Silencio Silencio {displaystyle {boldsymbol {theta } a algunos valores aleatorios.
2. Computar la probabilidad de cada valor posible Z{displaystyle mathbf {Z} dado Silencio Silencio {displaystyle {boldsymbol {theta }.
3. A continuación, utilice los valores justos de Z{displaystyle mathbf {Z} para calcular una mejor estimación de los parámetros Silencio Silencio {displaystyle {boldsymbol {theta }.
4. Itear los pasos 2 y 3 hasta la convergencia.

El algoritmo que acabamos de describir se aproxima monótonamente a un mínimo local de la función de costo.

Propiedades

Aunque una iteración EM aumenta los datos observados (es decir, marginal) función de probabilidad, no existe garantía de que la secuencia converge a un estimador de probabilidad máxima. Para las distribuciones multimodales, esto significa que un algoritmo EM puede converger a un máximo local de la función de probabilidad de datos observada, dependiendo de los valores iniciales. Existen una variedad de enfoques heurísticos o metaheurísticos para escapar de un máximo local, como la subida de colinas al azar (a partir de varias estimaciones iniciales al azar Silencio Silencio ()t){displaystyle {boldsymbol {theta} {} {}}}}} {fnMicrosoft Sans Serif}), o la aplicación de métodos de anealing simulados.

EM es especialmente útil cuando la probabilidad es una familia exponencial; consulte Sundberg (2019, capítulo 8) para un tratamiento integral: el paso E se convierte en la suma de expectativas de estadísticas suficientes y el paso M implica maximizar una función lineal.. En tal caso, normalmente es posible derivar actualizaciones de expresiones de forma cerrada para cada paso, utilizando la fórmula de Sundberg (probada y publicada por Rolf Sundberg, basada en resultados no publicados de Per Martin-Löf y Anders Martin-Löf).

El método EM fue modificado para calcular estimaciones máximas a posteriori (MAP) para la inferencia bayesiana en el artículo original de Dempster, Laird y Rubin.

Existen otros métodos para encontrar estimaciones de máxima verosimilitud, como el descenso de gradiente, el gradiente conjugado o variantes del algoritmo de Gauss-Newton. A diferencia de los EM, estos métodos normalmente requieren la evaluación de la primera y/o segunda derivada de la función de probabilidad.

Prueba de corrección

Expectativa-Maximización trabaja para mejorar Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t)){displaystyle Q({boldsymbol {theta }midboldsymbol {theta }}{(t)}}}} en lugar de mejorar directamente log⁡ ⁡ p()X▪ ▪ Silencio Silencio ){displaystyle log p(mathbf {X} mid {boldsymbol {theta}}}. Aquí se muestra que las mejoras en la primera implican mejoras en esta última.

Para cualquier Z{displaystyle mathbf {Z} con probabilidad no cero p()Z▪ ▪ X,Silencio Silencio ){displaystyle p(mathbf {Z} mid mathbf {X}{boldsymbol {theta }}}}, podemos escribir

log⁡ ⁡ p()X▪ ▪ Silencio Silencio )=log⁡ ⁡ p()X,Z▪ ▪ Silencio Silencio )− − log⁡ ⁡ p()Z▪ ▪ X,Silencio Silencio ).{displaystyle log p(mathbf {X} mid {boldsymbol {theta })=log p(mathbf {X}mathbf {Z} mid {boldsymbol {theta })-log p(mathbf {Z} mid mathbf {X}{boldsymbol {theta }}}}}}}}}

Tomamos la expectativa sobre posibles valores de los datos desconocidos Z{displaystyle mathbf {Z} de la estimación del parámetro actual Silencio Silencio ()t){displaystyle theta ^{(t)} multiplicando ambos lados por p()Z▪ ▪ X,Silencio Silencio ()t)){displaystyle p(mathbf {Z} mid mathbf {X}{boldsymbol {theta }}}}} {(t)}}} y sumar (o integrar) Z{displaystyle mathbf {Z}. El lado izquierdo es la expectativa de una constante, así que obtenemos:

log⁡ ⁡ p()X▪ ▪ Silencio Silencio )=.. Zp()Z▪ ▪ X,Silencio Silencio ()t))log⁡ ⁡ p()X,Z▪ ▪ Silencio Silencio )− − .. Zp()Z▪ ▪ X,Silencio Silencio ()t))log⁡ ⁡ p()Z▪ ▪ X,Silencio Silencio )=Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))+H()Silencio Silencio ▪ ▪ Silencio Silencio ()t)),{fnMicrosoft Sans Serif}

Donde H()Silencio Silencio ▪ ▪ Silencio Silencio ()t)){displaystyle ¿Qué? se define por la suma negada que está reemplazando. Esta última ecuación sostiene por cada valor Silencio Silencio {displaystyle {boldsymbol {theta } incluido Silencio Silencio =Silencio Silencio ()t){displaystyle {boldsymbol {theta - Sí.,

log⁡ ⁡ p()X▪ ▪ Silencio Silencio ()t))=Q()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t))+H()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t)),{displaystyle log p(mathbf {X} mid {boldsymbol {theta} {} {} {fn]}=Q({boldsymbol {theta }} {t)}mid {boldsymbol {theta } {} {} {} {} {fnun} {f}} {fnun} {fnun}}}f}}fnfnun} {fnfnun}fnun}fnun}fnun}fnun}fnun}}fnun}fnun}fnun}fnun}}fnun}fnun}fnun}fnun}fnun}fnun}fnun}fnun}fnun}fnun}}}fnun}fnun}fnun}fnun}fnun}fnun}fnun}fnun

y restando esta última ecuación de la ecuación anterior da

log⁡ ⁡ p()X▪ ▪ Silencio Silencio )− − log⁡ ⁡ p()X▪ ▪ Silencio Silencio ()t))=Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))− − Q()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t))+H()Silencio Silencio ▪ ▪ Silencio Silencio ()t))− − H()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t)).{displaystyle log p(mathbf {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {f} {fnMicrosoft} {fnMicrosoft Sans Serif} {fnMicrosoft} {fnMicrosoft} {fnMicrosoft {fnMicrosoft}fnMicrosoft} {fnMicrosoft {f}}fnMicrosoft}fnMicrosoft}fnMicrosoft {fnMicrosoft {fnMicrosoft}fnMicrosoft}fnMicrosoft {fnMicrosoft {fnMicrosoft}fnMicrosoft}fnMicrosoft {f}fnMicrosoft}fnMicrosoft}fnMicrosoft {fnMicrosoft}fnun}fnMicrosoft {fnMicrosoft}fnMicrosoft {fnMicrosoft {fnMicro

Sin embargo, la desigualdad de Gibbs nos dice que H()Silencio Silencio ▪ ▪ Silencio Silencio ()t))≥ ≥ H()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t)){displaystyle ¿Qué?, así podemos concluir que

log⁡ ⁡ p()X▪ ▪ Silencio Silencio )− − log⁡ ⁡ p()X▪ ▪ Silencio Silencio ()t))≥ ≥ Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))− − Q()Silencio Silencio ()t)▪ ▪ Silencio Silencio ()t)).{displaystyle log p(mathbf {X} mid {boldsymbol {theta })-log p(mathbf {X} mid {boldsymbol {theta} {} {} {}} {} {} {} {fnuncio} {fnun} {fnuncio} {fnuncio}fnun} {fnun}fnun}fnfnun}fnfnun}fnun}fnfnfnfnfnfnun}fnun}fnfnun}fnfnfnun}fnun}fnfnun}fnun}fnun}fnun}nun}fnun}fnun}fnhnun}fnun}fnun}fnun}fnun}fnun}fnun}nhnun}nun}

En palabras, elegir Silencio Silencio {displaystyle {boldsymbol {theta } para mejorar Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t)){displaystyle Q({boldsymbol {theta }midboldsymbol {theta }}{(t)}}}} causas log⁡ ⁡ p()X▪ ▪ Silencio Silencio ){displaystyle log p(mathbf {X} mid {boldsymbol {theta}}} para mejorar al menos tanto.

Como procedimiento de maximización-maximización

El algoritmo EM puede verse como dos pasos de maximización alternos, es decir, como un ejemplo de descenso de coordenadas. Considere la función:

F()q,Silencio Silencio ):=Eq⁡ ⁡ [log⁡ ⁡ L()Silencio Silencio ;x,Z)]+H()q),{displaystyle F(q,theta):=operatorname {E} _{q}[log L(theta;x,Z)]+H(q),}

donde q es una distribución de probabilidad arbitraria sobre los datos no observados z y H(q) es la entropía de la distribución q. Esta función se puede escribir como

F()q,Silencio Silencio )=− − DKL()q∥ ∥ pZ▪ ▪ X()⋅ ⋅ ▪ ▪ x;Silencio Silencio ))+log⁡ ⁡ L()Silencio Silencio ;x),{displaystyle F(q,theta)=-D_{mathrm {KL} } {big (}qparallel p_{Zmid X}(cdot mid x;theta){big)}+log L(theta;x),}

Donde pZ▪ ▪ X()⋅ ⋅ ▪ ▪ x;Silencio Silencio ){displaystyle p_{Zmid X}(cdot mid x;theta)} es la distribución condicional de los datos no guardados dados los datos observados x{displaystyle x} y DKL{displaystyle D_{KL} es la divergencia Kullback-Leibler.

Entonces los pasos del algoritmo EM pueden verse como:

Paso de espera: Elija q{displaystyle q} para maximizar F{displaystyle F}:

q()t)=argmaxq⁡ ⁡ F()q,Silencio Silencio ()t)){displaystyle q^{(t)}=operatorname {arg,max} _{q} F(q,theta ^{(t)}}}

Paso de maximización: Elija Silencio Silencio {displaystyle theta } para maximizar F{displaystyle F}:

Silencio Silencio ()t+1)=argmaxSilencio Silencio ⁡ ⁡ F()q()t),Silencio Silencio ){displaystyle theta ^{(t+1)}=operatorname {arg,max} _{theta } F(q^{(t)},theta)}

Aplicaciones

EM se utiliza con frecuencia para la estimación de parámetros de modelos mixtos, especialmente en genética cuantitativa.

En psicometría, EM es una herramienta importante para estimar los parámetros de los ítems y las habilidades latentes de los modelos de teoría de respuesta al ítem.

Con la capacidad de lidiar con datos faltantes y observar variables no identificadas, los mercados emergentes se están convirtiendo en una herramienta útil para valorar y gestionar el riesgo de una cartera.

El algoritmo EM (y su variante más rápida de maximización de expectativas de subconjunto ordenado) también se usa ampliamente en la reconstrucción de imágenes médicas, especialmente en tomografía por emisión de positrones, tomografía computarizada por emisión de fotón único y tomografía computarizada por rayos X. Consulte a continuación otras variantes más rápidas de EM.

En ingeniería estructural, el algoritmo de identificación estructural mediante maximización de expectativas (STRIDE) es un método de salida únicamente para identificar las propiedades de vibración natural de un sistema estructural utilizando datos de sensores (consulte Análisis modal operativo).

EM también se utiliza para agrupar datos. En el procesamiento del lenguaje natural, dos ejemplos destacados del algoritmo son el algoritmo de Baum-Welch para modelos ocultos de Markov y el algoritmo de adentro hacia afuera para la inducción no supervisada de gramáticas probabilísticas libres de contexto.

En el análisis de los tiempos de espera entre operaciones, es decir, el tiempo entre operaciones posteriores de acciones en una bolsa de valores, el algoritmo EM ha demostrado ser muy útil.

Algoritmos EM de filtrado y suavizado

Por lo general, se utiliza un filtro de Kalman para la estimación del estado en línea y se puede emplear un suavizador de varianza mínima para la estimación del estado fuera de línea o por lotes. Sin embargo, estas soluciones de varianza mínima requieren estimaciones de los parámetros del modelo de espacio de estados. Los algoritmos EM se pueden utilizar para resolver problemas de estimación de parámetros y estados conjuntos.

Los algoritmos EM de filtrado y suavizado surgen al repetir este procedimiento de dos pasos:

E-step

Operar un filtro Kalman o un liso mínimo de variación diseñado con estimaciones actuales del parámetro para obtener estimaciones actualizadas del estado.

M-step

Utilice las estimaciones de estado filtradas o suavizadas dentro de los cálculos de probabilidad máxima para obtener estimaciones actualizadas del parámetro.

Supongamos que un filtro de Kalman o un suavizador de varianza mínima funciona en mediciones de un sistema de una sola entrada y una sola salida que posee ruido blanco aditivo. Se puede obtener una estimación actualizada de la varianza del ruido de medición a partir del cálculo de máxima verosimilitud.

σ σ ^ ^ v2=1N.. k=1N()zk− − x^ ^ k)2,{displaystyle {widehat {sigma - ¿Qué? {1}{N}sum} {fnMicrosoft Sans Serif} {x}_{k}}} {2}}

Donde x^ ^ k{displaystyle {widehat {x}_{k}} son estimaciones de salida de escalar calculadas por un filtro o un más suave de las mediciones de escalar N zk{displaystyle z_{k}. La actualización anterior también se puede aplicar para actualizar una intensidad de ruido de medición Poisson. Del mismo modo, para un proceso auto-regresivo de primer orden, se puede calcular una estimación actualizada de la varianza de ruido del proceso mediante

σ σ ^ ^ w2=1N.. k=1N()x^ ^ k+1− − F^ ^ x^ ^ k)2,{displaystyle {widehat {sigma ¿Qué? {1}{N}sum} {fnMicrosoft Sans Serif} {x}_{k+1}-{widehat {F} {fnK} {fnMicrosoft}} {fnK}} {f}} {f}} {f}} {f}}}} {f}}}} {f}}}} {f}}}} {f}}}} {f}}}}}}}}} { {x}_{k}}} {2}}

Donde x^ ^ k{displaystyle {widehat {x}_{k}} y x^ ^ k+1{displaystyle {widehat {x}_{k+1} son estimaciones del estado del escalar calculadas por un filtro o un liso. El cálculo actualizado del coeficiente se obtiene mediante

F^ ^ =.. k=1N()x^ ^ k+1− − F^ ^ x^ ^ k)2.. k=1Nx^ ^ k2.{displaystyle {widehat {F}={frac} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fn} {fnK} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnK} {fnK} {fnKfnKfnKf}}fnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnKfnK}}}}fnK {x}_{k+1}-{widehat {F} {fnK} {fnMicrosoft}} {fnK}} {f}} {f}} {f}} {f}}}} {f}}}} {f}}}} {f}}}} {f}}}} {f}}}}}}}}} { {fnK}} {fnK}} {fnK}} {fnK}}} {f}}} {fn}}} {fn}}} {fnK}}}}}}}}} {fn}}}}}} {f}}}}}}}} {f}}}}}}} {f}}} {f}}}}}}}}}}}}}}}}} {f}}}}}}}}}}}}}}}}}}}} {}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} {}}}}}}}}}}}}} {}}}}}}} {f}} {}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} { ###{k=1} {fn} {fnMicrosoft Sans Serif} {x}_{k}}}}

La convergencia de estimaciones de parámetros como las anteriores está bien estudiada.

Variantes

Se han propuesto varios métodos para acelerar la convergencia a veces lenta del algoritmo EM, como los que utilizan gradiente conjugado y los métodos de Newton modificados (Newton-Raphson). Además, EM se puede utilizar con métodos de estimación restringidos.

maximización de expectativas ampliadas por parámetros (PX-EM) a menudo proporciona aceleración al "usar un `ajuste de covarianza&#39"; corregir el análisis del paso M, aprovechando la información adicional capturada en los datos completos imputados".

Maximización condicional de expectativas (ECM) reemplaza cada M paso con una secuencia de pasos de maximización condicional (CM) en los que cada parámetro θ_i se maximiza individualmente, condicionado a que los demás parámetros permanezcan fijos. Se puede ampliar al algoritmo Maximización condicional de expectativas (ECME).

Esta idea se amplía aún más en el algoritmo maximización de expectativas generalizada (GEM), en el que se busca sólo un aumento en la función objetivo F tanto para el paso E como para M. paso como se describe en la sección Como procedimiento de maximización-maximización. GEM se desarrolla aún más en un entorno distribuido y muestra resultados prometedores.

También es posible considerar el algoritmo EM como una subclase del algoritmo MM (Mayorizar/Minimizar o Minorizar/Maximizar, según el contexto) y, por tanto, utilizar cualquier maquinaria desarrollada en el algoritmo más general. caso.

Algoritmo Α-EM

La función Q utilizada en el algoritmo EM se basa en la probabilidad logarítmica. Por lo tanto, se considera el algoritmo log-EM. El uso del registro de probabilidad se puede generalizar al de la razón de probabilidad logarítmica α. Entonces, la razón de probabilidad logarítmica α de los datos observados se puede expresar exactamente como igualdad utilizando la función Q de la razón de probabilidad logarítmica α y la divergencia α. La obtención de esta función Q es un paso E generalizado. Su maximización es un paso M generalizado. Este par se llama algoritmo α-EM. que contiene el algoritmo log-EM como subclase. Por tanto, el algoritmo α-EM de Yasuo Matsuyama es una generalización exacta del algoritmo log-EM. No es necesario calcular el gradiente ni la matriz de Hesse. El α-EM muestra una convergencia más rápida que el algoritmo log-EM al elegir un α apropiado. El algoritmo α-EM conduce a una versión más rápida del algoritmo de estimación del modelo oculto de Markov α-HMM.

Relación con los métodos variacionales de Bayes

EM es un método de máxima verosimilitud parcialmente no bayesiano. Su resultado final proporciona una distribución de probabilidad sobre las variables latentes (en el estilo bayesiano) junto con una estimación puntual para θ (ya sea una estimación de máxima verosimilitud o una moda posterior). Es posible que desee una versión completamente bayesiana de esto, que proporcione una distribución de probabilidad sobre θ y las variables latentes. El enfoque bayesiano de la inferencia consiste simplemente en tratar θ como otra variable latente. En este paradigma, la distinción entre los pasos E y M desaparece. Si se utiliza la aproximación Q factorizada como se describe anteriormente (Bayes variacional), la resolución puede iterar sobre cada variable latente (que ahora incluye θ) y optimizarlas una a la vez. Ahora, se necesitan k pasos por iteración, donde k es el número de variables latentes. Para los modelos gráficos, esto es fácil de hacer ya que la nueva Q de cada variable depende solo de su manta de Markov, por lo que el paso de mensajes local se puede utilizar para una inferencia eficiente.

Interpretación geométrica

En geometría de la información, el paso E y el paso M se interpretan como proyecciones bajo conexiones afines duales, llamadas conexión e y conexión m; La divergencia Kullback-Leibler también puede entenderse en estos términos.

Ejemplos

Mezcla gaussiana

Vamos x=()x1,x2,...... ,xn){displaystyle mathbf {x} =(mathbf {x} _{1},mathbf {x} _{2},ldotsmathbf {x} ser una muestra de n{displaystyle n} observaciones independientes de una mezcla de dos distribuciones normales multivariadas de dimensión d{displaystyle d}, y dejar z=()z1,z2,...... ,zn){displaystyle mathbf {z} =(z_{1},z_{2},ldotsz_{n}} ser las variables latentes que determinan el componente desde el cual se origina la observación.

Xi▪ ▪ ()Zi=1)♪ ♪ Nd()μ μ 1,.. 1){displaystyle X_{i}mid (Z_{i}=1)sim {mathcal {N}_{d}({boldsymbol {mu} }_{1},Sigma _{1})} y Xi▪ ▪ ()Zi=2)♪ ♪ Nd()μ μ 2,.. 2),{displaystyle X_{i}mid (Z_{i}=2)sim {mathcal {N}_{d}({boldsymbol {mu} - Sí.

dónde

P⁡ ⁡ ()Zi=1)=τ τ 1{displaystyle operatorname {P} (Z_{i}=1)=tau _{1},} y P⁡ ⁡ ()Zi=2)=τ τ 2=1− − τ τ 1.{displaystyle operatorname (Z_{i}=2)= # {2}=1-tau _{1}

El objetivo es estimar los parámetros desconocidos que representan el valor de mezcla entre las gaussianas y las medias y covarianzas de cada una:

Silencio Silencio =()τ τ ,μ μ 1,μ μ 2,.. 1,.. 2),{displaystyle theta ={big (}{boldsymbol {tau },{boldsymbol {mu }_{1},{boldsymbol {mu - Sí.

donde está la función de probabilidad de datos incompletos

L()Silencio Silencio ;x)=∏ ∏ i=1n.. j=12τ τ jf()xi;μ μ j,.. j),{displaystyle L(theta;mathbf {x}=prod - ¿Qué? ¿Por qué? ##{i};{boldsymbol {mu} - Sí.

y la función de probabilidad de datos completos es

L()Silencio Silencio ;x,z)=p()x,z▪ ▪ Silencio Silencio )=∏ ∏ i=1n∏ ∏ j=12[f()xi;μ μ j,.. j)τ τ j]I()zi=j),{displaystyle L(theta;mathbf {x}mathbf {z}=p(mathbf {x}mathbf {z} mid theta)=prod ¿Qué? ¿Qué? _{i};{boldsymbol # ¿Qué? {I} (z_{i}=j)}}

o

L()Silencio Silencio ;x,z)=exp⁡ ⁡ {}.. i=1n.. j=12I()zi=j)[log⁡ ⁡ τ τ j− − 12log⁡ ⁡ Silencio.. jSilencio− − 12()xi− − μ μ j)⊤ ⊤ .. j− − 1()xi− − μ μ j)− − d2log⁡ ⁡ ()2π π )]},{displaystyle L(theta;mathbf {x}mathbf {z}=exp left{sum - ¿Qué? - ¿Qué? {I} (z_{i}=j){big [}log tau _{j}-{tfrac {1}{2}log Silencioso, cariño. ¿Por qué? _{i}-{boldsymbol # }_{j})-{tfrac {d}{2}log(2pi){big}right}}

Donde I{displaystyle mathbb {I} es una función indicadora y f{displaystyle f} es la función de densidad de probabilidad de una normalidad multivariada.

En la última igualdad, para cada uno i, un indicador I()zi=j){displaystyle mathbb {I} (z_{i}=j)} es igual a cero, y un indicador es igual a uno. La suma interna reduce así a un término.

Paso E

Dada nuestra estimación actual de los parámetros θ^(t), la distribución condicional de la Z_i está determinado por el teorema de Bayes como la altura proporcional de la densidad normal ponderada por τ:

Tj,i()t):=P⁡ ⁡ ()Zi=j▪ ▪ Xi=xi;Silencio Silencio ()t))=τ τ j()t)f()xi;μ μ j()t),.. j()t))τ τ 1()t)f()xi;μ μ 1()t),.. 1()t))+τ τ 2()t)f()xi;μ μ 2()t),.. 2()t)).{displaystyle T_{j,i}{(t)}:=operatorname [P] (Z_{i}=jmid X_{i}=mathbf {x} _{i};theta ^{(t)}={frac {tau _{j}{(t)} f(mathbf {x} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {f}} {fnMicrosoft Sans Serif} _{i},Sigma _{(t)})+tau _{2}{(t)} f(mathbf {x} - Sí.

These are called the "membership probabilities#34;, which are normally considered the output of the E step (although this is not the Q function of below).

Este paso E corresponde con la configuración de esta función para Q:

Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))=EZ▪ ▪ X=x;Silencio Silencio ()t)⁡ ⁡ [log⁡ ⁡ L()Silencio Silencio ;x,Z)]=EZ▪ ▪ X=x;Silencio Silencio ()t)⁡ ⁡ [log⁡ ⁡ ∏ ∏ i=1nL()Silencio Silencio ;xi,Zi)]=EZ▪ ▪ X=x;Silencio Silencio ()t)⁡ ⁡ [.. i=1nlog⁡ ⁡ L()Silencio Silencio ;xi,Zi)]=.. i=1nEZi▪ ▪ Xi=xi;Silencio Silencio ()t)⁡ ⁡ [log⁡ ⁡ L()Silencio Silencio ;xi,Zi)]=.. i=1n.. j=12P()Zi=j▪ ▪ Xi=xi;Silencio Silencio ()t))log⁡ ⁡ L()Silencio Silencio j;xi,j)=.. i=1n.. j=12Tj,i()t)[log⁡ ⁡ τ τ j− − 12log⁡ ⁡ Silencio.. jSilencio− − 12()xi− − μ μ j)⊤ ⊤ .. j− − 1()xi− − μ μ j)− − d2log⁡ ⁡ ()2π π )].{displaystyle {begin{aligned}Q(theta mid theta ^{(t)}) {E} _{mathbf {Z} mid mathbf {X} =mathbf {x};mathbf {theta } { {(t)} {log L(theta;mathbf {x}mathbf {Z})}\\cH00=operatorname {E} _{mathbf {Z} mid mathbf {X} =mathbf {x};mathbf {theta } {(t)}[log prod] ¿Por qué? ########################### {E} _{mathbf {Z} mid mathbf {X} =mathbf {x};mathbf {theta} {} {} {sum _{i=1}log L(theta;mathbf {x} - ¿Por qué? ################################################################################################################################################################################################################################################################ {E}{Z_{i}mid} [log L(theta] - ¿Por qué? - ¿Por qué? - ¿Por qué? X_{i}=mathbf {x} _{i};theta ^{(t)})log L(theta ¿Qué? - ¿Qué? - ¿Por qué? ¿Por qué? {1}{2}log Silencioso, cariño. ¿Por qué? _{i}-{boldsymbol # }_{j})-{tfrac {d} {2}log(2pi){big]}end{aligned}}

La expectativa de log⁡ ⁡ L()Silencio Silencio ;xi,Zi){displaystyle log L(theta;mathbf {x} ¿Qué? dentro de la suma se toma con respecto a la función de densidad de probabilidad P()Zi▪ ▪ Xi=xi;Silencio Silencio ()t)){displaystyle P(Z_{i}mid ¿Qué?, que puede ser diferente para cada xi{displaystyle mathbf {x} _{i} del conjunto de entrenamiento. Todo en el paso E se conoce antes de que se tome el paso excepto Tj,i{displaystyle T_{j,i}, que se calcula según la ecuación al comienzo de la sección E paso.

No es necesario calcular esta expectativa condicional completa en un solo paso, porque τ y μ/Σ aparecen en términos lineales separados y De este modo se puede maximizar de forma independiente.

M paso

Q(θ | θ^(t)) siendo cuadrático en forma significa que determinar los valores maximizadores de θ es relativamente sencillo. Además, τ, (μ₁,Σ₁) y (μ₂,Σ₂) pueden maximizarse independientemente ya que todos aparecen en términos lineales separados.

To begin, consider τ, which has the constraint τ₁ + τ₂=1:

τ τ ()t+1)=argmaxτ τ Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))=argmaxτ τ {}[.. i=1nT1,i()t)]log⁡ ⁡ τ τ 1+[.. i=1nT2,i()t)]log⁡ ⁡ τ τ 2}.{fnMicrosoft } {fnMicrosoft } {fnMicrosoft }} {fnMicrosoft}}} {fnunci}} {fnMicrosoft } {fnMicrosoft}}}\fnMicrosoft} Q(theta mid theta ^{(t)})\ tropezón {boldsymbol {tau {fnMicrosoft Sans Serif}left{left[sum] _{i=1}{n} {=t)}log tau _{1}+left[sum _{i=1}}{n} {n}T_{2,i}{(t)}rightlog tau _{2}right}end{aligned}}}}}}}}

This has the same form as the MLE for the binomial distribution, sop>

τ τ j()t+1)=.. i=1nTj,i()t).. i=1n()T1,i()t)+T2,i()t))=1n.. i=1nTj,i()t).{displaystyle tau _{(t+1)}={frac {sum ¿Qué? ¿Por qué? ¿Qué?

Para las siguientes estimaciones de (μ₁,Σ₁):

()μ μ 1()t+1),.. 1()t+1))=argmaxμ μ 1,.. 1Q()Silencio Silencio ▪ ▪ Silencio Silencio ()t))=argmaxμ μ 1,.. 1.. i=1nT1,i()t){}− − 12log⁡ ⁡ Silencio.. 1Silencio− − 12()xi− − μ μ 1)⊤ ⊤ .. 1− − 1()xi− − μ μ 1)}.{displaystyle {begin{aligned}({boldsymbol {mu} ################################################################################################################################################################################################################################################################ }_{1},Sigma ¿Por qué? Q(theta mid theta ^{(t)})\ tropezón {{boldsymbol {mu] }_{1},Sigma ################################################################################################################################################################################################################################################################ {1}{2}log Silencio _{1} tortura-{2} {x} _{i}-{boldsymbol {fnMicrosoft Sans Serif} _{i}-{boldsymbol # - Bien.

Esto tiene la misma forma que un MLE ponderado para una distribución normal, por lo que

μ μ 1()t+1)=.. i=1nT1,i()t)xi.. i=1nT1,i()t){displaystyle {boldsymbol {mu} {fnMicrosoft _{i=1} {fn} {fn} {fn} {fn} {cH00}}}mtbf {x} ¿Qué? ¿Qué? y .. 1()t+1)=.. i=1nT1,i()t)()xi− − μ μ 1()t+1))()xi− − μ μ 1()t+1))⊤ ⊤ .. i=1nT1,i()t){displaystyle Sigma _{1}{(t+1)}={frac {sum _{i=1}^{n}T_{1,i}{(t)}(mathbf {x} ¿Por qué? _{i}-{boldsymbol {fnK} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} ¿Qué?

y, por simetría,

μ μ 2()t+1)=.. i=1nT2,i()t)xi.. i=1nT2,i()t){displaystyle {boldsymbol {mu} {fnMicroc {fnMicrosoft {fnMicroc} ¿Qué? ¿Qué? ¿Qué? y .. 2()t+1)=.. i=1nT2,i()t)()xi− − μ μ 2()t+1))()xi− − μ μ 2()t+1))⊤ ⊤ .. i=1nT2,i()t).{displaystyle Sigma _{2}{(t+1)}={frac {sum _{i=1}^{n}T_{2,i}{(t)} {mathbf {x} ¿Por qué? _{i}-{boldsymbol # {fnMicrosoft Sans Serif} - ¿Qué?

Terminación

Conclusión del proceso iterativo si EZ▪ ▪ Silencio Silencio ()t),x[log⁡ ⁡ L()Silencio Silencio ()t);x,Z)]≤ ≤ EZ▪ ▪ Silencio Silencio ()t− − 1),x[log⁡ ⁡ L()Silencio Silencio ()t− − 1);x,Z)]+ε ε {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {f} {f} {fnMicrosoft}f}cH0}f}cH00} {cH0}cH00}cH00} para ε ε {displaystyle varepsilon } debajo de un umbral preestablecido.

Generalización

El algoritmo ilustrado anteriormente se puede generalizar para mezclas de más de dos distribuciones normales multivariadas.

Regresión truncada y censurada

El algoritmo EM se ha implementado en el caso en el que existe un modelo de regresión lineal subyacente que explica la variación de alguna cantidad, pero donde los valores realmente observados son versiones censuradas o truncadas de los representados en el modelo. Los casos especiales de este modelo incluyen observaciones censuradas o truncadas de una distribución normal.

Alternativas

Los EM típicamente convergen a un óptimo local, no necesariamente al óptimo global, sin límite en la tasa de convergencia en general. Es posible que pueda ser arbitrariamente pobre en dimensiones altas y que pueda haber un número exponencial de óptimos locales. Por lo tanto, existe la necesidad de métodos alternativos para garantizar el aprendizaje, especialmente en entornos de alta dimensión. Existen alternativas a los EM con mejores garantías de coherencia, que se denominan enfoques basados en momentos o las llamadas técnicas espectrales. Los enfoques basados en momentos para aprender los parámetros de un modelo probabilístico son de creciente interés recientemente, ya que disfrutan de garantías como la convergencia global bajo ciertas condiciones, a diferencia de los mercados emergentes, que a menudo se ven afectados por el problema de quedarse estancados en óptimos locales. Se pueden derivar algoritmos con garantías de aprendizaje para una serie de modelos importantes, como modelos de mezcla, HMM, etc. Para estos métodos espectrales, no se producen óptimos locales espurios y los parámetros verdaderos se pueden estimar consistentemente bajo algunas condiciones de regularidad.