Capítulo 33 de 37 · avanzado

Modelos de tópicos

De qué trata este capítulo

Aquí tienes un problema sin una sola etiqueta a la vista. Tienes una pila de documentos — correos, resúmenes de artículos, reseñas de productos, lo que sea — y sospechas que se agrupan en un puñado de temas. Unos hablan de animales, otros de computadoras, otros de comida. Nadie los etiquetó. Nadie te dijo cuántos temas hay ni qué palabras pertenecen a cuál. Solo tienes las palabras. ¿Puede un algoritmo leer la pila y devolverte los temas, junto con qué documento pertenece a cuál?

Eso es lo que hace un modelo de tópicos, y este capítulo construye el más simple que funciona. Modelamos cada documento como generado por un solo tópico oculto, y cada tópico como un sesgo sobre el vocabulario — una probabilidad que se carga fuerte hacia unas palabras y se aleja de otras. Luego lo ajustamos con EM, el mismo ciclo paso-E/paso-M del capítulo de EM, solo que ahora lo oculto no es cuál gaussiana generó un número sino cuál tópico escribió un documento. Armamos todo el ciclo a mano en NumPy sobre un corpus sintético con temas que nosotros mismos plantamos, para poder verificar al final si EM los encontró. Después lo ponemos frente a frente contra el LDA de scikit-learn, el famoso descendiente bayesiano del modelo que construimos.

Lo que hay que esperar es la animación. Arranca con cada tópico como ruido plano, cada palabra igual de probable, sin tema alguno, y puedes ver cómo los tópicos se afilan, iteración tras iteración, hasta que el tópico de animales concentra toda su masa en dog, cat y tiger, y el de comida en pizza y bread. Temas emergiendo del ruido, sin nadie supervisando. Esa es toda la idea de la página.

Un poco de historia

El linaje arranca con un artículo de 1999 de Thomas Hofmann, "Probabilistic Latent Semantic Indexing". La recuperación de texto en ese entonces se apoyaba en el análisis semántico latente, un truco puramente de álgebra lineal — tomas la matriz documento-término y la factorizas con un SVD. Funcionaba, pero no tenía probabilidades adentro ni nada que pudieras llamar una historia generativa. Hofmann le dio una. Su LSA probabilístico, pLSA, decía que un documento es una mezcla sobre tópicos latentes y que cada tópico es una distribución sobre palabras, y lo ajustó con EM. Por primera vez los "conceptos" latentes eran distribuciones de probabilidad honestas sobre las que podías razonar. Ese modelo — una mezcla de multinomiales sobre palabras, ajustada con EM — es exactamente el que construimos aquí.

pLSA tenía un hoyo, y era profundo. Aprendía una mezcla de tópicos para cada documento del conjunto de entrenamiento pero no tenía forma principiada de asignarle una a un documento que nunca había visto, porque las mezclas eran parámetros, no muestras de nada. El arreglo llegó en 2003 de la mano de David Blei, Andrew Ng y Michael Jordan en "Latent Dirichlet Allocation". Le pusieron un prior Dirichlet sobre las mezclas de tópicos, lo que convirtió a pLSA en un modelo bayesiano completamente generativo: ahora la mezcla de un documento también se muestrea, así que un documento nuevo es simplemente otra muestra de la misma historia. LDA se volvió uno de los artículos más citados del machine learning y el caballo de batalla del modelado de tópicos por una década — todo pipeline de "cuáles son los temas de este corpus" corría sobre alguna variante suya. Construimos el modelo más simple de Hofmann porque deja el mecanismo al descubierto; LDA es lo que usas en la práctica, y lo comparamos al final.

La intuición

Piensa cómo escribirías un documento si fueras un autor muy simple. Primero eliges un tema — digamos que lanzas un dado cargado y cae en "animales". Luego metes la mano en la bolsa de animales, que está llena sobre todo de las palabras dog, cat, tiger, wolf, rabbit y solo unas cuantas de todo lo demás, y sacas palabras una por una hasta que el documento tenga la longitud suficiente. El documento que sale es una bolsa de palabras cargada de palabras de animales. Quien lo lea no puede ver el tiro del dado ni la bolsa — solo ve las palabras — pero el tema se transparenta clarito en cuáles palabras aparecieron y con qué frecuencia.

Ahora corre eso al revés. Te entregan un montón de estas bolsas y no te dicen nada. No conoces los pesos del dado, no sabes qué hay en cada bolsa y no sabes de cuál bolsa salió cada documento. Pero fíjate en la misma estructura de huevo o gallina que tiene todo problema de EM. Si alguien te dijera qué tópico escribió cada documento, llenar las bolsas sería trivial — juntas todos los documentos de animales y cuentas sus palabras. Y si alguien te dijera qué hay en cada bolsa, adivinar cuál bolsa escribió un documento dado sería trivial — evalúas el documento bajo cada bolsa y eliges la que mejor le queda. No tienes ninguna de las dos, así que arrancas de la nada: empiezas con bolsas aleatorias, las usas para adivinar suavemente qué tópico escribió cada documento, rellenas las bolsas a partir de esas adivinanzas suaves, y repites. Ese ciclo es EM, y esta es la estructura que intenta encontrar: los tres tópicos plantados como un heatmap sobre el vocabulario de 15 palabras:

Tres bloques brillantes sobre la diagonal — cada tópico se concentra en sus propias cinco palabras y se queda apagado en las otras diez. Esa estructura de bloques es la señal. El trabajo de EM es recuperarla solo a partir de los conteos, sin que jamás le muestren esta imagen.

Las matemáticas

Escribe cada documento como un renglón de conteos de palabras $x_i = (x_{i1}, \dots, x_{iV})$ , donde $x_{iw}$ es cuántas veces apareció la palabra $w$ y $V$ es el tamaño del vocabulario. Modelamos el corpus con $K$ tópicos. El tópico $k$ tiene un peso a priori $\pi_k$ (el dado, con $\sum_k \pi_k = 1$ ) y una distribución de palabras $\phi_k = (\phi_{k1}, \dots, \phi_{kV})$ sobre el vocabulario (la bolsa, con $\sum_w \phi_{kw} = 1$ ). Una sola etiqueta oculta $z_i \in \{1, \dots, K\}$ dice qué tópico escribió el documento $i$ . La historia generativa — eliges un tópico, luego sacas sus palabras — da la probabilidad de un documento como una mezcla sobre esa elección oculta:

p(x_i \mid \theta) = \sum_{k=1}^{K} \pi_k \prod_{w=1}^{V} \phi_{kw}^{\,x_{iw}}

El producto interno es la verosimilitud multinomial de la bolsa bajo el tópico $k$ : cada palabra aporta su probabilidad en el tópico elevada al número de veces que ocurrió. El paso E hace la pregunta inversa — dados los tópicos actuales, ¿cuál es la probabilidad posterior de que el tópico $k$ haya escrito el documento $i$ ? Es la regla de Bayes, la verosimilitud ponderada del tópico sobre el total, y le llamamos la responsabilidad $W_{ik}$ :

W_{ik} = \frac{\pi_k \prod_{w} \phi_{kw}^{\,x_{iw}}}{\sum_{j=1}^{K} \pi_j \prod_{w} \phi_{jw}^{\,x_{iw}}}

Cada renglón $W_{i\cdot}$ suma 1 — una asignación suave del documento $i$ entre los tópicos. El paso M entonces trata esas responsabilidades como membresías fraccionarias y reajusta cada tópico juntando conteos de palabras ponderados por responsabilidad. La probabilidad actualizada de la palabra $w$ en el tópico $k$ es el conteo ponderado de esa palabra en todos los documentos, normalizado para que la distribución del tópico sume uno:

\phi_{kw} = \frac{\sum_{i} W_{ik}\, x_{iw}}{\sum_{i} W_{ik} \sum_{w'} x_{iw'}}

La actualización del prior es igual de sencilla — el nuevo peso del tópico $k$ es su responsabilidad promedio a lo largo del corpus, $\pi_k = \frac{1}{N}\sum_i W_{ik}$ . Y la cantidad que ambos pasos mejoran calladamente es la log-verosimilitud de los datos, la log-probabilidad de todos los documentos con los tópicos ocultos marginalizados:

\ell(\theta) = \sum_{i=1}^{N} \log \sum_{k=1}^{K} \pi_k \prod_{w=1}^{V} \phi_{kw}^{\,x_{iw}}

Este es el número que EM escala, y por el mismo argumento del capítulo anterior, donde el paso E aprieta una cota inferior y el paso M la sube, nunca baja de una iteración a la siguiente. Una nota práctica que le da forma al código: esos productos sobre el vocabulario son productos de muchas probabilidades chiquitas, que en punto flotante caen en underflow a cero casi de inmediato. Así que hacemos todo en espacio logarítmico, donde el producto se vuelve una suma $\sum_w x_{iw} \log \phi_{kw}$ y la normalización se vuelve un log-sum-exp. Cada función de la construcción carga logaritmos, no probabilidades crudas.

En qué es bueno y en qué no

Lo atractivo de un modelo de tópicos es que encuentra estructura en texto sin ninguna etiqueta y te la devuelve en un formato que puedes leer. La salida no es un embedding opaco ni un id de cluster. Es una lista de palabras por tópico, ordenadas por probabilidad, que una persona puede ver y nombrar. "Este tópico es dog, cat, tiger, wolf: es sobre animales". Esa interpretabilidad es rara y valiosa, y sale gratis de la estructura multinomial-sobre-palabras. Encima de eso heredas todo lo que EM te da: asignaciones suaves, que son la respuesta honesta cuando un documento queda a caballo entre dos temas, y una subida monótona sin learning rate que andar cuidando.

Los costos son los que carga todo modelo de variable latente, más filosos en texto. EM encuentra un óptimo local, y con tópicos eso es un riesgo real. Un mal arranque puede fusionar dos temas genuinos en un tópico borroso o partir un tema en dos, y la log-verosimilitud se va a estancar ahí felizmente, con cara de haber convergido. Tienes que decirle cuántos tópicos buscar, y elegir $K$ es más arte que ciencia. Y todo el modelo descansa en bag-of-words: tira el orden de las palabras por completo, así que "the dog bit the man" y "the man bit the dog" le son idénticos, y no tiene noción de que "cat" y "kitten" estén relacionadas a menos que ambas coocurran con la misma compañía. Esa última limitación es justo la brecha que los word embeddings y, más tarde, los large language models vinieron a cerrar, que es donde aterrizamos al final.

Los datos

Armé un corpus sintético para poder calificar la recuperación contra una verdad que nosotros controlamos. El vocabulario son 15 palabras en tres bloques de cinco: un bloque de animales (dog, cat, tiger, wolf, rabbit), un bloque de tech (cpu, laptop, server, code, network) y un bloque de comida (pizza, bread, cheese, coffee, sugar). Cada bloque es un tópico plantado. Cada documento se genera con la historia de la sección de intuición: sacas un tópico del prior $[0.40, 0.33, 0.27]$ — deliberadamente desigual, para que los tópicos no sean del mismo tamaño — y luego sacas una bolsa de 24 a 60 palabras de la multinomial de ese tópico, que pone peso fuerte y variado sobre sus propias cinco palabras y solo una fuga pequeña sobre las otras diez. El resultado son 150 documentos, 6,038 tokens de palabra en total, repartidos 61 / 44 / 45 entre los tres tópicos. El true_topic de cada documento se registra y luego se sella: EM nunca lo ve, y abrimos el sobre solo al final para calificar qué tan bien coinciden los tópicos recuperados con los plantados.

Aquí está la matriz cruda de conteos documento-término — un renglón por documento, una columna por palabra, ordenada para que los documentos se agrupen por su tópico verdadero. El color es el conteo de la palabra, y el brillo en la diagonal por bloques es la misma estructura plantada de antes, solo que ahora son los datos reales que EM va a leer, con ruido y todo:

Léela de arriba abajo y puedes ver las tres bandas: el grupo de documentos de arriba enciende las cinco columnas de la izquierda, el grupo de en medio las siguientes cinco, el de abajo las últimas cinco. Ese es el ground truth. Pero recuerda que EM ve estos renglones sin ningún orden particular y sin etiquetas de banda — para el algoritmo esto son 150 renglones de 15 números, y tiene que descubrir que caen en tres grupos y de qué trata cada grupo.

Constrúyelo, una función a la vez

Seis funciones, todas en espacio logarítmico. De dónde partimos, los dos pasos de EM, el score que escalan y el ciclo. Todo empieza con la inicialización, y para un modelo de tópicos el init es la razón entera por la que vale la pena ver la animación:

def init_params(X, k, rng):
    """The starting point: nearly flat topics and a uniform prior.

    Each topic-word distribution is uniform over the vocabulary jittered by a
    little noise, so no two topics are identical (EM can't split them if they
    are) but every one starts almost featureless. This is what makes the
    animation worth watching — the topics emerge from that flat noise.
    """
    d = X.shape[1]
    P = rng.uniform(0.95, 1.05, size=(k, d))
    P /= P.sum(axis=1, keepdims=True)
    log_pi = np.log(np.full(k, 1.0 / k))
    return np.log(P), log_pi

Cada tópico arranca casi exactamente como la distribución uniforme — cada palabra alrededor de $1/15$ — sacudido por un pelín de ruido. El jitter importa: si dos tópicos empezaran idénticos byte por byte, EM jamás podría separarlos, porque las responsabilidades también serían idénticas y la simetría nunca se rompería. Así que los arrancamos casi planos pero no del todo, y dejamos que EM encuentre la asimetría. El prior arranca uniforme. Desde ese inicio sin rasgos, los tópicos tienen que hacer crecer sus temas.

El paso E calcula las responsabilidades. Para cada documento y cada tópico, el prior más la log-verosimilitud de la bolsa bajo ese tópico, el log-conjunto, y normalizar cada renglón sobre los tópicos con un log-sum-exp lo convierte en el log posterior:

def e_step(X, log_P, log_pi):
    """E-step: the responsibility of each topic for each document.

    X is (N, d) integer counts, log_P is (k, d) log topic-word probabilities,
    log_pi is (k,) log topic priors. For document i and topic j the unnormalized
    log-joint is log(pi_j) + sum_w x_iw * log(P_jw) — the prior plus the
    log-likelihood of the bag under that topic's multinomial. Normalizing each
    row over topics (subtract the log-sum-exp) turns it into log W_ij, the log
    posterior probability that topic j generated document i. Returns (N, k).
    """
    log_joint = X @ log_P.T + log_pi[None, :]              # (N, k)
    log_evidence = logsumexp(log_joint, axis=1, keepdims=True)  # (N, 1)
    return log_joint - log_evidence

Todo es una sola multiplicación de matrices: X @ log_P.T calcula, para cada par documento-tópico de un jalón, la suma $\sum_w x_{iw} \log \phi_{kw}$ — la log-verosimilitud multinomial de cada bolsa bajo cada tópico. Súmale el log del prior, réstale el log-sum-exp por renglón, y tienes log $W_{ik}$ . El paso M se parte entonces en dos piezas. Primero los tópicos: ponderas los conteos de cada documento por la responsabilidad que cada tópico toma sobre él, los juntas y normalizas cada tópico para que sume uno:

def m_step_topics(X, log_W, eps=1e-12):
    """M-step, part one: re-estimate each topic's word distribution.

    Weight every document's word counts by the responsibility the topic takes for
    it and pool them: E_jw = sum_i W_ij * x_iw (a tiny eps keeps the log finite
    for words no document credited to the topic). Normalizing each topic's row to
    sum to one gives the updated multinomial. Returns log_P, shape (k, d).
    """
    W = np.exp(log_W)                                     # (N, k) responsibilities
    weighted = W.T @ X + eps                              # (k, d) pooled counts
    P = weighted / weighted.sum(axis=1, keepdims=True)
    return np.log(P)

W.T @ X es el conteo agrupado y ponderado por responsabilidad de cada palabra en cada tópico, exactamente el numerador de la actualización de $\phi_{kw}$ , y dividir cada renglón entre su suma hace la normalización. El eps minúsculo mantiene el log finito para una palabra que ningún documento le acreditó a un tópico. Segundo los priors, la responsabilidad promedio que cargó cada tópico, hecho en espacio logarítmico:

def m_step_priors(log_W):
    """M-step, part two: re-estimate the topic priors.

    The updated prior for topic j is the average responsibility it carried across
    all documents: pi_j = (1/N) * sum_i W_ij. Done in log-space with logsumexp
    down the document axis. Returns log_pi, shape (k,).
    """
    N = log_W.shape[0]
    return logsumexp(log_W, axis=0) - np.log(N)           # (k,)

Ahora el score que ambos pasos escalan, la log-verosimilitud de los datos — marginaliza el tópico de cada documento con un log-sum-exp, luego suma sobre los documentos:

def log_likelihood(X, log_P, log_pi):
    """Data log-likelihood under the current parameters.

    Marginalize the latent topic out of each document — sum_j pi_j * P(doc|topic
    j) — in log-space with logsumexp, then sum over documents. This is the number
    EM is climbing; it must never go down from one full iteration to the next.
    """
    log_joint = X @ log_P.T + log_pi[None, :]              # (N, k)
    return float(logsumexp(log_joint, axis=1).sum())

Este es el número que nunca debe caer; si alguna vez cae, hay un bug en alguno de los pasos. Por último el ciclo que los amarra:

def em(X, k, rng, max_iter=100, tol=1e-4):
    """The EM loop: E-step, M-step, repeat until the log-likelihood stops rising.

    Records a snapshot every iteration — the full topic-word table, the priors,
    and the data log-likelihood — which the chapter replays frame by frame as the
    topics sharpen out of the flat start. Returns the hard document assignments
    (argmax responsibility), the recovered topics P (k, d), the priors pi (k,),
    the final log-likelihood, and the snapshot history.
    """
    log_P, log_pi = init_params(X, k, rng)
    history = [_snapshot(X, log_P, log_pi)]
    prev_ll = history[0]["ll"]
    for _ in range(max_iter):
        log_W = e_step(X, log_P, log_pi)          # responsibilities
        log_P = m_step_topics(X, log_W)           # re-estimate topics
        log_pi = m_step_priors(log_W)             # re-estimate priors
        snap = _snapshot(X, log_P, log_pi)
        history.append(snap)
        if snap["ll"] - prev_ll < tol:            # converged — no more to gain
            break
        prev_ll = snap["ll"]
    log_W = e_step(X, log_P, log_pi)
    labels = np.argmax(log_W, axis=1)
    return labels, np.exp(log_P), np.exp(log_pi), history[-1]["ll"], history

Paso E, paso M de tópicos, paso M de priors, score, repetir hasta que la log-verosimilitud deje de subir más que un suspiro. Toma una foto de la tabla completa de tópicos, de los priors y de la log-verosimilitud en cada iteración — esa historia es exactamente lo que reproduce la animación. Al final regresa las asignaciones duras de documentos (el tópico más responsable de cada documento), los tópicos recuperados, los priors, la log-verosimilitud final y la historia.

Míralo funcionar

Esta es la pieza central. Abajo hay una corrida real de la función em de arriba, un cuadro por iteración, arrancando desde el init aleatorio casi plano. Tres paneles, uno por tópico recuperado, acomodados de modo que después de converger el panel uno se lea como el tópico de animales, el dos como tech y el tres como comida. Dentro de cada panel cada barra es una palabra, y la barra va coloreada por el tópico que realmente plantó esa palabra — cian para palabras de animales, morado para palabras de tech, naranja para palabras de comida. Así que un panel que se enganchó al tema de animales debería terminar con puras barras cian altas y todo lo demás chiquito.

Mira el primer cuadro: cada panel es una cerca plana de barras cortas de colores mezclados, porque en el arranque aleatorio ningún tópico favorece nada, cada palabra sentada cerca de la tasa de azar de $1/15$ . Luego dale play y mira cómo se afilan las barras. Los colores se ordenan solos: las barras cian de un panel crecen mientras sus barras moradas y naranjas se encogen hasta desaparecer, y lo mismo para los otros dos paneles con sus colores. Para el último cuadro cada panel es un pico limpio de un solo color — el tópico concentró toda su masa de probabilidad en sus propias cinco palabras. El pie de la animación va marcando la iteración y la log-verosimilitud de los datos subiendo debajo de todo.

La emergencia es rápida porque los tópicos plantados están bien separados, y eso vale la pena verlo claro. En el arranque aleatorio la log-verosimilitud es −16,365 y cada tópico es ruido. Una pasada E/M la lleva a −15,015 y ya se alcanza a ver que las barras empiezan a inclinarse hacia sus colores. La segunda pasada salta a −12,688, el cuadro donde los temas se vuelven inconfundibles, cada panel ya visiblemente dominado por un color. La tercera llega a −11,648, y la cuarta lo confirma: ya no hay cambio, los tópicos quedaron fijos. Cuatro iteraciones del ruido a la respuesta.

Ahora busca un cuadro donde la log-verosimilitud baje. No hay ninguno. La misma garantía de todo EM, hecha visible otra vez. Aquí está esa subida por su cuenta:

La forma familiar de EM: una subida empinada al principio, luego un codo, luego plano. Aquí se aplana después de cuatro iteraciones en lugar de las veinte del capítulo de gaussianas, porque los tópicos están limpiamente separados y las responsabilidades se pegan a la casi-certeza casi de inmediato. En un corpus más sucio, con temas que comparten vocabulario, verías la cola larga y lenta en su lugar.

La implementación completa

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

"""Mixture-of-multinomials topic model, fit by EM — built from scratch.

Each document is a bag of words generated by ONE latent topic. A topic is a
multinomial distribution over the vocabulary. Given the count matrix X we don't
know which topic wrote which document, nor what the topics are — so we alternate:

  E-step: given the current topics and priors, compute the posterior probability
          that each document came from each topic (its responsibilities).
  M-step: given those responsibilities, re-estimate each topic's word
          distribution and the topic priors from responsibility-weighted counts.

Everything runs in log-space with scipy's logsumexp for numerical stability —
the per-word products become sums of log-probabilities and never underflow.
Pure NumPy for the model; the `# region:` markers are what the chapter includes.
"""

import numpy as np
from scipy.special import logsumexp


# region: e_step
def e_step(X, log_P, log_pi):
    """E-step: the responsibility of each topic for each document.

    X is (N, d) integer counts, log_P is (k, d) log topic-word probabilities,
    log_pi is (k,) log topic priors. For document i and topic j the unnormalized
    log-joint is log(pi_j) + sum_w x_iw * log(P_jw) — the prior plus the
    log-likelihood of the bag under that topic's multinomial. Normalizing each
    row over topics (subtract the log-sum-exp) turns it into log W_ij, the log
    posterior probability that topic j generated document i. Returns (N, k).
    """
    log_joint = X @ log_P.T + log_pi[None, :]              # (N, k)
    log_evidence = logsumexp(log_joint, axis=1, keepdims=True)  # (N, 1)
    return log_joint - log_evidence
# endregion


# region: log_likelihood
def log_likelihood(X, log_P, log_pi):
    """Data log-likelihood under the current parameters.

    Marginalize the latent topic out of each document — sum_j pi_j * P(doc|topic
    j) — in log-space with logsumexp, then sum over documents. This is the number
    EM is climbing; it must never go down from one full iteration to the next.
    """
    log_joint = X @ log_P.T + log_pi[None, :]              # (N, k)
    return float(logsumexp(log_joint, axis=1).sum())
# endregion


# region: m_step_topics
def m_step_topics(X, log_W, eps=1e-12):
    """M-step, part one: re-estimate each topic's word distribution.

    Weight every document's word counts by the responsibility the topic takes for
    it and pool them: E_jw = sum_i W_ij * x_iw (a tiny eps keeps the log finite
    for words no document credited to the topic). Normalizing each topic's row to
    sum to one gives the updated multinomial. Returns log_P, shape (k, d).
    """
    W = np.exp(log_W)                                     # (N, k) responsibilities
    weighted = W.T @ X + eps                              # (k, d) pooled counts
    P = weighted / weighted.sum(axis=1, keepdims=True)
    return np.log(P)
# endregion


# region: m_step_priors
def m_step_priors(log_W):
    """M-step, part two: re-estimate the topic priors.

    The updated prior for topic j is the average responsibility it carried across
    all documents: pi_j = (1/N) * sum_i W_ij. Done in log-space with logsumexp
    down the document axis. Returns log_pi, shape (k,).
    """
    N = log_W.shape[0]
    return logsumexp(log_W, axis=0) - np.log(N)           # (k,)
# endregion


# region: init_params
def init_params(X, k, rng):
    """The starting point: nearly flat topics and a uniform prior.

    Each topic-word distribution is uniform over the vocabulary jittered by a
    little noise, so no two topics are identical (EM can't split them if they
    are) but every one starts almost featureless. This is what makes the
    animation worth watching — the topics emerge from that flat noise.
    """
    d = X.shape[1]
    P = rng.uniform(0.95, 1.05, size=(k, d))
    P /= P.sum(axis=1, keepdims=True)
    log_pi = np.log(np.full(k, 1.0 / k))
    return np.log(P), log_pi
# endregion


# region: em
def em(X, k, rng, max_iter=100, tol=1e-4):
    """The EM loop: E-step, M-step, repeat until the log-likelihood stops rising.

    Records a snapshot every iteration — the full topic-word table, the priors,
    and the data log-likelihood — which the chapter replays frame by frame as the
    topics sharpen out of the flat start. Returns the hard document assignments
    (argmax responsibility), the recovered topics P (k, d), the priors pi (k,),
    the final log-likelihood, and the snapshot history.
    """
    log_P, log_pi = init_params(X, k, rng)
    history = [_snapshot(X, log_P, log_pi)]
    prev_ll = history[0]["ll"]
    for _ in range(max_iter):
        log_W = e_step(X, log_P, log_pi)          # responsibilities
        log_P = m_step_topics(X, log_W)           # re-estimate topics
        log_pi = m_step_priors(log_W)             # re-estimate priors
        snap = _snapshot(X, log_P, log_pi)
        history.append(snap)
        if snap["ll"] - prev_ll < tol:            # converged — no more to gain
            break
        prev_ll = snap["ll"]
    log_W = e_step(X, log_P, log_pi)
    labels = np.argmax(log_W, axis=1)
    return labels, np.exp(log_P), np.exp(log_pi), history[-1]["ll"], history
# endregion


def _snapshot(X, log_P, log_pi):
    """One frame of the run: the topics, the priors, and the log-likelihood."""
    return {
        "P": np.exp(log_P),
        "pi": np.exp(log_pi),
        "ll": log_likelihood(X, log_P, log_pi),
    }


def load_counts(path="../data/counts.csv"):
    """The committed document-term matrix plus the true topic per document.

    Returns (X, true_topic, vocab). The true_topic column is ground truth we
    NEVER fit on — EM sees only the counts. It exists to grade recovery.
    """
    import pandas as pd
    df = pd.read_csv(path)
    vocab = [c for c in df.columns if c != "true_topic"]
    X = df[vocab].to_numpy(float)
    return X, df["true_topic"].to_numpy(int), vocab

La versión con librería

Nadie escribe a mano un modelo de tópicos en producción, y la herramienta estándar no es exactamente el modelo que construimos. El LatentDirichletAllocation de scikit-learn es LDA, el descendiente bayesiano de la sección de historia. La diferencia es real y vale la pena decirla: nuestro modelo asigna cada documento a un solo tópico latente, mientras que LDA deja que un documento sea una mezcla de tópicos, cada tópico sacado de un prior Dirichlet, y lo ajusta por inferencia variacional en lugar de EM a secas. Es el modelo más rico. En un corpus como el nuestro, donde cada documento sí fue escrito por un solo tópico, ambos deberían recuperar la misma estructura plantada — pero no esperes números idénticos, porque están optimizando objetivos distintos bajo supuestos distintos.

def sklearn_lda(counts, k, seed=0):
    """Fit LDA with k topics on the integer count matrix.

    Returns the hard document assignments (the topic each document loads on most)
    and the normalized topic-word distributions (k, d) — the same two things we
    grade our own topics on.
    """
    lda = LatentDirichletAllocation(
        n_components=k, learning_method="batch", max_iter=50, random_state=seed
    )
    doc_topic = lda.fit_transform(counts)                 # (N, k) topic mixtures
    topic_word = lda.components_ / lda.components_.sum(axis=1, keepdims=True)
    labels = doc_topic.argmax(axis=1)
    return labels, topic_word

Le pedimos las mismas dos cosas con las que calificamos nuestro propio modelo: un tópico duro por documento (aquel en el que el documento carga más) y las distribuciones tópico-palabra normalizadas. Todo lo demás, los priors Dirichlet y la cota variacional, es asunto de LDA, y lo tratamos como una referencia de caja negra.

Desde cero contra la librería

Ambos modelos ven los mismos 150 documentos y se les piden tres tópicos. Para calificarlos abrimos el sobre sellado de etiquetas verdaderas y usamos dos medidas. La primera es el índice de Rand ajustado entre la asignación dura de documentos de cada modelo y los tópicos verdaderos — acuerdo entre dos etiquetados, corregido por azar, donde 1.0 es perfecto. La segunda es qué tan cerca está cada distribución tópico-palabra recuperada de la que la plantó, medida por similitud coseno después de emparejar los tópicos recuperados con los plantados (ambos modelos devuelven los tópicos en orden arbitrario, así que primero los apareamos por máxima similitud).

Ambos modelos sacan un índice de Rand ajustado perfecto de 1.00: cada documento acomodado en el tópico correcto, los 150. Y ambos recuperan las distribuciones de palabras casi exactamente: coseno promedio de 0.998 para nuestro EM desde cero, 0.998 para LDA. Nuestro modelo hasta recupera el prior desigual del que nunca le dijeron nada, aterrizando en pesos de tópico de 0.41, 0.29, 0.30 contra los plantados 0.40, 0.33, 0.27. En este corpus limpio los dos modelos son indistinguibles al nivel de quién tiene la razón, que es justamente el punto de un corpus limpio. Confirma que ambas implementaciones son correctas antes de que confíes en cualquiera de las dos sobre datos donde la respuesta se desconoce. Para ver la recuperación palabra por palabra, aquí están los tópicos de nuestro modelo encimados sobre la verdad plantada, un panel por tópico:

Las barras recuperadas caen justo encima de las plantadas en cada panel. EM, sin nada más que 150 bolsas de palabras, reconstruyó los tres tópicos que escondimos adentro: las palabras correctas, los pesos correctos, los documentos correctos en cada uno, sin supervisión en ningún paso. Ese es el truco completo funcionando.

Conclusiones

Un modelo de tópicos es EM disfrazado de texto. Cambia las gaussianas del capítulo anterior por multinomiales sobre palabras y exactamente el mismo ciclo paso-E/paso-M descubre temas en un corpus sin etiquetar: el paso E adivina suavemente qué tópico escribió cada documento, el paso M rellena los tópicos con esas adivinanzas, y la log-verosimilitud de los datos sube monótonamente hasta un punto estacionario. Constrúyelo una vez y ya tienes el mecanismo detrás de pLSA y, con un prior bayesiano atornillado encima, de LDA, el modelo que llevó la exploración de corpus por una década. La recuperación en nuestro corpus plantado fue perfecta porque lo hicimos limpio; en texto real peleas con los demonios de siempre de EM, los óptimos locales que fusionan o parten temas genuinos, y la eterna pregunta de cuántos tópicos pedir.

El lugar honesto para terminar es dónde queda esto hoy. Las dos familias no siguieron parejas. pLSA es la estimación puntual de máxima verosimilitud, rápida y simple, exactamente lo que construimos, pero en el fondo es un modelo solo del conjunto de entrenamiento y hace overfitting sin un prior. LDA respondió a eso con tratamiento bayesiano completo, y por años fue la opción por defecto. Luego se movió el piso. Ambos modelos comparten el supuesto lisiado que ya señalamos: bag-of-words, que no sabe nada del orden de las palabras ni de que dos palabras puedan significar lo mismo. Los word embeddings empezaron a cerrar esa brecha, y los large language models la cerraron por completo. Pídele a un modelo moderno que agrupe o resuma un corpus por tema y lee los documentos, con significado y sintaxis incluidos, de una forma que ninguna multinomial podría. Los modelos de tópicos clásicos no desaparecieron; siguen ganándose el pan donde necesitas velocidad, listas de palabras interpretables o un modelo lo bastante chico como para razonar sobre él, y los pipelines de clustering de texto corto todavía los corren. Pero la frontera de "de qué trata esta pila de texto" se movió hacia representaciones que entienden las palabras en lugar de solo contarlas. Conocer la versión con EM es cómo entiendes qué fue lo que esas representaciones reemplazaron, y por qué el reemplazo fue un salto tan grande.