El algoritmo EM

Qué cubre este capítulo

k-means tomaba una decisión dura por cada punto: perteneces a este cluster, y punto. Eso funcionaba porque los blobs estaban bien separados y la decisión era obvia. Pero la mayoría de los grupos reales se traslapan, y cerca de la costura entre dos de ellos una decisión dura es una adivinanza disfrazada de hecho. Este capítulo trata del algoritmo que se niega a adivinar. Le entrega a cada punto una probabilidad de pertenecer a cada grupo, ajusta los grupos usando esas probabilidades como pesos, y repite — suavizando cada decisión hasta que los números dejan de moverse.

Ese algoritmo es Expectation-Maximization, EM, y es una de las ideas que cargan el peso de toda la estadística. El planteamiento es siempre el mismo: tienes datos, y falta alguna pieza de información que necesitarías para ajustarlos limpiamente. Aquí la pieza faltante es cuál de dos gaussianas produjo cada número. Si lo supieras, ajustar las dos gaussianas sería cosa de una línea. Si conocieras las dos gaussianas, adivinar las etiquetas sería cosa de una línea. No sabes ni una cosa ni la otra, así que EM alterna: adivina las etiquetas de forma suave a partir de las gaussianas actuales (el paso E), reajusta las gaussianas a partir de esas etiquetas suaves (el paso M), y vuelve a empezar. Construimos el ciclo completo a mano en NumPy sobre el caso más limpio posible — una sola columna de números extraídos de dos campanas que se traslapan — y luego lo verificamos contra scikit-learn.

Lo que hay que vigilar es la log-verosimilitud, un solo número que califica qué tan bien los parámetros actuales explican los datos. La garantía central de EM es que ese número nunca baja. No "casi siempre sube" — nunca baja, en cada iteración, de forma demostrable. Esa subida monótona es la razón entera por la que se confía en el método, y la animación te hace verla suceder. Este capítulo trata de la receta general; la historia completa en 2-D, donde cada componente es una elipse de covarianza inclinada, tiene su propio capítulo sobre modelos de mezclas gaussianas.

Un poco de historia

El nombre y la forma general vienen de un solo paper famoso: Arthur Dempster, Nan Laird y Donald Rubin, "Maximum Likelihood from Incomplete Data via the EM Algorithm," leído ante la Royal Statistical Society en 1977. Ellos no inventaron los trucos individuales — los casos especiales llevaban años flotando por ahí, desde el trabajo de H.O. Hartley de 1958 sobre conteos hasta el procedimiento Baum-Welch para modelos ocultos de Markov a finales de los sesenta. Lo que hicieron Dempster, Laird y Rubin fue darse cuenta de que todos esos eran el mismo algoritmo. Cualquier problema que pudieras plantear como "los datos que tengo son una versión incompleta de datos que podría ajustar fácilmente" encajaba en una receta maestra, y ellos le dieron un nombre, un enunciado general y el argumento de por qué funciona.

Ese argumento tenía un hueco. El paper de 1977 afirmaba que la log-verosimilitud convergía de una manera ligeramente demasiado fuerte, y en 1983 C.F. Jeff Wu publicó la corrección, precisando exactamente qué promete EM y qué no — sube hasta un punto estacionario de la verosimilitud, que no tiene por qué ser el máximo global. Esa salvedad no es una nota al pie; es lo que sientes cada vez que corres EM en serio, y volveremos a ella. Casi cincuenta años después, el paper es uno de los más citados en la historia de la estadística, porque el patrón que nombró está en todas partes: mezclas, respuestas faltantes en encuestas, estados ocultos, temas latentes. Aprende EM una vez y lo sigues reconociendo por todos lados.

La intuición

Imagina los números sobre una línea. En algún lugar a la izquierda hay un montón, en algún lugar a la derecha hay un montón más grande, y en medio se difuminan uno con el otro de modo que no puedes decir dónde termina uno y empieza el otro. Crees que dos procesos generaron estos datos, pero la etiqueta que te diría qué proceso generó cada punto nunca se registró. Esa es la información faltante, y es toda la dificultad.

Aquí está el truco que la disuelve. Supón que alguien te susurrara al oído las gaussianas verdaderas — la de la izquierda centrada en -2, la de la derecha en 3, estas dispersiones, esta mezcla. Entonces, para cualquier punto, podrías calcular qué tan probable era que cada gaussiana lo hubiera producido, y la proporción te da una etiqueta suave: 90% izquierda, 10% derecha, digamos. Ahora supón, en cambio, que alguien te entregara esas etiquetas suaves para cada punto. Entonces ajustar las gaussianas es solo un promedio ponderado — cada punto contribuye a la media de la gaussiana izquierda en proporción a qué tanto parece de la izquierda. Cualquiera de las dos mitades es fácil dada la otra. No tienes ninguna, así que haces bootstrap: empiezas con una conjetura de las gaussianas, por mala que sea, calculas las etiquetas suaves que implican, reajustas las gaussianas a partir de esas etiquetas, y repites el ciclo.

Ese ciclo es EM. El paso E es "dados los parámetros actuales, qué tanto pertenece cada punto a cada componente". El paso M es "dadas esas membresías suaves, cuál es el componente que mejor ajusta". Cada pasada afina a la otra, y lo hermoso es que una sola cantidad — la log-verosimilitud de los datos — sube cada vez, así que el proceso no puede dar bandazos ni vagar sin rumbo. Aquí están los datos sobre los que va a trabajar, como histograma de densidad, dos campanas embarradas en una sola forma llena de bultos:

Casi puedes ver dos jorobas — una más baja cerca de -2, una más ancha cerca de 3 — pero el centro es genuinamente ambiguo, y ningún corte duro ahí es honesto. El trabajo de EM es sacar las dos campanas de vuelta de ese bulto.

Las matemáticas

Escribe los datos como $x_1, \dots, x_N$, un número cada uno. Los modelamos como una mezcla de $K$ gaussianas (aquí $K = 2$), donde el componente $k$ tiene peso de mezcla $\pi_k$ (con $\sum_k \pi_k = 1$), media $\mu_k$ y varianza $\sigma_k^2$. Reúne todo esto en $\theta$. La probabilidad que el modelo asigna a un punto es la suma ponderada sobre los componentes:

La información faltante es una etiqueta oculta $z_i \in \{1, \dots, K\}$ por punto que dice qué componente lo generó. El paso E calcula su posterior — la responsabilidad $\gamma_{ik}$, la probabilidad de que el componente $k$ haya producido el punto $i$ dados los parámetros actuales. Es la regla de Bayes: la densidad ponderada del componente sobre el total.

Cada fila $\gamma_{i\cdot}$ suma 1 — una asignación suave completa. El paso M entonces trata estas responsabilidades como membresías fraccionarias y reajusta cada componente mediante promedios ponderados. Sea $N_k = \sum_i \gamma_{ik}$ el conteo suave de puntos en el componente $k$:

Estas son las fórmulas ordinarias de máxima verosimilitud para un peso, una media y una varianza, con $\gamma_{ik}$ representando "qué tanto del punto $i$ cuenta para el componente $k$". La cantidad que ambos pasos van mejorando en silencio es la log-verosimilitud de datos incompletos — la log-probabilidad de los números solos, con las etiquetas ocultas sumadas fuera:

Por qué la subida está garantizada requiere una idea más. Para cualquier elección de etiquetas suaves $q(z_i)$, la log-verosimilitud se divide exactamente en una cota inferior más una brecha:

El primer término $\mathcal{L}(q, \theta)$ es el ELBO, la cota inferior de la evidencia; el segundo es una divergencia KL, que nunca es negativa. El paso E fija $q$ en el posterior verdadero, lo que hace cero la brecha KL y vuelve ajustada la cota. El paso M maximiza la cota sobre $\theta$. Subir una cota inferior ajustada solo puede subir lo que está encima de ella, así que:

Esa desigualdad es la promesa entera de EM. No dice nada sobre alcanzar la mejor respuesta — solo que nunca empeoras. Que es exactamente la cosa de doble filo que construimos a continuación.

En qué es bueno, en qué no

El don de EM es que vuelve tratables verosimilitudes intratables. Maximizar $\ell(\theta)$ directamente significa derivar un log-de-una-suma, y la suma dentro del log enreda todos los componentes en algo sin solución de forma cerrada. EM se salta eso por completo: el paso E y el paso M tienen cada uno formas cerradas y limpias, y alternarlos optimiza el objetivo difícil sin tocarlo nunca de frente. Obtienes asignaciones suaves gratis, que es la salida honesta cuando los grupos se traslapan, y la subida monótona significa que el ciclo es estable — sin learning rate que afinar, sin paso que pueda explotar. Cuando tu problema genuinamente tiene estructura oculta o valores faltantes, EM es muy a menudo lo primero que funciona.

El detalle es el que Wu tuvo que corregir en la demostración: EM encuentra un óptimo local, no el global, y cuál óptimo local depende por completo de dónde empiezas. Un mal init puede dejarlo varado en un mal ajuste que es perfectamente estable — la log-verosimilitud se aplana, el algoritmo reporta éxito, y la respuesta está mal. También puede hacer cosas genuinamente feas que una media nunca hace: dejar que un componente colapse sobre un solo punto, llevando su varianza hacia cero y la verosimilitud hacia infinito, una solución degenerada contra la que tienes que protegerte. Y puede ser lento, arrastrándose hacia el óptimo en pasos diminutos una vez que está cerca. Las defensas estándar son todas sobre el arranque: córrelo desde varios inits aleatorios y quédate con el mejor, o siémbralo con k-means, que es exactamente lo que hace scikit-learn.

Los datos

Una columna de 300 números. La construí muestreando de dos gaussianas y tirando las etiquetas: 40% de los puntos vienen de $\mathcal{N}(-2, 1)$, el otro 60% de $\mathcal{N}(3, 1.5)$ — centros distintos, dispersiones distintas, mezcla desigual, y suficiente traslape en el centro como para que ninguna línea limpia los separe. El generador registra de qué gaussiana vino cada punto en una columna comp, pero eso es un sobre sellado: EM solo ve los números. Abrimos el sobre una vez al final, para calificar qué tan bien el clustering suave recuperado coincide con la verdad. Dos componentes en una dimensión es la mezcla más pequeña que aun así muestra todo lo que hace EM, y por eso es el lugar correcto para observar el mecanismo en vez del paisaje.

Constrúyelo, una función a la vez

Cinco funciones. Una gaussiana, los dos pasos, el puntaje que optimizan y el ciclo que los une — ese es el algoritmo entero. Todo empieza con la densidad de una sola gaussiana 1-D, evaluada sobre todos los puntos a la vez:

def gaussian(x, mu, var):
    """1-D Gaussian density N(x | mu, var), evaluated elementwise over x.

    var is the variance (sigma squared), not the standard deviation. This is
    the single building block the whole algorithm leans on — every step below
    is just weighted sums of this number.
    """
    return np.exp(-0.5 * (x - mu) ** 2 / var) / np.sqrt(2.0 * np.pi * var)

Nota que recibe la varianza, no la desviación estándar — ese es el parámetro que el paso M realmente actualiza, así que cargarla directamente mantiene el código honesto. Encima de esa función está el paso E. Para cada punto, pondera la densidad de cada componente por su peso de mezcla y normaliza entre componentes; el resultado es la matriz de responsabilidades, una fila por punto que suma 1:

def e_step(x, weights, means, variances):
    """E-step: responsibilities gamma[i, k] = P(component k produced x_i).

    For each point we weight every component's density by its mixing weight and
    normalize across components, so each row sums to 1. gamma[i, k] near 1 means
    point i almost certainly belongs to component k; 0.5 means it's on the fence.
    Returns an (N, K) matrix — the soft assignment the M-step averages over.
    """
    weighted = np.array([w * gaussian(x, m, v)
                         for w, m, v in zip(weights, means, variances)]).T  # (N, K)
    return weighted / weighted.sum(axis=1, keepdims=True)

Eso es la regla de Bayes en tres líneas. Un punto que cae bajo la campana izquierda pone casi toda su masa en el componente 0; un punto en el centro turbio se reparte casi 50/50. El paso M toma esas membresías suaves y reajusta cada componente mediante promedios ponderados:

def m_step(x, gamma):
    """M-step: re-estimate weights, means, variances as weighted averages.

    N_k is the soft count of points in component k — the column sum of the
    responsibilities. Everything is the ordinary formula for a weighted mean and
    weighted variance, with gamma[:, k] as the weights: a point that half-belongs
    to a component contributes half of itself to that component's estimates.
    """
    Nk = gamma.sum(axis=0)                                   # (K,) soft counts
    weights = Nk / len(x)
    means = (gamma * x[:, None]).sum(axis=0) / Nk
    variances = (gamma * (x[:, None] - means) ** 2).sum(axis=0) / Nk
    return weights, means, np.maximum(variances, 1e-6)       # floor guards /0

Nk es el conteo suave — el número de puntos que posee cada componente, fraccionario porque la posesión se comparte. El peso es ese conteo sobre N; la media y la varianza son la media y la varianza ponderadas por responsabilidad. La única salvaguarda es un piso sobre la varianza, para que un componente no pueda colapsar en un pico y mandar la verosimilitud a infinito. Ahora el puntaje que ambos pasos van escalando, la log-verosimilitud de datos incompletos:

def log_likelihood(x, weights, means, variances):
    """Incomplete-data log-likelihood: sum_i log sum_k pi_k N(x_i | mu_k, var_k).

    "Incomplete" because it's the likelihood of the x's alone, with the hidden
    component labels summed out. This is the single number EM optimizes, and the
    one guarantee is that it never goes down from one iteration to the next.
    """
    per = np.array([w * gaussian(x, m, v)
                    for w, m, v in zip(weights, means, variances)]).T   # (N, K)
    return float(np.log(per.sum(axis=1)).sum())

Directo de las matemáticas: para cada punto, el log de su densidad total de mezcla; suma sobre los puntos. Este es el número que nunca debe caer. Vigilarlo es cómo sabes que el ciclo está funcionando y cómo atrapas un bug — si alguna vez baja, algo está mal en el paso E o en el paso M. Por último, el ciclo mismo:

def em(x, init_weights, init_means, init_variances, max_iter=200, tol=1e-6):
    """EM loop: alternate E-step and M-step until the log-likelihood plateaus.

    Starts from whatever parameters you hand it — EM only finds a local optimum,
    so the init matters. Each pass computes responsibilities (E), re-estimates the
    components from them (M), and checks the log-likelihood; when it stops rising
    by more than tol, we've converged. Returns the fitted parameters, the final
    log-likelihood, and the full per-iteration history the chapter replays.
    """
    weights = np.array(init_weights, float)
    means = np.array(init_means, float)
    variances = np.array(init_variances, float)

    ll = log_likelihood(x, weights, means, variances)
    history = [_snapshot(x, weights, means, variances, ll)]
    for _ in range(max_iter):
        gamma = e_step(x, weights, means, variances)            # E-step
        weights, means, variances = m_step(x, gamma)            # M-step
        new_ll = log_likelihood(x, weights, means, variances)
        history.append(_snapshot(x, weights, means, variances, new_ll))
        if new_ll - ll < tol:                                   # plateaued
            ll = new_ll
            break
        ll = new_ll
    return weights, means, variances, ll, history

Paso E, paso M, revisar el puntaje, parar cuando deja de subir más que un pelito. Registra una instantánea en cada iteración — los parámetros, la log-verosimilitud y la responsabilidad actual de cada punto — que es lo que la animación reproduce. La convergencia es solo la log-verosimilitud aplanándose, porque una secuencia monótona plana no tiene a dónde más ir.

Míralo trabajar

Este es el premio por usar una dimensión. Abajo hay una corrida real de la función em de arriba, un frame por iteración, y la arranqué a propósito desde un init deliberadamente terrible: ambas medias de los componentes amontonadas cerca del centro en -0.5 y 0.5, pesos iguales, varianza unitaria. Desde ahí te toca ver cómo EM separa las dos gaussianas a la fuerza.

Tres paneles, todos en vivo. El de arriba superpone las dos curvas actuales de los componentes — cada gaussiana escalada por su peso de mezcla — sobre el histograma fijo de los datos; observa la curva cian deslizarse a la izquierda y la naranja deslizarse a la derecha y engordar. La franja de en medio son los 300 puntos, cada uno coloreado por su responsabilidad actual: cian puro significa "con certeza el componente izquierdo", naranja puro significa "con certeza el derecho", y una mezcla lodosa en medio significa que el punto sigue indeciso. Observa cómo el centro se endurece de papilla a dos colores comprometidos. El panel de abajo traza la log-verosimilitud, el único número que solo sabe subir.

Dale play y el primer frame hace casi todo el trabajo. La log-verosimilitud salta de -1,386.5 en el mal init a -703.9 después de una sola pasada E/M — las dos medias se separan de golpe desde su arranque amontonado hacia -2 y 3 en un solo paso. De ahí en adelante es refinamiento: -689.7, luego -682.8, luego un arrastre largo y lento mientras las curvas se acomodan en el histograma y los indecisos del centro escogen un lado. Para la iteración 8 el puntaje ya es -681.5 y la imagen ha dejado de cambiar visiblemente; el ciclo corre hasta 20 iteraciones solo porque sigue exprimiendo mejoras más pequeñas de lo que puedes ver, y se detiene cuando caen por debajo de la tolerancia. Converge a una log-verosimilitud de -681.48.

Ahora mira el panel de abajo a lo largo de toda la corrida y encuentra un solo frame donde la línea baje. No hay ninguno. Eso no es suerte de esta semilla — es la garantía de las matemáticas hecha visible, y es la razón por la que puedes confiar en que EM va subiendo incluso cuando, como aquí, arranca de un ajuste francamente ridículo. Aquí está esa misma subida por sí sola, cada iteración:

Un acantilado en el primer paso, luego una rodilla, luego plano. Monótono todo el camino. Esa forma — ganancias enormes al principio, luego una cola larga de rendimientos decrecientes — es la firma de EM, y es la razón por la que la gente pone un tope de iteraciones y fija una tolerancia en vez de esperar a que los números dejen de moverse por completo.

La implementación completa

El archivo entero, sin librería, de arriba a abajo. Esto es exactamente lo que corrió la animación:

"""Expectation-Maximization for a 1-D Gaussian mixture, built from scratch.

The data is one column of numbers that came from two overlapping Gaussians, and
nobody wrote down which point came from which. EM recovers the two Gaussians —
their means, their spreads, and how much of the data each one produced — by
alternating two steps: an E-step that computes, for every point, how strongly it
belongs to each component (soft assignment), and an M-step that re-estimates the
components as responsibility-weighted averages. Pure NumPy — no ML library in
this file.

Every function below appears in the chapter one step at a time (the
`# region:` markers are what the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: gaussian
def gaussian(x, mu, var):
    """1-D Gaussian density N(x | mu, var), evaluated elementwise over x.

    var is the variance (sigma squared), not the standard deviation. This is
    the single building block the whole algorithm leans on — every step below
    is just weighted sums of this number.
    """
    return np.exp(-0.5 * (x - mu) ** 2 / var) / np.sqrt(2.0 * np.pi * var)
# endregion


# region: e_step
def e_step(x, weights, means, variances):
    """E-step: responsibilities gamma[i, k] = P(component k produced x_i).

    For each point we weight every component's density by its mixing weight and
    normalize across components, so each row sums to 1. gamma[i, k] near 1 means
    point i almost certainly belongs to component k; 0.5 means it's on the fence.
    Returns an (N, K) matrix — the soft assignment the M-step averages over.
    """
    weighted = np.array([w * gaussian(x, m, v)
                         for w, m, v in zip(weights, means, variances)]).T  # (N, K)
    return weighted / weighted.sum(axis=1, keepdims=True)
# endregion


# region: m_step
def m_step(x, gamma):
    """M-step: re-estimate weights, means, variances as weighted averages.

    N_k is the soft count of points in component k — the column sum of the
    responsibilities. Everything is the ordinary formula for a weighted mean and
    weighted variance, with gamma[:, k] as the weights: a point that half-belongs
    to a component contributes half of itself to that component's estimates.
    """
    Nk = gamma.sum(axis=0)                                   # (K,) soft counts
    weights = Nk / len(x)
    means = (gamma * x[:, None]).sum(axis=0) / Nk
    variances = (gamma * (x[:, None] - means) ** 2).sum(axis=0) / Nk
    return weights, means, np.maximum(variances, 1e-6)       # floor guards /0
# endregion


# region: log_likelihood
def log_likelihood(x, weights, means, variances):
    """Incomplete-data log-likelihood: sum_i log sum_k pi_k N(x_i | mu_k, var_k).

    "Incomplete" because it's the likelihood of the x's alone, with the hidden
    component labels summed out. This is the single number EM optimizes, and the
    one guarantee is that it never goes down from one iteration to the next.
    """
    per = np.array([w * gaussian(x, m, v)
                    for w, m, v in zip(weights, means, variances)]).T   # (N, K)
    return float(np.log(per.sum(axis=1)).sum())
# endregion


# region: em
def em(x, init_weights, init_means, init_variances, max_iter=200, tol=1e-6):
    """EM loop: alternate E-step and M-step until the log-likelihood plateaus.

    Starts from whatever parameters you hand it — EM only finds a local optimum,
    so the init matters. Each pass computes responsibilities (E), re-estimates the
    components from them (M), and checks the log-likelihood; when it stops rising
    by more than tol, we've converged. Returns the fitted parameters, the final
    log-likelihood, and the full per-iteration history the chapter replays.
    """
    weights = np.array(init_weights, float)
    means = np.array(init_means, float)
    variances = np.array(init_variances, float)

    ll = log_likelihood(x, weights, means, variances)
    history = [_snapshot(x, weights, means, variances, ll)]
    for _ in range(max_iter):
        gamma = e_step(x, weights, means, variances)            # E-step
        weights, means, variances = m_step(x, gamma)            # M-step
        new_ll = log_likelihood(x, weights, means, variances)
        history.append(_snapshot(x, weights, means, variances, new_ll))
        if new_ll - ll < tol:                                   # plateaued
            ll = new_ll
            break
        ll = new_ll
    return weights, means, variances, ll, history
# endregion


def _snapshot(x, weights, means, variances, ll):
    """One frame of the run: the current parameters, the incomplete-data
    log-likelihood, and each point's responsibility to component 1 (the value
    the animation colors points by so you watch the soft assignment harden)."""
    resp1 = e_step(x, weights, means, variances)[:, 1]
    return {
        "weights": weights.copy(),
        "means": means.copy(),
        "variances": variances.copy(),
        "ll": ll,
        "resp1": resp1.copy(),
    }


def load_data(path="../data/mixture.csv"):
    """The 1-D mixture snapshot: one column x of numbers, plus the true source
    component of each. EM sees only x; the component column is ground truth held
    back to grade the recovered clustering at the end."""
    return pd.read_csv(path)

La versión de librería

Nadie implementa EM a mano para una mezcla gaussiana en producción, y una vez que lo construiste no te hace falta. GaussianMixture de scikit-learn es el mismo ciclo de paso E/paso M con los filos limados — una inicialización con k-means en vez de una aleatoria, varios reinicios para que una sola mala semilla no pueda decidir la respuesta, y un manejo de covarianza que generaliza más allá de nuestra varianza única. Sobre datos 1-D con dos componentes optimiza la log-verosimilitud idéntica, así que debería aterrizar donde aterrizamos nosotros:

def sklearn_gmm(x, k=2, seed=0):
    """Fit a k-component 1-D Gaussian mixture with scikit-learn's EM.

    Returns the same four things our em() returns: mixing weights, means,
    variances, and the total incomplete-data log-likelihood. score() reports the
    mean log-likelihood per sample, so we multiply by N to match our total.
    """
    X = np.asarray(x, float).reshape(-1, 1)
    gm = GaussianMixture(n_components=k, covariance_type="full",
                         max_iter=200, tol=1e-6, random_state=seed)
    gm.fit(X)
    weights = gm.weights_
    means = gm.means_.ravel()
    variances = gm.covariances_.ravel()
    ll = float(gm.score(X) * len(X))
    return weights, means, variances, ll

El único detalle es de contabilidad: score devuelve la log-verosimilitud media por punto, así que multiplicamos por N para recuperar el total que reporta nuestro em. Todo lo demás cuadra. Aquí están ambos ajustes encimados sobre los datos — las dos curvas de componente ajustadas, escaladas por sus pesos, sobre el mismo histograma del inicio:

Las dos campanas quedan justo donde tu ojo adivinó que estaban, y su mezcla traza el histograma lleno de bultos. Los parámetros recuperados están cerca de la verdad con la que construimos: medias -1.92 y 3.09 contra las verdaderas -2 y 3, pesos 0.43 y 0.57 contra los verdaderos 0.40 y 0.60. No exactos, porque 300 muestras es una extracción finita, pero eso es el ajuste, no el algoritmo — encontró las mejores gaussianas para estos datos en particular.

Desde cero contra librería

Mismos datos, mismos dos componentes, ambos corridos hasta converger — nuestro EM desde cero contra el GaussianMixture de sklearn. Cada parámetro ajustado, lado a lado:

Los pares son indistinguibles. Ambos alcanzan una log-verosimilitud de datos incompletos de -681.48; los pesos coinciden a tres decimales (0.431 y 0.569), las medias a dos (-1.92 y 3.09), las varianzas a dos (1.13 y 2.03). Eso no es coincidencia ni es un empate que yo haya redondeado hasta existir — en una mezcla limpia de dos gaussianas hay un óptimo dominante, y cualquier EM correcto que suba hasta él aterriza en los mismos parámetros, ya sea que haya arrancado de k-means como sklearn o de mi init deliberadamente malo. La única nota de contabilidad: EM es libre de llamar a la gaussiana izquierda "componente 0" o "componente 1", así que ordeno ambos ajustes por media antes de comparar. Alinéalos y son el mismo modelo.

Ahora abre el sobre sellado. Nunca ajustamos sobre las etiquetas verdaderas de componente, pero podemos calificar contra ellas: toma la asignación dura de cada punto (el componente que tenga la responsabilidad más alta) y puntúala con el adjusted Rand index, que mide el acuerdo entre dos etiquetados después de corregir por azar. Tanto nuestro clustering como el de sklearn sacan 0.87 — no un 1.0 perfecto, y no deberían, porque las dos gaussianas genuinamente se traslapan y los puntos del centro de verdad son ambiguos. EM recuperó la estructura casi tan bien como se puede recuperar; el 13% que "falla" es incertidumbre honesta que los datos mismos contienen, no una falla del ajuste.

Conclusiones

EM es la herramienta que agarras cuando lo que haría fácil tu problema es lo que no tienes. Etiquetas faltantes, mediciones faltantes, un estado oculto, una causa latente — plantea el hueco como "datos incompletos" y EM te da un ciclo que lo rellena suavemente y reajusta, una y otra vez, con la garantía de que nunca vas hacia atrás. Esa garantía es la razón entera para sabérselo de memoria: sin learning rate, sin divergencia, solo una subida monótona hasta un punto estacionario. En esta mezcla de juguete la subida fue dramática y el final fue correcto, pero guarda la corrección de Wu en el bolsillo — el final es un óptimo local, y en un problema más duro un mal arranque te deja en uno malo que se ve igual de convergido. Córrelo desde varios inits, o siémbralo con k-means, y quédate con la mejor log-verosimilitud. Ese hábito no es opcional.

La razón por la que EM se gana un capítulo propio, aparte de cualquier modelo en particular, es que está debajo de todo un estante de métodos que vas a volver a encontrar. Ajusta esta receta exacta en dos o más dimensiones con matrices de covarianza completas y obtienes modelos de mezclas gaussianas, donde cada componente es una elipse inclinada y estirada y k-means resulta ser el límite de asignación dura del mismo ciclo. Cambia las gaussianas por multinomiales sobre palabras y obtienes modelos de temas, EM descubriendo temas en documentos que nadie etiquetó. Córrelo sobre secuencias con una estructura de transición y es Baum-Welch, el algoritmo de entrenamiento de los modelos ocultos de Markov. Parecen algoritmos distintos hasta que notas que todos son el mismo paso E y paso M, escalando el mismo tipo de cota. Construye el ciclo a mano una vez, sobre dos campanas en una dimensión, y habrás construido el motor de todos ellos.