Traducción al Castellano

Capítulo 3 - Sección 3

SIN EDICION (Ng)

Introducción

En estas notas, discutimos el algoritmo EM (Expectation-Maximization) para la estimación de densidad. Este algoritmo es ampliamente utilizado en problemas de estimación donde las variables latentes dificultan la optimización directa. El EM permite manejar este tipo de problemas dividiendo el proceso en dos etapas iterativas: el paso de “esperanza” (E-step) y el paso de “maximización” (M-step).

FALTA MOTIVACION. El da uno lindo…Mide cosas, y quiere detectar anomalias, entonces vos queres estimar la densidad para saber si tu punto esta en un lugar de la densidad es muy baja. Pero no es normal, pero se puede pensar como una mezcla de gaussiana. Un grafico seria bueno. Es muy lindocomo cuneta que cadda valor puede no ser raro pero la combinacion si (eso me parece que ilumina).

EM para Mezcla de Gaussianas

ACA MOTIVA de nuevo lindo, con datos unidimensionales viniendo de dos gaussianas. Supone que sabe los puntos de que gausiana entonces luego estima…esto hasta seria lindo hacer las cuentas para dos conjuntos de puntos unidimemsnionales que sabemos de donde vienen.

Supongamos que tenemos un conjunto de entrenamiento $\{x^{(1)}, \ldots, x^{(n)}\}$. Dado que estamos en el contexto de aprendizaje no supervisado, estos puntos no tienen etiquetas.

Para ello suponemos que hay variables $z^{(i)}$ no observables y modelamos los datos especificando una distribución conjunta $p(x^{(i)}, z^{(i)}) = p(x^{(i)}|z^{(i)})p(z^{(i)})$ . Aquí:

$z^{(i)} \sim \text{Multinomial}(\phi)$ , donde:
- $\phi_j \geq 0$ ,
- $\sum_{j=1}^k \phi_j = 1$ , y
- el parámetro $\phi_j=p(z^{(i)} = j)$ .
Condicional a $z^{(i)} = j$ , se tiene $x^{(i)}|z^{(i)} = j \sim \mathcal{N}(\mu_j, \Sigma_j)$.

Dejamos que $k$ sea el número de valores que $z^{(i)}$ puede tomar. Así, nuestro modelo establece que cada $x^{(i)}$ fue generado eligiendo aleatoriamente $z^{(i)}$ de $\{1, \ldots, k$ , y luego $x^{(i)}$ fue extraido de una de $k$ gaussianas dependiendo de $z^{(i)}$ . Esto se denomina mezcla de gaussianas. También, note que las $z^{(i)}$ son variables aleatorias latentes, es decir, no son observadas directamente. Este hecho introduce complejidad en el proceso de estimación de parámetros.

Parece Análisis Gausiano discriminante, pero tiene dos diferencias: - $z^{(i)} \sim \text{Multinomial}(\phi)$, mixto de gausiana con dos clases es Bernulli, pero AGD tambien hay con mas clases, no seria lo mismo? deberiamos dar AGD con mas de dos clases si no esta dado. HABLAR. Para mi no hay diferenica. - Las gausianas tienen diferente varianza, pero de nuevo en analisis discriminante gausiano cuadratico es asi.

Para mi la diferencia es que aca no sabes de donde vienen los casos, punto.

Los parámetros de nuestro modelo son $\phi$, $\mu$ y $\Sigma$ . Para estimarlos, escribimos la verosimilitud de nuestros datos:

\[ \begin{align*} \ell(\phi, \mu, \Sigma) &= \sum_{i=1}^n \log p(x^{(i)}; \phi, \mu, \Sigma)\\ & = \sum_{i=1}^n \log \sum_{z^{(i)}=1}^k p(x^{(i)}|z^{(i)}; \mu, \Sigma)p(z^{(i)}; \phi). \end{align*} \]

No es posible encontrar estimaciones de máxima verosimilitud en forma cerrada resolviendo las derivadas de esta ecuación. Sin embargo, si supiéramos las $z^{(i)}$, la verosimilitud sería fácil de calcular:

\[ \ell(\phi, \mu, \Sigma) = \sum_{i=1}^n \log p(x^{(i)}|z^{(i)}; \mu, \Sigma) + \log p(z^{(i)}; \phi). \]

Maximizando con respecto a $\phi$, $\mu$ y $\Sigma$ da los parámetros:

\[ \begin{align*} \phi_j &= \frac{1}{n} \sum_{i=1}^n 1\{z^{(i)} = j\}, \\ \mu_j &= \frac{\sum_{i=1}^n 1\{z^{(i)} = j\}x^{(i)}}{\sum_{i=1}^n 1\{z^{(i)} = j\}}, \\ \Sigma_j &= \frac{\sum_{i=1}^n 1\{z^{(i)} = j\}(x^{(i)} - \mu_j)(x^{(i)} - \mu_j)^T}{\sum_{i=1}^n 1\{z^{(i)} = j\}}. \end{align*} \]

Estas expresiones permiten entender que si las variables latentes $z^{(i)}$ fueran conocidas, la estimación de los parámetros sería directa, similar a otros problemas supervisados como el Análisis Discriminante Gaussiano.

Sin embargo, dado que las $z^{(i)}$ no son conocidas, necesitamos un método iterativo para estimarlas y, a su vez, refinar los parámetros del modelo. Aquí es donde entra en juego el algoritmo EM.

El Algoritmo EM

El algoritmo EM es una solución iterativa para problemas con variables latentes. Tiene dos pasos principales:

Paso E (Expectation): “Adivinar” los valores de las $z^{(i)}$ calculando su valor esperado dado los parámetros actuales.
Paso M (Maximization): Actualizar los parámetros $\phi$ , $\mu$ y $\Sigma$ maximizando la verosimilitud con base en las “adivinanzas” del paso anterior.

Procedimiento Detallado

Repetir hasta la convergencia:

Paso E: Para cada $i, j$ , calcular:

\[ \begin{align*} w_j^{(i)} &:= p(z^{(i)} = j|x^{(i)}; \phi, \mu, \Sigma) \\&= \frac{p(x^{(i)}|z^{(i)} = j; \mu, \Sigma)p(z^{(i)} = j; \phi)}{\sum_{l=1}^k p(x^{(i)}|z^{(i)} = l; \mu, \Sigma)p(z^{(i)} = l; \phi)}. \end{align*} \]

Estos valores representan la probabilidad a posteriori de que el dato $x^{(i)}$ provenga de la componente $j$ de la mezcla.

Paso M: Actualizar los parámetros usando los pesos $w_j^{(i)}$:

\[ \begin{align*} \phi_j &= \frac{1}{n} \sum_{i=1}^n w_j^{(i)}, \\ \mu_j &= \frac{\sum_{i=1}^n w_j^{(i)}x^{(i)}}{\sum_{i=1}^n w_j^{(i)}}, \\ \Sigma_j &= \frac{\sum_{i=1}^n w_j^{(i)}(x^{(i)} - \mu_j)(x^{(i)} - \mu_j)^T}{\sum_{i=1}^n w_j^{(i)}}. \end{align*} \]

El algoritmo se detiene cuando las actualizaciones de los parámetros son suficientemente pequeñas, indicando convergencia.

Interpretación del Algoritmo

El paso E permite calcular las probabilidades a posteriori de las variables latentes $z^{(i)}$ usando la regla de Bayes. Estas probabilidades “suaves” reemplazan las asignaciones “duras” del algoritmo $k$-medias, lo que hace al EM más flexible.

El paso M es una optimización directa de los parámetros $\phi$, $\mu$, y $\Sigma$ asumiendo que las probabilidades $w_j^{(i)}$ son correctas. Además, deberías comparar las actualizaciones en el paso M con las fórmulas que teníamos cuando los valores de eran conocidos exactamente. Son idénticas, excepto que en lugar de las funciones indicadoras , que indicaban de qué distribución gaussiana provenía cada punto de datos, ahora tenemos en su lugar los valores .

Entonces, una intuición sobre este algoritmo de mezcla de Gaussianas es que es un poco como el k-means, pero con asignaciones suaves. En k-means, en el primer paso, tomaríamos cada punto y lo asignaríamos a uno de los k centroides del clúster, Y si estuviera un poco más cerca del centroide del clúster rojo que del centroide del clúster azul, simplemente lo asignaríamos al centroide del clúster rojo. Entonces, incluso si estuviera solo un poco más cerca de un centroide que de otro, k-means simplemente hace lo que se llama una asignación dura. Es decir, asigna el punto al 100% al centroide del clúster más cercano.

Por otro lado, EM implementa una forma más suave de asignar puntos a los diferentes centroides de los clústeres. En lugar de elegir simplemente el centro de la Gaussiana más cercana y asignarlo ahí, utiliza estas probabilidades y les da un peso, en términos de cuánto se asigna a la Gaussiana 1 frente a la Gaussiana 2.

En el segundo paso, actualiza las medias en consecuencia. Suma todos los $x^{(i)}$ según la medida en que se asignan a ese centroide del clúster, dividido por el número de ejemplos asignados a ese centroide. Entonces, esa es una intuición sobre la relación entre EM y k-means.

Cuando ejecutas este algoritmo, resulta que este algoritmo convergerá, con algunas salvedades que mencionaré más adelante, y encontrará una estimación bastante decente de los parámetros, por ejemplo, al ajustar un modelo de mezcla de dos Gaussianas.

Propiedades del Algoritmo

Garantiza convergencia a un máximo local de la verosimilitud.
Puede ser sensible a las condiciones iniciales (similar a $k$-medias), por lo que se recomienda inicializarlo varias veces con diferentes valores.

El algoritmo EM proporciona una metodología general para problemas de estimación donde hay variables latentes y es ampliamente aplicable en diversas áreas como mezcla de gaussianas, inferencia en modelos bayesianos y estimación de parámetros ocultos en cadenas de Markov.

VUELVE ACA AL EJEMPLO MOTIVADOR DE DETECCION DE ANOMALIDADES Y CUENTA COSAS (Esta bueno)

Funciona? (esto habria que dar algo no se si esto esta bien)

Medir si el algoritmo EM para mezcla de Gaussianas “funciona” depende del contexto, los datos y el objetivo del análisis. Aquí hay varias métricas y métodos para evaluar su desempeño:

Visualización de Clusters

Método: Observa si los clusters formados reflejan estructuras coherentes en los datos, especialmente si las distribuciones gaussianas se ajustan bien.
Herramienta: Graficar los datos con colores según los clusters asignados y las gaussianas ajustadas (como en los gráficos de elipses).
Interpretación: Si los puntos en un cluster están cerca entre sí y las gaussianas parecen ajustarse a la densidad real de los datos, el modelo probablemente funciona.

Log-Verosimilitud

Método: Observa la evolución de la log-verosimilitud durante las iteraciones del algoritmo EM.
Interpretación: Aumento constante: Indica que el algoritmo está ajustándose progresivamente a los datos. Convergencia estable: Indica que el modelo ha alcanzado un máximo local de verosimilitud. Si la log-verosimilitud disminuye o no converge, puede ser un signo de problemas en los datos o en la inicialización.

Validación Cruzada (Evaluación Cuantitativa)

Método: Divide los datos en un conjunto de entrenamiento y otro de prueba. Ajusta el modelo con el conjunto de entrenamiento y evalúa la log-verosimilitud en el conjunto de prueba.
Interpretación: Una log-verosimilitud alta en los datos de prueba sugiere que el modelo generaliza bien y capta la estructura subyacente.

Coherencia Interna

Método: Mide si los puntos dentro de un cluster son más similares entre sí que con puntos de otros clusters. Métrica:
Índice de Silueta: Mide qué tan cerca están los puntos de su cluster frente a otros clusters. Valores cercanos a 1 indican buenos clusters.
Inercia o SSE (Suma de Errores Cuadrados): Mide la compactación de los clusters (similar a k-means).

Evaluación Supervisada (Si Hay Etiquetas)

Si los datos tienen etiquetas reales (como en el conjunto de Iris):

Método: Compara las asignaciones de clusters del modelo EM con las etiquetas reales.
Métricas: Precisión y Recall: Proporción de etiquetas correctas y falsos positivos. ARI (Adjusted Rand Index): Mide qué tan similares son las asignaciones del modelo con las etiquetas reales. NMI (Normalized Mutual Information): Mide qué tan dependientes son los clusters generados de las etiquetas.

Robustez ante Diferentes Inicializaciones

Método: Ejecuta el algoritmo varias veces con diferentes inicializaciones y compara los resultados (clusters formados, log-verosimilitud final).
Interpretación: Si los resultados son consistentes, el modelo es robusto. Si los resultados varían significativamente, puede indicar que el modelo está atrapado en óptimos locales.

Evaluación Manual o Contextual

Método: Evalúa si los clusters obtenidos tienen sentido en el contexto del problema (por ejemplo, categorías de usuarios o segmentos de mercado).
Herramienta: Consulta con expertos en el dominio o analiza patrones específicos en los clusters.

Implementación del Algoritmo EM y aplicación a un ejemplo simulado y un ejemplo Real

Convergencia a un mínimo local del algoritmo EM

En el caso particular de las mezclas de Gaussianas con el algoritmo EM, hay una forma más formal de derivar el algoritmo EM que demuestra que este es un algoritmo de estimación de máxima verosimilitud y que converge al menos en un óptimo local.

En particular, lo que haremos será mostrar que, si tu objetivo es dado un modelo $P(x,z)$ parametrizado por $\theta$, maximizar $ $P(x)$.

Esto es lo que la máxima verosimilitud se supone que debe hacer. Y el algoritmo EM está exactamente intentando hacer eso, ¿de acuerdo?

Entonces, presentaremos esta derivación más general, la forma general de derivar el algoritmo EM que no depende de este argumento intuitivo de que es más fácil usar máxima verosimilitud con valores iniciales estimados.

Desigualdad de Jensen (capaz se enseño antes, igual pongamoslo aca por las duads y luego lo pasamos)

Comenzamos nuestra discusión con un resultado muy útil llamado Desigualdad de Jensen.

Sea $f$ una función cuyo dominio sea el conjunto de los números reales. Recordemos que $f$ es una función convexa si su segunda derivada $f''(x) \geq 0$ (para todo $x \in \mathbb{R}$). En el caso de que $f$ tome entradas vectoriales, esto se generaliza a la condición de que su hessiano $H$ sea semidefinido positivo ($H \geq 0$). Si $f''(x) > 0$ para todo $x$, decimos que $f$ es estrictamente convexa (en el caso vectorial, esto corresponde a que $H$ debe ser definido positivo, denotado como $H > 0$). La desigualdad de Jensen puede ser enunciada como sigue:

Teorema. Sea $f$ una función convexa y sea $X$ una variable aleatoria. Entonces:

\[ \mathbb{E}[f(X)] \geq f(\mathbb{E}[X]). \]

Además, si $f$ es estrictamente convexa, entonces $\mathbb{E}[f(X)] = f(\mathbb{E}[X])$ sólo si y sólo si $X = \mathbb{E}[X]$ con probabilidad 1 (es decir, si $X$ es constante).

Recordemos nuestra convención de omitir ocasionalmente los paréntesis al escribir esperanzas, por lo que en el teorema anterior, $f(\mathbb{E}[X]) = f(\mathbb{E}[X])$.

Para interpretar el teorema, consideremos la figura a continuación:

Aquí, $f$ es una función convexa representada por la línea sólida. Además, $X$ es una variable aleatoria que tiene una probabilidad de 0.5 de tomar el valor $a$ y una probabilidad de 0.5 de tomar el valor $b$ (indicados en el eje $x$). Así, el valor esperado de $X$ está dado por el punto medio entre $a$ y $b$.

También vemos los valores $f(a)$, $f(b)$ y $f(\mathbb{E}[X])$ indicados en el eje $y$. Además, el valor $\mathbb{E}[f(X)]$ ahora es el punto medio en el eje $y$ entre $f(a)$ y $f(b)$. A partir de este ejemplo, vemos que, debido a que $f$ es convexa, debe cumplirse que $\mathbb{E}[f(X)] \geq f(\mathbb{E}[X])$.

Incidentalmente, muchas personas tienen problemas para recordar en qué dirección va la desigualdad, y recordar un gráfico como este es una buena manera de resolver rápidamente la duda.

Observación. Recordemos que $f$ es [estrictamente] cóncava si y sólo si $-f$ es [estrictamente] convexa (es decir, $f''(x) \leq 0$ o $H \leq 0$). La desigualdad de Jensen también se cumple para funciones cóncavas, pero con la dirección de todas las desigualdades invertida:

\[ \mathbb{E}[f(X)] \leq f(\mathbb{E}[X]), \text{ etc.} \]

La línea recta es convexa y cóncava (no estricta). La igualdad en Jensen se da si y solo si $X$ es casi siempre igual a constante. (Ponerlo como ejercicio?). O $X=E(X)$ con probabilidad 1.

Algoritmos Generales EM

Introducción

Supongamos que tenemos un problema de estimación en el cual disponemos de un conjunto de entrenamiento $\{x^{(1)}, \ldots, x^{(n)}\}$ que consiste en $n$ ejemplos independientes. Este conjunto puede ser observado, pero también existe una variable latente $z$, no observable, que dificulta la estimación directa de los parámetros. Este tipo de problemas son comunes en muchas áreas, como el modelado de mezclas gaussianas o modelos ocultos de Markov.

Tenemos un modelo de variable latente $p(x, z; \theta)$, donde $z$ es la variable latente (que, por simplicidad, se asume que toma un número finito de valores). La densidad marginal de $x$ se puede obtener marginalizando sobre la variable latente $z$:

\[ p(x; \theta) = \sum_z p(x, z; \theta). \]

El objetivo es ajustar los parámetros $\theta$ de tal manera que se maximice la verosimilitud de los datos observados. Para esto, definimos la función de verosimilitud como:

\[ \begin{align*} \ell(\theta) &= \sum_{i=1}^n \log p(x^{(i)}; \theta). \end{align*} \]

Esto implica maximizar la probabilidad conjunta de los datos observados. Sin embargo, el cálculo directo de esta función puede ser complicado debido a la marginalización sobre las variables latentes.

Podemos reescribir la función objetivo en términos de la densidad conjunta $p(x, z; \theta)$:

\[ \ell(\theta) = \sum_{i=1}^n \log \sum_{z^{(i)}} p(x^{(i)}, z^{(i)}; \theta). \]

La suma interna hace que esta expresión sea computacionalmente costosa. Esto se debe a que la maximización directa de la verosimilitud en presencia de variables latentes puede conducir a problemas de optimización no convexa y difícil de resolver.

Motivación para el Algoritmo EM

En muchos casos, si las variables latentes $z^{(i)}$ fueran conocidas, el problema de estimación de máxima verosimilitud sería mucho más sencillo. En este contexto, el algoritmo EM (Expectation-Maximization) proporciona una solución eficiente al iterar entre dos pasos principales:

Paso E (Expectation): Calcular las expectativas de las variables latentes dado el conjunto de datos y los parámetros actuales.
Paso M (Maximization): Actualizar los parámetros maximizando la verosimilitud condicionada en las expectativas calculadas.

Idea general EM (ver las imagines)

Construcción de una cota inferior con el algoritmo EM

El algoritmo EM, en el paso E, construye una cota inferior (representada aquí por una curva verde) para la log-verosimilitud.

Propiedades de la cota inferior

Cota inferior: En todos los valores de (), la curva verde está por debajo de la curva azul, que representa la log-verosimilitud. Por lo tanto, esta curva verde es efectivamente una cota inferior.
Igualdad en el valor actual de (): En el valor actual de (), la curva verde es igual a la curva azul.

Estas propiedades son fundamentales para entender cómo funciona el paso E.

Representación del paso E

El paso E construye esta cota inferior (la curva verde), y esta construcción utiliza un adendum de la desigualdad de Jensen. En particular, bajo ciertas condiciones, la desigualdad de Jensen se convierte en una igualdad, es decir:

\[\mathbb{E}[f(x)] = f(\mathbb{E}[x]).\]

Usamos este principio para garantizar que la curva verde sea igual a la curva azul en el valor actual de (). El paso E dibuja esta curva verde.

Representación del paso M

Por otro lado, el paso M toma la curva verde y encuentra su máximo. Esto se traduce en mover el valor de () desde su posición actual (el valor inicial representado en verde) a un nuevo valor (representado en rojo).

Iteraciones del algoritmo EM

Primera iteración: El paso E construye la curva verde, y el paso M encuentra su máximo, moviendo () del valor verde al valor rojo.
Segunda iteración: En el nuevo valor de () (el valor rojo), el paso E construye una nueva cota inferior, una nueva curva verde. Luego, el paso M maximiza esta nueva curva verde, moviendo () nuevamente.
Iteraciones subsiguientes: Este proceso se repite, con el algoritmo EM construyendo y maximizando nuevas cotas inferiores, aumentando constantemente el valor de (L()) (la log-verosimilitud) hasta que converge a un óptimo local.

Convergencia del algoritmo EM

El algoritmo EM siempre converge a un óptimo local, pero no necesariamente a un óptimo global. Por ejemplo, si existiera un óptimo mejor en otra región del espacio de parámetros, el algoritmo podría no encontrarlo si no comienza cerca de esa región. Sin embargo, al repetir este proceso iterativo, el algoritmo EM converge, en general, a un buen óptimo local.

Idea EM Idea EM Idea EM

Reformulación del Problema

Para simplificar, consideremos primero la optimización de la verosimilitud $\log p(x; \theta)$ para un solo ejemplo $x$. Posteriormente, generalizaremos el resultado para $n$ ejemplos.

La función de log-verosimilitud para un ejemplo $x$ se puede reescribir como:

\[ \log p(x; \theta) = \log \sum_z p(x, z; \theta). \]

Aquí introducimos una distribución auxiliar $Q(z)$ sobre los valores posibles de $z$. Esta distribución cumple que $\sum_z Q(z) = 1$ y $Q(z) \geq 0$. Utilizando $Q(z)$, podemos reescribir la ecuación como:

\[ \begin{align*} \log p(x; \theta) &= \log \sum_z Q(z) \frac{p(x, z; \theta)}{Q(z)} \\ & \geq \sum_z Q(z) \log \frac{p(x, z; \theta)}{Q(z)}. \end{align*} \]

Esta desigualdad proviene de la desigualdad de Jensen, que establece que la función $f(x) = \log x$ es cóncava.

Prueba: consideremos el término:

\[ \sum_z Q(z) \log \frac{p(x, z; \theta)}{Q(z)} \] donde la suma representa una expectativa de la cantidad $\frac{p(x, z; \theta)}{Q(z)}$, calculada con respecto a $z$ según la distribución definida por $Q(z)$. Usando la desigualdad de Jensen, tenemos:

\[f\left( \mathbb{E}_{z \sim Q} \left[ \frac{p(x, z; \theta)}{Q(z)} \right] \right) \geq \mathbb{E}_{z \sim Q} \left[ f\left( \frac{p(x, z; \theta)}{Q(z)} \right) \right], \]

En este caso, los subíndices $z \sim Q$ indican que las expectativas se calculan respecto a $z$ distribuidos según $Q$.

Elección de $Q(z)$

Para cualquier distribución $Q$, la ecuación de arriba proporciona una cota inferior para $\log p(x; \theta)$. Sin embargo, hay muchas elecciones posibles para $Q$. ¿Cuál deberíamos elegir? Si tenemos una estimación inicial de $\theta$, es natural intentar ajustar la cota inferior de manera que sea más ajustada para ese valor particular de $\theta$. En otras palabras, buscamos que la desigualdad sea una igualdad para nuestro valor particular de $\theta$.

Recordemos esta figura: Idea EM

Para que esto sea cierto, sabemos que es suficiente que la esperanza sea tomada con respecto a una variable aleatoria “constante”. Es decir, necesitamos que:

\[ \frac{p(x,z;\theta)}{Q(z)}= c\]

para alguna constante $c$ que no dependa de $z$. Esto se logra fácilmente eligiendo (recordar que $\sum_z Q(z) = 1$ por ser una distribución):

\[ \begin{align*} Q(z) &= \frac{p(x, z; \theta)}{\sum_{z'} p(x, z'; \theta)} \\ &= p(z|x; \theta). \end{align*} \]

Así, simplemente establecemos que los $Q$ sean la distribución posterior de las $z$ dado $x$ y los parámetros $\theta$.

De hecho, podemos verificar directamente que cuando $Q(z) = p(z|x; \theta)$, entonces tenemos:

\[ \begin{align*} \sum_z Q(z) \log \frac{p(x, z; \theta)}{Q(z)} &= \sum_z p(z|x; \theta) \log \frac{p(x, z; \theta)}{p(z|x; \theta)} \\& = \sum_z p(z|x; \theta) \log \frac{p(x, z; \theta) p(z|x; \theta)}{p(z|x; \theta)} \\ &= \sum_z p(z|x; \theta) \log p(x; \theta) \\ &= \log p(x; \theta) \sum_z p(z|x; \theta) \\ &= \log p(x; \theta) \quad (\text{porque } \sum_z p(z|x; \theta) = 1). \end{align*} \]

Con esta elección, podemos definir la cota inferior de la evidencia (ELBO, por sus siglas en inglés) como:

\[ \text{ELBO}(Q, \theta) = \sum_z Q(z) \log \frac{p(x, z; \theta)}{Q(z)}. \tag{11.9} \]

Con esta ecuación, podemos reescribir:

\[ \forall Q, \theta, x, \quad \log p(x; \theta) \geq \text{ELBO}(x; Q, \theta). \tag{11.10} \]

Intuitivamente, el algoritmo EM alterna entre actualizar $Q$ y $\theta$ mediante:

Establecer $Q(z) = p(z|x; \theta)$ según la Ecuación (11.8), de forma que $\text{ELBO}(x; Q, \theta) = \log p(x; \theta)$ para un $x$ y el $\theta$ actual.
Maximizar $\text{ELBO}(x; Q, \theta)$ con respecto a $\theta$, manteniendo $Q$ fijo.

OPINO: no seria mejor escribirlo luego como algoritmo que quede claro que $$ es cada vez? ya pongo donde…

Generalización a Múltiples Ejemplos

Hasta ahora, toda la discusión se centraba en optimizar $\log p(x; \theta)$ para un solo ejemplo $x$. Sin embargo, cuando tenemos múltiples ejemplos de entrenamiento, la idea básica sigue siendo la misma. Sólo necesitamos sumar los ejemplos en los lugares relevantes. A continuación, construiremos la cota inferior de la evidencia para múltiples ejemplos de entrenamiento y formalizaremos el algoritmo EM.

Supongamos que tenemos un conjunto de entrenamiento $\{x^{(1)}, \ldots, x^{(n)}\}$. Observemos que la elección óptima de $Q$ es $p(z|x; \theta)$ y depende del ejemplo particular $x^{(i)}$. Por lo tanto, introduciremos $n$ distribuciones $Q_1, \ldots, Q_n$, una para cada ejemplo $x^{(i)}$. Para cada ejemplo $x^{(i)}$, podemos construir la cota inferior de la evidencia como:

\[ \log p(x^{(i)}; \theta) \geq \text{ELBO}(x^{(i)}; Q_i, \theta) = \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \theta)}{Q_i(z^{(i)})}. \]

Construcción de la Cota Inferior

Sumando sobre todos los ejemplos, obtenemos una cota inferior para la verosimilitud total:

\[ \begin{align*} \ell(\theta) &\geq \sum_i \text{ELBO}(x^{(i)}; Q_i, \theta) \\&= \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \theta)}{Q_i(z^{(i)})}. \end{align*} \]

Para cualquier conjunto de distribuciones $Q_1, \ldots, Q_n$, la Ecuación (11.11) proporciona una cota inferior para $\ell(\theta)$. De manera análoga al argumento de la Ecuación (11.8), la $Q_i$ que satisface la igualdad es:

\[ Q_i(z^{(i)}) = p(z^{(i)}|x^{(i)}; \theta). \]

Por lo tanto, simplemente establecemos que los $Q_i$ sean la distribución posterior de las $z^{(i)}$ dado $x^{(i)}$ con los parámetros actuales $\theta$. Esto asegura que la cota inferior se convierta en una igualdad, lo que permite que el algoritmo EM optimice de manera eficiente la verosimilitud total.

Implementación del Algoritmo EM

Ahora, para esta elección de $Q_i$, la Ecuación (11.11) proporciona una cota inferior para la verosimilitud logarítmica $\ell(\theta)$ que estamos tratando de maximizar. Este es el Paso E. En el Paso M del algoritmo, maximizamos la Ecuación (11.11) con respecto a los parámetros para obtener un nuevo ajuste de $\theta$. Repetir estos dos pasos nos da el algoritmo EM, que se describe a continuación:

Repetir hasta la convergencia:

Paso E: Para cada $i$, establecer:

\[ Q_i(z^{(i)}) := p(z^{(i)}|x^{(i)}; \theta). \]

Paso M: Maximizar:

\[ \theta := \arg\max_{\theta} \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log p(x^{(i)}, z^{(i)}; \theta). \]

Esta cota inferior nos permite maximizar la verosimilitud de manera eficiente mediante un enfoque iterativo.

##LO ESCRIBIRIA ASI:

Repetir hasta la convergencia:

Paso E: Para cada $i$, y el actual valor de $\theta_t$ establecer:

\[ Q_i(z^{(i)}) := p(z^{(i)}|x^{(i)}; \theta_t). \]

Paso M: Maximizar:

\[ \theta_{t+1} := \arg\max_{\theta} \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log p(x^{(i)}, z^{(i)}; \theta_t). \]

Esta cota inferior nos permite maximizar la verosimilitud de manera eficiente mediante un enfoque iterativo.

Reescritura con ELBO

Con esta ecuación, podemos reescribir la Ecuación (11.7) como:

\[ \forall Q, \theta, x, \quad \log p(x; \theta) \geq \text{ELBO}(x; Q, \theta). \tag{11.10} \]

Intuitivamente, el algoritmo EM alterna entre actualizar $Q$ y $\theta$ mediante:

Establecer $Q(z) = p(z|x; \theta)$ según la Ecuación (11.8), de forma que $\text{ELBO}(x; Q, \theta) = \log p(x; \theta)$ para un $x$ y el $\theta$ actual.
Maximizar $\text{ELBO}(x; Q, \theta)$ con respecto a $\theta$, manteniendo $Q$ fijo.

Generalización a Múltiples Ejemplos

\[ \log p(x^{(i)}; \theta) \geq \text{ELBO}(x^{(i)}; Q_i, \theta) = \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \theta)}{Q_i(z^{(i)})}. \]

Construcción de la Cota Inferior

Sumando sobre todos los ejemplos, obtenemos una cota inferior para la verosimilitud total:

\[ \ell(\theta) \geq \sum_i \text{ELBO}(x^{(i)}; Q_i, \theta) = \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \theta)}{Q_i(z^{(i)})}. \tag{11.11} \]

\[ Q_i(z^{(i)}) = p(z^{(i)}|x^{(i)}; \theta). \]

Implementación del Algoritmo EM

Repetir hasta la convergencia:

Paso E: Para cada $i$, establecer:

\[ Q_i(z^{(i)}) := p(z^{(i)}|x^{(i)}; \theta). \]

Paso M: Maximizar:

\[ \theta := \arg\max_{\theta} \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log p(x^{(i)}, z^{(i)}; \theta). \]

Convergencia del Algoritmo EM

¿Cómo sabemos que este algoritmo convergerá? Supongamos que $\theta^{(t)}$ y $\theta^{(t+1)}$ son los parámetros de dos iteraciones sucesivas del EM. Ahora probaremos que $\ell(\theta^{(t)}) \leq \ell(\theta^{(t+1)})$, lo que demuestra que EM siempre mejora monótonamente la verosimilitud logarítmica. La clave para mostrar este resultado radica en nuestra elección de los $Q_i$.

Específicamente, en la iteración del EM en la que los parámetros comenzaron como $\theta^{(t)}$, habríamos elegido $Q_i^{(t)}(z^{(i)}) := p(z^{(i)}|x^{(i)}; \theta^{(t)})$. Como vimos anteriormente, esta elección asegura que la desigualdad de Jensen, aplicada para obtener la Ecuación (11.11), se cumpla con igualdad, y por lo tanto:

\[ \ell(\theta^{(t)}) = \sum_{i=1}^n \text{ELBO}(x^{(i)}; Q_i^{(t)}, \theta^{(t)}). \tag{11.13} \]

Los parámetros $\theta^{(t+1)}$ se obtienen maximizando el lado derecho de la ecuación anterior. Por lo tanto:

\[ \begin{align*} \ell(\theta^{(t+1)})& \geq \sum_{i=1}^n \text{ELBO}(x^{(i)}; Q_i^{(t)}, \theta^{(t+1)})\\ & \geq \sum_{i=1}^n \text{ELBO}(x^{(i)}; Q_i^{(t)}, \theta^{(t)}) \\ &= \ell(\theta^{(t)}). \end{align*} \]

Donde la última desigualdad se deriva del hecho de que $\theta^{(t+1)}$ se elige explícitamente como:

\[ \theta := \arg\max_{\theta} \sum_{i=1}^n \text{ELBO}(x^{(i)}; Q_i^{(t)}, \theta). \]

Por lo tanto, EM hace que la verosimilitud converja de manera monótona. En nuestra descripción del algoritmo EM, dijimos que lo ejecutaríamos hasta la convergencia. Dados los resultados que acabamos de mostrar, una prueba razonable de convergencia sería verificar si el incremento en $\ell(\theta)$ entre iteraciones sucesivas es menor que algún parámetro de tolerancia, y declarar la convergencia si EM mejora $\ell(\theta)$ demasiado lentamente.

Observación sobre ELBO

Si definimos (sobrecargando ELBO):

\[ \text{ELBO}(Q, \theta) = \sum_{i=1}^n \text{ELBO}(x^{(i)}; Q_i, \theta) = \sum_i \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \theta)}{Q_i(z^{(i)})}, \]

entonces sabemos que $\ell(\theta) \geq \text{ELBO}(Q, \theta)$ por nuestra derivación anterior. El EM también puede verse como un algoritmo de maximización alternada sobre $\text{ELBO}(Q, \theta)$, donde el Paso E lo maximiza con respecto a $Q$ (verifícalo por ti mismo) y el Paso M lo maximiza con respecto a $\theta$.

Propiedades del Algoritmo EM

Convergencia Monótona: El algoritmo garantiza que la verosimilitud logarítmica $\ell(\theta)$ aumenta en cada iteración.
Criterio de Parada: Se puede detener cuando el incremento en $\ell(\theta)$ entre iteraciones sea menor que un umbral predefinido.

Otras Interpretaciones de ELBO

Sea \[ \text{ELBO}(x; Q, \theta) = \sum_z Q(z) \log \frac{p(x,z;\theta)}{Q(z)}. \] Existen varias otras formas de ELBO. Primero, podemos reescribirlo como:

\[ \begin{align*} \text{ELBO}(x; Q, \theta) &= \mathbb{E}_{z\sim Q}[\log p(x, z; \theta)] - \mathbb{E}_{z\sim Q}[\log Q(z)]. \\ &= \mathbb{E}_{z\sim Q}[\log p(x \vert z; \theta)] - D_{\text{KL}}(Q \| p_z), \end{align*} \]

donde usamos $p_z $ para denotar la distribución marginal de $z $ (bajo la distribución $p(x, z; ) $), y $D_{} $ denota la divergencia de Kullback-Leibler:

\[ D_{\text{KL}}(Q \| p_z) = \sum_z Q(z) \log \frac{Q(z)}{p(z)}. \]

En muchos casos, la distribución marginal de $z$ no depende del parámetro $\theta$. En este caso, podemos ver que maximizar respecto a $\theta$ equivale a maximizar el primer término. Esto corresponde a maximizar la verosimilitud condicional de $x$ dado $z$, lo que a menudo es una pregunta más simple que la pregunta original.

Otra forma de {ELBO} es:

\[ \text{ELBO}(x; Q, \theta) = \log p(x) - D_{\text{KL}}(Q \| p_{z \vert x}), \]

donde $p_{z \vert x}$ es la distribución condicional de $z $ dado $x$ bajo el parámetro $\theta$. Esta forma muestra que el maximizador de $\text{ELBO}(Q, \theta)$ sobre $Q $ se obtiene cuando $Q = p_{z \vert x}$, lo cual se mostró en la ecuación (11.8) anteriormente.

Mezclas de Gaussianas revisadas

Armados con nuestra definición general del algoritmo EM, volvamos a nuestro ejemplo anterior de ajuste de los parámetros $\phi$, $\mu$ y $\Sigma$ en una mezcla de Gaussianas.

Para simplificar, llevaremos a cabo las derivaciones para las actualizaciones del paso M solo para $\phi$ y $\mu_j$, y dejaremos las actualizaciones para $\Sigma_j$ como un ejercicio para el lector.

Paso E (Expectation)

El paso E es sencillo. Siguiendo nuestra derivación anterior del algoritmo, simplemente calculamos:

\[ w_j^{(i)} = Q_i(z^{(i)} = j) = P(z^{(i)} = j|x^{(i)}; \phi, \mu, \Sigma). \]

Aquí, “$Q_i(z^{(i)} = j)$” denota la probabilidad de que $z^{(i)}$ tome el valor $j$ bajo la distribución $Q_i$.

Paso M (Maximization)

En el paso M, necesitamos maximizar, con respecto a nuestros parámetros $\phi$, $\mu$, $\Sigma$, la cantidad:

\[ \sum_{i=1}^n \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \phi, \mu, \Sigma)}{Q_i(z^{(i)})}. \]

Reescribiendo esto, tenemos:

\[ \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log \frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}} \exp\left(-\frac{1}{2}(x^{(i)} - \mu_j)^T \Sigma_j^{-1} (x^{(i)} - \mu_j)\right) \cdot \phi_j/ w_j^{(i)}. \]

La expresión se puede descomponer en términos independientes, uno para cada parámetro que deseamos actualizar.

Actualización de $\mu_j$

Maximicemos con respecto a $\mu$. Tomando la derivada con respecto a $\mu$, obtenemos:

\[ \nabla_{\mu_j} \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log \phi_j \exp\left(-\frac{1}{2}(x^{(i)} - \mu_j)^T \Sigma_j^{-1}(x^{(i)} - \mu_j)\right)/ w_j^{(i)} = 0. \]

Resolviendo, encontramos: (puede dejarse como ejercicio)

\[ \mu_j = \frac{\sum_{i=1}^n w_j^{(i)} x^{(i)}}{\sum_{i=1}^n w_j^{(i)}}. \]

Actualización de $\phi_j$

Para $\phi_j$, debemos maximizar:

\[ \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log \phi_j. \]

Con la restricción de que $\sum_{j=1}^k \phi_j = 1$, construimos el Lagrangiano:

\[ \mathcal{L}(\phi) = \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log \phi_j + \beta \left(\sum_{j=1}^k \phi_j - 1\right), \]

donde $\beta$ es el multiplicador de Lagrange. Resolviendo, y usando la restriccion, obtenemos:

\[ \phi_j = \frac{1}{n} \sum_{i=1}^n w_j^{(i)}. \]

La derivación para las actualizaciones del paso M de $\Sigma_j$ también es directa.

Ejemplo Simulado

Ejemplo con Datos Reales

Usaremos el conjunto de datos de iris para ilustrar la aplicación del algoritmo EM en un caso real.

Factor análisis

Introducción

Cuando tenemos datos $x^{(i)} \in \mathbb{R}^d$ que provienen de una mezcla de varias distribuciones Gaussianas, el algoritmo EM se puede aplicar para ajustar un modelo de mezcla. En este contexto, generalmente imaginamos problemas donde disponemos de suficientes datos como para discernir la estructura de múltiples Gaussianas en los datos. Por ejemplo, este sería el caso si el tamaño de nuestro conjunto de entrenamiento $n$ fuera significativamente mayor que la dimensión $p$ de los datos.

Ahora, consideremos un escenario en el que $p > n$. En tal problema, podría ser difícil modelar los datos incluso con una sola Gaussiana, mucho menos con una mezcla de Gaussianas. Específicamente, dado que los $n$ puntos de datos solo abarcan un subespacio de baja dimensión de $\mathbb{R}^p$, si modelamos los datos como una Gaussiana y estimamos la media y la covarianza utilizando los estimadores de máxima verosimilitud habituales, tenemos:

\[ \mu = \frac{1}{n} \sum_{i=1}^{n} x^{(i)} \]

\[ \Sigma = \frac{1}{n} \sum_{i=1}^{n} (x^{(i)} - \mu)(x^{(i)} - \mu)^\top, \]

encontraríamos que la matriz $\Sigma$ es singular. Esto significa que $\Sigma^{-1}$ no existe, y que $1/|\Sigma|^{1/2} = 1/0$. Pero ambos términos son necesarios para calcular la densidad habitual de una distribución Gaussiana multivariada. Otra manera de expresar esta dificultad es que las estimaciones de máxima verosimilitud de los parámetros resultan en una Gaussiana que coloca toda su probabilidad en el espacio afín abarcado por los datos, y esto corresponde a una matriz de covarianza singular.

Más generalmente, a menos que $n$ exceda a $p$ por una cantidad razonable, las estimaciones de máxima verosimilitud de la media y la covarianza pueden ser bastante deficientes. No obstante, aún nos gustaría poder ajustar un modelo Gaussiano razonable a los datos y, tal vez, capturar alguna estructura de covarianza interesante en los datos. ¿Cómo podemos lograr esto?

En la siguiente sección, comenzamos revisando dos posibles restricciones en $\Sigma$ que nos permiten ajustarla con pequeñas cantidades de datos, aunque ninguna proporcionará una solución satisfactoria a nuestro problema. Luego, discutimos algunas propiedades de las Gaussianas que serán necesarias más adelante; específicamente, cómo encontrar distribuciones marginales y condicionales de Gaussianas. Finalmente, presentamos el modelo de análisis factorial y el algoritmo EM para este modelo.

Restricciones de $\Sigma$

Si no tenemos suficientes datos para ajustar una matriz de covarianza completa, podemos imponer algunas restricciones en el espacio de matrices $\Sigma$ que consideraremos. Por ejemplo, podríamos optar por ajustar una matriz de covarianza $\Sigma$ que sea diagonal. En este caso, el lector puede verificar fácilmente que la estimación de máxima verosimilitud de la matriz de covarianza está dada por la matriz diagonal $\Sigma$ que satisface:

\[ \Sigma_{jj} = \frac{1}{n} \sum_{i=1}^n (x_j^{(i)} - \mu_j)^2. \]

Así, $\Sigma_{jj}$ es simplemente la estimación empírica de la varianza de la $j$-ésima coordenada de los datos.

Recuerda que los contornos de una densidad Gaussiana son elipses. Una $\Sigma$ diagonal corresponde a una Gaussiana donde los ejes principales de estas elipses están alineados con los ejes.

A veces, podemos imponer una restricción adicional a la matriz de covarianza: no solo debe ser diagonal, sino que todas sus entradas diagonales deben ser iguales. En este caso, tenemos $\Sigma = \sigma^2 I$, donde $\sigma^2$ es el parámetro bajo nuestro control. La estimación de máxima verosimilitud de $\sigma^2$ se encuentra como:

\[ \sigma^2 = \frac{1}{nd} \sum_{j=1}^d \sum_{i=1}^n (x_j^{(i)} - \mu_j)^2. \]

Este modelo corresponde a usar distribuciones Gaussianas cuyas densidades tienen contornos que son círculos (en dos dimensiones) o esferas/hiperesferas en dimensiones superiores.

Si estamos ajustando una matriz de covarianza $\Sigma$ completa y sin restricciones a los datos, es necesario que $n \geq p + 1$ para que la estimación de máxima verosimilitud de $\Sigma$ no sea singular. Bajo cualquiera de las dos restricciones anteriores, podemos obtener una $\Sigma$ no singular cuando $n \geq 2$.

Sin embargo, restringir $\Sigma$ a ser diagonal también implica modelar las diferentes coordenadas $x_i$, $x_j$ de los datos como no correlacionadas e independientes. A menudo, sería deseable capturar alguna estructura de correlación interesante en los datos. Si utilizáramos cualquiera de las restricciones descritas anteriormente sobre $\Sigma$, no podríamos hacerlo. En este conjunto de notas, describiremos el modelo de análisis factorial, que utiliza más parámetros que la $\Sigma$ diagonal y captura algunas correlaciones en los datos, pero sin necesidad de ajustar una matriz de covarianza completa.

Marginales y Condicionales de Gaussianas (esto no va aca sino con la normal pero bueno, para que ya este)

Antes de describir el análisis factorial, nos detenemos brevemente para hablar sobre cómo encontrar las distribuciones marginales y condicionales de variables aleatorias con una distribución conjunta Gaussiana multivariada.

Supongamos que tenemos una variable aleatoria vectorial:

\[ x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \]

donde $x_1 \in \mathbb{R}^r$, $x_2 \in \mathbb{R}^s$, y $x \in \mathbb{R}^{r+s}$. Supongamos que $x \sim \mathcal{N}(\mu, \Sigma)$, donde:

\[ \mu = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \quad \Sigma = \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}. \]

Aquí, $\mu_1 \in \mathbb{R}^r$, $\mu_2 \in \mathbb{R}^s$, $\Sigma_{11} \in \mathbb{R}^{r \times r}$, $\Sigma_{22} \in \mathbb{R}^{s \times s}$, y así sucesivamente. Nota que las matrices de covarianza son simétricas, lo que implica que $\Sigma_{12} = \Sigma_{21}^\top$.

Bajo nuestras suposiciones, $x_1$ y $x_2$ son conjuntamente Gaussianas multivariadas. ¿Cuál es la distribución marginal de $x_1$? No es difícil ver que $\mathbb{E}[x_1] = \mu_1$ y que $\text{Cov}[x_1] = \Sigma_{11}$. Para ver que esto último es cierto, observemos que, por definición de la covarianza conjunta de $x_1$ y $x_2$, tenemos:

\[ \text{Cov}[x] = \mathbb{E}\left[ \begin{bmatrix} x_1 - \mu_1 \\ x_2 - \mu_2 \end{bmatrix} \begin{bmatrix} x_1 - \mu_1 \\ x_2 - \mu_2 \end{bmatrix}^\top \right]. \]

Expandiendo y emparejando los términos correspondientes en las matrices, podemos demostrar que la distribución marginal de $x_1$ está dada por $x_1 \sim \mathcal{N}(\mu_1, \Sigma_{11})$.

Ahora, ¿cuál es la distribución condicional de $x_1$ dado $x_2$? Usando el hecho de que la distribución conjunta de $x_1$ y $x_2$ es Gaussiana multivariada, podemos derivar que la distribución condicional de $x_1$ dado $x_2$ está dada por:

\[ x_1 \mid x_2 \sim \mathcal{N}(\mu_{1 \mid 2}, \Sigma_{1 \mid 2}), \]

donde:

\[ \mu_{1 \mid 2} = \mu_1 + \Sigma_{12} \Sigma_{22}^{-1} (x_2 - \mu_2), \]

\[ \Sigma_{1 \mid 2} = \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21}. \]

Cuando trabajamos con el modelo de análisis factorial en la próxima sección, estas fórmulas para encontrar distribuciones marginales y condicionales de Gaussianas serán muy útiles.

Modelo de Análisis de Factores

En este modelo, los datos $x$ son generados por una combinación lineal de variables latentes $z$ y ruido Gaussiano:

\[ \begin{aligned} &z \sim \mathcal{N}(0, I), \\ &x \vert z \sim \mathcal{N}(\mu + \Lambda z, \Psi). \end{aligned} \]

Donde:

$\mu$: Vector de medias de las variables observadas (dimensión $\mathbb{R}^p$).
$\Lambda$: Matriz de cargas factoriales, que define la relación lineal entre las variables latentes y las observadas (dimensión $\mathbb{R}^{p \times k}$, con $k < p$).
$\Psi$: Matriz diagonal que representa la varianza específica o ruido asociado a cada variable observada.

El modelo supone que las variables latentes $z$ tienen una distribución normal estándar y que el ruido Gaussiano es independiente y tiene una estructura diagonal.

Imaginemos que cada punto de datos $x^{(i)}$ es generado al muestrear una Gaussiana multivariada de dimensión $k$, denotada como $z^{(i)}$. Luego, $z^{(i)}$ es mapeado a un espacio afín de $\mathbb{R}^d$ al computar $\mu + \Lambda z^{(i)}$. Finalmente, $x^{(i)}$ es generado al agregar ruido con covarianza $\Psi$ a $\mu + \Lambda z^{(i)}$.

De manera equivalente (puedes verificar que este es el caso), también podemos definir el modelo de análisis factorial como:

\[ z \sim \mathcal{N}(0, I) \]

\[ \epsilon \sim \mathcal{N}(0, \Psi) \]

\[ x = \mu + \Lambda z + \epsilon, \]

donde $\epsilon$ y $z$ son independientes.

Ahora, calculemos exactamente qué distribución define nuestro modelo. Las variables aleatorias $z$ y $x$ tienen una distribución conjunta Gaussiana:

\[ \begin{bmatrix} z \\ x \end{bmatrix} \sim \mathcal{N}(\mu_{zx}, \Sigma). \]

Vamos a encontrar $\mu_{zx}$ y $\Sigma$.

Sabemos que $\mathbb{E}[z] = 0$, debido al hecho de que $z \sim \mathcal{N}(0, I)$. Además, tenemos que:

\[ \begin{align*} \mathbb{E}[x] &= \mathbb{E}[\mu + \Lambda z + \epsilon]\\ &= \mu + \Lambda \mathbb{E}[z] + \mathbb{E}[\epsilon]\\ &= \mu. \end{align*} \]

Combinando estos resultados, obtenemos:

\[ \mu_{zx} = \begin{bmatrix} 0 \\ \mu \end{bmatrix}. \]

A continuación, para encontrar $\Sigma$, necesitamos calcular los siguientes bloques:

$\Sigma_{zz} = \mathbb{E}[(z - \mathbb{E}[z])(z - \mathbb{E}[z])^\top]$ (el bloque superior izquierdo de $\Sigma$),
$\Sigma_{zx} = \mathbb{E}[(z - \mathbb{E}[z])(x - \mathbb{E}[x])^\top]$ (el bloque superior derecho de $\Sigma$),
$\Sigma_{xx} = \mathbb{E}[(x - \mathbb{E}[x])(x - \mathbb{E}[x])^\top]$ (el bloque inferior derecho de $\Sigma$).

Dado que $z \sim \mathcal{N}(0, I)$, encontramos fácilmente que $\Sigma_{zz} = \text{Cov}(z) = I$. Además:

\[ \begin{align*} \mathbb{E}[(z - \mathbb{E}[z])(x - \mathbb{E}[x])^\top] &= \mathbb{E}[z(\mu + \Lambda z + \epsilon - \mu)^\top]\\ &= \mathbb{E}[zz^\top \Lambda^\top + z \epsilon^\top]\\ &= \Lambda^\top. \end{align*} \]

En el último paso, usamos el hecho de que $\mathbb{E}[zz^\top] = \text{Cov}(z)$ (ya que $z$ tiene media cero), y $\mathbb{E}[z\epsilon^\top] = \mathbb{E}[z]\mathbb{E}[\epsilon^\top] = 0$ (ya que $z$ y $\epsilon$ son independientes).

De manera similar, podemos encontrar $\Sigma_{xx}$ de la siguiente manera:

\[ \begin{align*} \mathbb{E}[(x - \mathbb{E}[x])(x - \mathbb{E}[x])^\top] &= \mathbb{E}[(\mu + \Lambda z + \epsilon - \mu)(\mu + \Lambda z + \epsilon - \mu)^\top] \\ &= \mathbb{E}[\Lambda zz^\top \Lambda^\top + \epsilon z^\top \Lambda^\top + \Lambda z \epsilon^\top + \epsilon \epsilon^\top] \\ &= \Lambda \mathbb{E}[zz^\top] \Lambda^\top + \mathbb{E}[\epsilon \epsilon^\top]\\ &= \Lambda \Lambda^\top + \Psi. \end{align*} \]

Combinando todo, obtenemos lo siguiente:

\[ \begin{bmatrix} z \\ x \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mathbf{0} \\ \mu \end{bmatrix}, \begin{bmatrix} I & \Lambda^\top \\ \Lambda & \Lambda \Lambda^\top + \Psi \end{bmatrix} \right). \]

Distribución Marginal: La distribución marginal de las variables observadas $x$ es una distribución Gaussiana con media $\mu$ y matriz de covarianza $\Sigma = \Lambda\Lambda^T + \Psi$:

\[ x \sim \mathcal{N}(\mu, \Lambda\Lambda^T + \Psi). \]
Distribución Condicional: La distribución de las variables latentes $z$ dado $x$también es Gaussiana, con media y covarianza dadas por:

\[ \begin{aligned} \mu_{z \vert x} &= \Lambda^T(\Lambda\Lambda^T + \Psi)^{-1}(x - \mu), \\ \Sigma_{z \vert x} &= I - \Lambda^T(\Lambda\Lambda^T + \Psi)^{-1}\Lambda. \end{aligned} \]

Estas propiedades son esenciales para la implementación del algoritmo EM, ya que podemos escribir la verosimilitud usando lo anterior: \[ \ell(\mu, \Lambda, \Psi) = \log \prod_{i=1}^n \frac{1}{(2\pi)^{d/2} |\Lambda \Lambda^\top + \Psi|^{1/2}} \exp\left(-\frac{1}{2} (x^{(i)} - \mu)^\top (\Lambda \Lambda^\top + \Psi)^{-1} (x^{(i)} - \mu)\right). \]

Para realizar la estimación de máxima verosimilitud, nos gustaría maximizar esta cantidad con respecto a los parámetros. Sin embargo, maximizar esta fórmula explícitamente es difícil (inténtalo tú mismo), y no conocemos ningún algoritmo que lo haga en forma cerrada. Por lo tanto, en su lugar, utilizaremos el algoritmo EM. En la próxima sección, derivamos EM para el análisis factorial.o

Algoritmo EM para Análisis de Factores

El algoritmo EM optimiza los parámetros del modelo iterativamente mediante dos pasos:

Paso E (Esperanza)

Síntesis: en este paso, calculamos la distribución condicional de las variables latentes $z$ dado $x$ y los parámetros actuales del modelo. Esto implica calcular:

\[ \begin{aligned} \mu_{z \vert x} &= \Lambda^T(\Lambda\Lambda^T + \Psi)^{-1}(x - \mu), \\ \Sigma_{z \vert x} &= I - \Lambda^T(\Lambda\Lambda^T + \Psi)^{-1}\Lambda. \end{aligned} \]

Estos valores representan la media y la covarianza condicionales de $z$, que se utilizan para estimar las contribuciones de las variables latentes a los datos observados.

Derivación del Paso E:*

La derivación del paso E es sencilla. Necesitamos calcular $Q_i(z^{(i)}) = p(z^{(i)} \mid x^{(i)}; \mu, \Lambda, \Psi)$. Sustituyendo la distribución dada en la Ecuación (3) en las fórmulas (1-2) utilizadas para encontrar la distribución condicional de una Gaussiana, encontramos que:

\[ z^{(i)} \mid x^{(i)}; \mu, \Lambda, \Psi \sim \mathcal{N}(\mu_{z^{(i)} \mid x^{(i)}}, \Sigma_{z^{(i)} \mid x^{(i)}}), \]

donde:

\[ \begin{align*} \mu_{z^{(i)} \mid x^{(i)}} &= \Lambda^\top (\Lambda \Lambda^\top + \Psi)^{-1} (x^{(i)} - \mu) \\ \Sigma_{z^{(i)} \mid x^{(i)}} &= I - \Lambda^\top (\Lambda \Lambda^\top + \Psi)^{-1} \Lambda. \end{align*} \]

Por lo tanto, usando estas definiciones para $\mu_{z^{(i)} \mid x^{(i)}}$ y $\Sigma_{z^{(i)} \mid x^{(i)}}$, tenemos:

\[ Q_i(z^{(i)}) = \frac{1}{(2\pi)^{k/2} |\Sigma_{z^{(i)} \mid x^{(i)}}|^{1/2}} \exp \left( -\frac{1}{2} (z^{(i)} - \mu_{z^{(i)} \mid x^{(i)}})^\top \Sigma_{z^{(i)} \mid x^{(i)}}^{-1} (z^{(i)} - \mu_{z^{(i)} \mid x^{(i)}}) \right). \]

Paso M (Maximización)

Síntesis: en este paso, actualizamos los parámetros $\mu$, $\Lambda$ y $\Psi$ maximizando la esperanza de la log-verosimilitud:

Actualización de $\mu$:

\[ \mu = \frac{1}{n} \sum_{i=1}^n x^{(i)}. \]
Actualización de $\Lambda$:

\[ \Lambda = \left( \sum_{i=1}^n (x^{(i)} - \mu)\mu_{z \vert x}^{(i)T} \right) \left( \sum_{i=1}^n \mu_{z \vert x}^{(i)}\mu_{z \vert x}^{(i)T} + \Sigma_{z \vert x}^{(i)} \right)^{-1}. \]
Actualización de $Psi$:

\[ \Psi = \text{diag}\left( \frac{1}{n} \sum_{i=1}^n \left[ x^{(i)}x^{(i)T} - \Lambda \mu_{z \vert x}^{(i)}x^{(i)T} \right] \right). \]

Este proceso se repite hasta que los parámetros convergen, indicando que se ha encontrado un óptimo local.

Derivación del paso M:

Trabajemos ahora en el paso M. Aquí, necesitamos maximizar:

\[ \sum_{i=1}^n \int_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)}, z^{(i)}; \mu, \Lambda, \Psi)}{Q_i(z^{(i)})} dz^{(i)} \]

con respecto a los parámetros $\mu$, $\Lambda$, $\Psi$. Resolveremos únicamente la optimización con respecto a $\Lambda$ y dejaremos como ejercicio para el lector las derivaciones de las actualizaciones para $\mu$ y $\Psi$.

Podemos simplificar la función a maximizar como:

\[ \begin{align*} & \sum_{i=1}^n \int_{z^{(i)}} Q_i(z^{(i)}) \big[ \log p(x^{(i)} \mid z^{(i)}; \mu, \Lambda, \Psi) + \log p(z^{(i)}) - \log Q_i(z^{(i)}) \big] dz^{(i)} \\& = \sum_{i=1}^n \mathbb{E}_{z^{(i)} \sim Q_i} \big[ \log p(x^{(i)} \mid z^{(i)}; \mu, \Lambda, \Psi) + \log p(z^{(i)}) - \log Q_i(z^{(i)}) \big] \end{align*} \]

Aquí, el subíndice “$z^{(i)} \sim Q_i$” indica que la esperanza se toma con respecto a $z^{(i)}$ muestreado de $Q_i$. En el desarrollo subsiguiente, omitiremos este subíndice cuando no haya riesgo de ambigüedad. Al eliminar términos que no dependen de los parámetros, encontramos que necesitamos maximizar:

\[ \begin{align*} &\sum_{i=1}^n \mathbb{E} \big[ \log p(x^{(i)} \mid z^{(i)}; \mu, \Lambda, \Psi) \big] \\ &=\sum_{i=1}^n \mathbb{E} \bigg[ \log \frac{1}{(2\pi)^{d/2} |\Psi|^{1/2}} \exp \left( -\frac{1}{2} (x^{(i)} - \mu - \Lambda z^{(i)})^\top \Psi^{-1} (x^{(i)} - \mu - \Lambda z^{(i)}) \right) \bigg]. \\ &=\sum_{i=1}^n \mathbb{E} \bigg[ -\frac{1}{2} \log |\Psi| - \frac{n}{2} \log (2\pi) - \frac{1}{2} (x^{(i)} - \mu - \Lambda z^{(i)})^\top \Psi^{-1} (x^{(i)} - \mu - \Lambda z^{(i)}) \bigg]. \end{align*} \]

Maximización con respecto a $\Lambda$

Maximicemos esta expresión con respecto a $\Lambda$. Solo el último término en la ecuación anterior depende de $\Lambda$. Tomando las derivadas y utilizando los hechos de que $\text{tr}(AB) = \text{tr}(BA)$ y $\nabla_{\text{tr}}(ABAC) = CAB + C^\top AB^\top$, obtenemos:

\[ \begin{align*} &\nabla_\Lambda \sum_{i=1}^n \mathbb{E} \bigg[ -\frac{1}{2} (x^{(i)} - \mu - \Lambda z^{(i)})^\top \Psi^{-1} (x^{(i)} - \mu - \Lambda z^{(i)}) \bigg]\\ &= \sum_{i=1}^n \nabla_\Lambda \mathbb{E} \bigg[ -\frac{1}{2} \text{tr} \big( z^{(i)} z^{(i)\top} \Lambda^\top \Psi^{-1} \Lambda \big) + \text{tr} \big( z^{(i)} (x^{(i)} - \mu)^\top \Psi^{-1} \Lambda^\top \big) \bigg] \\ &= \sum_{i=1}^n \nabla_\Lambda \mathbb{E} \bigg[ -\frac{1}{2} \text{tr} \big( \Lambda^\top \Psi^{-1} \Lambda z^{(i)} z^{(i)\top} \big) + \text{tr} \big( \Psi^{-1} (x^{(i)} - \mu) z^{(i)\top} \Lambda^\top \big) \bigg] \\ &= \sum_{i=1}^n \mathbb{E} \bigg[ -\Psi^{-1} \Lambda z^{(i)} z^{(i)\top} + \Psi^{-1} (x^{(i)} - \mu) z^{(i)\top} \bigg]. \end{align*} \]

Al igualar esta derivada a cero y simplificar, obtenemos:

\[ \sum_{i=1}^n \Lambda \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)} z^{(i)\top} \big] = \sum_{i=1}^n (x^{(i)} - \mu) \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)\top} \big]. \]

Resolviendo para $\Lambda$, encontramos:

\[ \Lambda = \bigg( \sum_{i=1}^n (x^{(i)} - \mu) \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)\top} \big] \bigg) \bigg( \sum_{i=1}^n \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)} z^{(i)\top} \big] \bigg)^{-1}. \tag{7} \]

Relación con las ecuaciones normales

Es interesante notar la relación entre esta ecuación y la ecuación normal que derivamos para la regresión por mínimos cuadrados ordinarios:

\[ \theta^\top = (y^\top X)(X^\top X)^{-1}. \]

La analogía aquí es que los $x$ son una función lineal de los $z$ (más ruido). Dadas las “estimaciones” para $z$ encontradas en el Paso E, ahora intentaremos estimar la linealidad desconocida $\Lambda$ que relaciona los $x$ y los $z$. No es sorprendente que obtengamos algo similar a la ecuación normal. Sin embargo, hay una diferencia importante entre este procedimiento y un algoritmo que realiza mínimos cuadrados usando solo las “mejores estimaciones” de los $z$. Veremos esta diferencia en breve.

Para completar la actualización del Paso M, trabajemos los valores de las esperanzas en la Ecuación (7). A partir de nuestra definición de $Q_i$ como Gaussiana con media $\mu_{z^{(i)} \mid x^{(i)}}$ y covarianza $\Sigma_{z^{(i)} \mid x^{(i)}}$, fácilmente encontramos:

\[ \begin{align*} \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)} \big] &= \mu_{z^{(i)} \mid x^{(i)}} \\ \mathbb{E}_{z^{(i)} \sim Q_i} \big[ z^{(i)} z^{(i)\top} \big] &= \mu_{z^{(i)} \mid x^{(i)}} \mu_{z^{(i)} \mid x^{(i)}\top} + \Sigma_{z^{(i)} \mid x^{(i)}}. \end{align*} \]

El último resultado proviene del hecho de que, para una variable aleatoria $Y$, $\text{Cov}(Y) = \mathbb{E}[YY^\top] - \mathbb{E}[Y]\mathbb{E}[Y^\top]$, y por lo tanto $\mathbb{E}[YY^\top] = \mathbb{E}[Y]\mathbb{E}[Y^\top] + \text{Cov}(Y)$. Sustituyendo esto de nuevo en la Ecuación (7), obtenemos la actualización del Paso M para $\Lambda$:

\[ \Lambda = \bigg( \sum_{i=1}^n (x^{(i)} - \mu) \mu_{z^{(i)} \mid x^{(i)}\top} \bigg) \bigg( \sum_{i=1}^n \big( \mu_{z^{(i)} \mid x^{(i)}} \mu_{z^{(i)} \mid x^{(i)}\top} + \Sigma_{z^{(i)} \mid x^{(i)}} \big) \bigg)^{-1}. \tag{8} \]

Es importante notar la presencia de $\Sigma_{z^{(i)} \mid x^{(i)}}$ en el lado derecho de esta ecuación. Esta es la covarianza en la distribución posterior $p(z^{(i)} \mid x^{(i)})$, y el Paso M debe tener en cuenta esta incertidumbre sobre $z^{(i)}$ en el posterior. Un error común al derivar el EM es asumir que, en el Paso E, necesitamos calcular únicamente la esperanza $\mathbb{E}[z]$ de la variable aleatoria latente $z$, y luego usarla en la optimización en el Paso M donde aparece $z$. Aunque esto funcionó para problemas simples como mezclas de Gaussianas, en nuestra derivación para análisis factorial necesitamos $\mathbb{E}[zz^\top]$ y $\Sigma_{z^{(i)} \mid x^{(i)}}$.

Por último, también podemos encontrar las optimizaciones del Paso M para los parámetros $\mu$ y $\Psi$. No es difícil demostrar que la primera está dada por:

\[ \mu = \frac{1}{n} \sum_{i=1}^n x^{(i)}. \]

Dado que este término no cambia a medida que los parámetros se varían (es decir, a diferencia de la actualización para $\Lambda$, el lado derecho no depende de $Q_i(z^{(i)}) = p(z^{(i)} \mid x^{(i)}; \mu, \Lambda, \Psi)$), puede calcularse una sola vez y no necesita ser actualizado continuamente durante el algoritmo. De manera similar, la matriz diagonal $\Psi$ puede encontrarse calculando:

\[ \Phi = \frac{1}{n} \sum_{i=1}^n x^{(i)} x^{(i)\top} - x^{(i)} \mu_{z^{(i)} \mid x^{(i)}}^\top \Lambda^\top - \Lambda \mu_{z^{(i)} \mid x^{(i)}} x^{(i)\top} + \Lambda \big( \mu_{z^{(i)} \mid x^{(i)}} \mu_{z^{(i)} \mid x^{(i)}\top} + \Sigma_{z^{(i)} \mid x^{(i)}} \big) \Lambda^\top, \]

y estableciendo $\Psi_{ii} = \Phi_{ii}$ (es decir, haciendo que $\Psi$ sea la matriz diagonal que contiene únicamente las entradas diagonales de $\Phi$).

Ejemplo Simulado

Generemos datos siguiendo un modelo de análisis de factores con dos dimensiones latentes:

Ejemplo del Mundo Real

Usaremos el conjunto de datos Iris para aplicar el modelo y evaluar su rendimiento:

Conclusiones

El análisis de factores es una técnica poderosa para modelar datos complejos y de alta dimensionalidad. Su implementación mediante el algoritmo EM permite:

Capturar la estructura latente de los datos.
Reducir la dimensionalidad de forma eficiente.
Manejar matrices de covarianza problemáticas en escenarios con pocos datos.

Los ejemplos muestran su aplicación en contextos simulados y reales, evaluando su rendimiento con errores de reconstrucción y visualizaciones.

Introducción

EM para Mezcla de Gaussianas

El Algoritmo EM

Procedimiento Detallado

Interpretación del Algoritmo

Propiedades del Algoritmo

Funciona? (esto habria que dar algo no se si esto esta bien)

Visualización de Clusters

Log-Verosimilitud

Validación Cruzada (Evaluación Cuantitativa)

Coherencia Interna

Evaluación Supervisada (Si Hay Etiquetas)

Robustez ante Diferentes Inicializaciones

Evaluación Manual o Contextual

Implementación del Algoritmo EM y aplicación a un ejemplo simulado y un ejemplo Real

Convergencia a un mínimo local del algoritmo EM

Desigualdad de Jensen (capaz se enseño antes, igual pongamoslo aca por las duads y luego lo pasamos)

Algoritmos Generales EM

Introducción

Motivación para el Algoritmo EM

Idea general EM (ver las imagines)

Construcción de una cota inferior con el algoritmo EM

Propiedades de la cota inferior

Representación del paso E

Representación del paso M

Iteraciones del algoritmo EM

Convergencia del algoritmo EM

Reformulación del Problema

Elección de \(Q(z)\)

Generalización a Múltiples Ejemplos

Construcción de la Cota Inferior

Implementación del Algoritmo EM

Reescritura con ELBO

Generalización a Múltiples Ejemplos

Construcción de la Cota Inferior

Implementación del Algoritmo EM

Convergencia del Algoritmo EM

Observación sobre ELBO

Propiedades del Algoritmo EM

Otras Interpretaciones de ELBO

Mezclas de Gaussianas revisadas

Paso E (Expectation)

Paso M (Maximization)

Actualización de \(\mu_j\)

Actualización de \(\phi_j\)

Ejemplo Simulado

Ejemplo con Datos Reales

Factor análisis

Introducción

Restricciones de \(\Sigma\)

Marginales y Condicionales de Gaussianas (esto no va aca sino con la normal pero bueno, para que ya este)

Modelo de Análisis de Factores

Algoritmo EM para Análisis de Factores

Paso E (Esperanza)

Paso M (Maximización)

Ejemplo Simulado

Ejemplo del Mundo Real

Conclusiones