¿Por qué se necesita Forward-Backward por separado?

En M18, utilizamos Viterbi para extraer una única secuencia de estados más probable dada la secuencia observada. La secuencia óptima para GCAT fue H H L L.

Sin embargo, en la práctica también surgen frecuentemente otras preguntas:

• ¿Cuál es exactamente la probabilidad de que el estado en el tiempo 3 sea H?
• ¿Cuál es la probabilidad total de que toda la secuencia observada se genere con este HMM?
• ¿Es la secuencia óptima de Viterbi abrumadoramente superior, o es una victoria ajustada?

Para responder a estas preguntas, debemos conocer la probabilidad posterior de cada estado en cada instante. El algoritmo Forward-Backward calcula exactamente esto. Además, el algoritmo Baum-Welch EM, que aprenderemos en la siguiente sección (M20), utiliza los resultados de Forward-Backward como insumo. Esta sección constituye la columna vertebral de la serie sobre HMM.

Tres nuevas definiciones

Utilizamos la misma tabla de probabilidades de M18: estado {H, L}, observación {A, C, G, T}, y las tablas de probabilidad inicial, de transición y de emisión.

Probabilidad Forward α_t(s):

\alpha_t(s) = P(o_1, o_2, ..., o_t, \text{estado}_t = s)

Es la probabilidad conjunta de generar las observaciones desde el tiempo 1 hasta t, y que el estado en el tiempo t sea s. Se rellena mediante programación dinámica (DP) hacia adelante.

Probabilidad Backward β_t(s):

\beta_t(s) = P(o_{t+1}, o_{t+2}, ..., o_T | \text{estado}_t = s)

Es la probabilidad condicional de que se generen las observaciones posteriores o_{t+1}~o_T, dado que el estado en el tiempo t es s. Se rellena mediante DP hacia atrás.

Probabilidad Posterior γ_t(s):

\gamma_t(s) = P(\text{estado}_t = s | o_1, ..., o_T) = \frac{\alpha_t(s) \cdot \beta_t(s)}{P(O)}

Es la probabilidad de que el estado en el tiempo t sea s, condicionada a haber observado ya toda la secuencia completa. Este es el valor más frecuentemente necesario en la práctica.

Forward — DP hacia adelante

Ecuación de recurrencia (donde el max de Viterbi se convierte en sum):

\alpha_t(s) = \left[ \sum_{s'} \alpha_{t-1}(s') \cdot a(s' \to s) \right] \cdot b(s, o_t)

Condiciones iniciales:

\alpha_1(s) = \pi(s) \cdot b(s, o_1)

Solo difiere de Viterbi en un punto: almacena la suma de probabilidades de todas las rutas posibles que pueden llegar, en lugar de una única ruta óptima. El resto del flujo es completamente idéntico.

Calculemos manualmente con la observación GCAT (usando la misma tabla de probabilidades que en M18).

t=1, G: \alpha_1(H) = 0.5 \times 0.4 = 0.20 \alpha_1(L) = 0.5 \times 0.1 = 0.05

t=2, C: \alpha_2(H) = (0.20 \times 0.7 + 0.05 \times 0.4) \times 0.4 = 0.16 \times 0.4 = 0.064 \alpha_2(L) = (0.20 \times 0.3 + 0.05 \times 0.6) \times 0.1 = 0.09 \times 0.1 = 0.009

Se observa la diferencia de que Viterbi eligió 0.14 con el máximo, mientras que Forward utilizó la suma de las dos trayectorias, 0.16.

t=3, A: \alpha_3(H) = (0.064 \times 0.7 + 0.009 \times 0.4) \times 0.1 = 0.0484 \times 0.1 = 0.00484 \alpha_3(L) = (0.064 \times 0.3 + 0.009 \times 0.6) \times 0.4 = 0.0246 \times 0.4 = 0.00984

t=4, T: \alpha_4(H) = (0.00484 \times 0.7 + 0.00984 \times 0.4) \times 0.1 = 0.007324 \times 0.1 = 0.0007324 \alpha_4(L) = (0.00484 \times 0.3 + 0.00984 \times 0.6) \times 0.4 = 0.007356 \times 0.4 = 0.002942

Resumen de la tabla:

Probabilidad de la observación completa: P(O) = \sum_s \alpha_T(s) = 0.0007324 + 0.002942 = 0.003676

Este valor representa "la probabilidad total de que este HMM genere GCAT". Se utiliza directamente en la comparación de modelos y las pruebas de razón de verosimilitud.

Backward — Programación dinámica hacia atrás

Ecuación de recurrencia:

\beta_t(s) = \sum_{s'} a(s \to s') \cdot b(s', o_{t+1}) \cdot \beta_{t+1}(s')

Condición de frontera (se llena desde el final):

\beta_T(s) = 1

Después del tiempo T no hay observaciones, por lo que la probabilidad es 1. Se llena retrocediendo en el tiempo.

t=T=4: \beta_4(H) = \beta_4(L) = 1

t=3, siguiente observación T: \beta_3(H) = a(H\to H) b(H, T) \beta_4(H) + a(H\to L) b(L, T) \beta_4(L) = 0.7 \times 0.1 \times 1 + 0.3 \times 0.4 \times 1 = 0.07 + 0.12 = 0.19 \beta_3(L) = 0.4 \times 0.1 \times 1 + 0.6 \times 0.4 \times 1 = 0.04 + 0.24 = 0.28

t=2, siguiente observación A: \beta_2(H) = 0.7 \times 0.1 \times 0.19 + 0.3 \times 0.4 \times 0.28 = 0.0133 + 0.0336 = 0.0469 \beta_2(L) = 0.4 \times 0.1 \times 0.19 + 0.6 \times 0.4 \times 0.28 = 0.0076 + 0.0672 = 0.0748

t=1, siguiente observación C: \beta_1(H) = 0.7 \times 0.4 \times 0.0469 + 0.3 \times 0.1 \times 0.0748 = 0.013132 + 0.002244 = 0.015376 \beta_1(L) = 0.4 \times 0.4 \times 0.0469 + 0.6 \times 0.1 \times 0.0748 = 0.007504 + 0.004488 = 0.011992

Verificación: P(O) = Σ_s π(s) · b(s, o_1) · β_1(s) ≠ α_1(H)·β_1(H)/… sino, más intuitivamente:

P(O) = \sum_s \pi(s) \cdot b(s, o_1) \cdot \beta_1(s) = 0.5 \times 0.4 \times 0.015376 + 0.5 \times 0.1 \times 0.011992 = 0.003075 + 0.000600 = 0.003675

Coincide con 0.003676 obtenido por el método Forward dentro del error de redondeo. Esta es la primera comprobación de sanidad del método Forward-Backward: los valores de P(O) calculados por ambos métodos deben ser iguales.

Posterior — Producto de ambas direcciones

\gamma_t(s) = \frac{\alpha_t(s) \cdot \beta_t(s)}{P(O)}

Extraigamos el posterior para cada estado en el instante 3 (P(O) = 0.003676).

\gamma_3(H) = \frac{0.00484 \times 0.19}{0.003676} = \frac{0.0009196}{0.003676} = 0.250 \gamma_3(L) = \frac{0.00984 \times 0.28}{0.003676} = \frac{0.002755}{0.003676} = 0.749

La suma es 1.000 (redondeado). En el instante 3, hay un 75% de probabilidad de estar en el estado L y un 25% de probabilidad de estar en el estado H. Viterbi determinó que este instante corresponde al estado L, mientras que Forward-Backward proporciona una distribución de probabilidades como respuesta.

Si extraemos el posterior de cada instante y estado, obtenemos la decodificación posterior, que consiste en un arreglo formado por los estados con mayor probabilidad posterior en cada instante. Puede diferir del resultado de Viterbi; en la práctica, se suelen consultar ambos para tomar decisiones.

Escala logarítmica y el truco log-sum-exp

Viterbi utiliza directamente el máximo en el dominio logarítmico. Forward debe manejar las sumas en el dominio logarítmico, y la herramienta necesaria para ello es log-sum-exp.

\log(e^a + e^b) = \max(a, b) + \log(1 + e^{-|a-b|})

En Python, se puede usar scipy.special.logsumexp o implementarlo directamente.

import math

def logsumexp(log_values):
    m = max(log_values)
    if m == float('-inf'):
        return float('-inf')
    return m + math.log(sum(math.exp(v - m) for v in log_values))

Si multiplicamos las probabilidades directamente sin este truco, el resultado se desbordará hacia abajo (underflow) y será 0 incluso con longitudes de secuencia ligeramente mayores.

Implementación en Python (~40 líneas)

import math

def logsumexp(vals):
    m = max(vals)
    return m + math.log(sum(math.exp(v - m) for v in vals)) if m != float('-inf') else float('-inf')

def forward_backward(obs, N, start_p, trans_p, emit_p):
    T = len(obs)
    log = math.log

    # Forward
    alpha = [[0.0] * N for _ in range(T)]
    for s in range(N):
        alpha[0][s] = log(start_p[s]) + log(emit_p[s][obs[0]])
    for t in range(1, T):
        for s in range(N):
            terms = [alpha[t-1][sp] + log(trans_p[sp][s]) for sp in range(N)]
            alpha[t][s] = logsumexp(terms) + log(emit_p[s][obs[t]])

    log_P_O = logsumexp(alpha[T-1])

    # Backward
    beta = [[0.0] * N for _ in range(T)]
    for s in range(N):
        beta[T-1][s] = 0.0   # log(1) = 0
    for t in range(T-2, -1, -1):
        for s in range(N):
            terms = [log(trans_p[s][sp]) + log(emit_p[sp][obs[t+1]]) + beta[t+1][sp]
                     for sp in range(N)]
            beta[t][s] = logsumexp(terms)

    # Posterior
    gamma = [[math.exp(alpha[t][s] + beta[t][s] - log_P_O) for s in range(N)] for t in range(T)]

    return alpha, beta, gamma, log_P_O

# Ejemplo anterior
obs_map = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
obs = [obs_map[c] for c in 'GCAT']
start_p = [0.5, 0.5]
trans_p = [[0.7, 0.3], [0.4, 0.6]]
emit_p  = [[0.1, 0.4, 0.4, 0.1], [0.4, 0.1, 0.1, 0.4]]

alpha, beta, gamma, log_P_O = forward_backward(obs, 2, start_p, trans_p, emit_p)
print('P(O) =', math.exp(log_P_O))
for t, g in enumerate(gamma):
    print(f't={t+1}: P(H)={g[0]:.3f}  P(L)={g[1]:.3f}')

Al ejecutarlo se obtiene la probabilidad posterior de cada instante. Los resultados coincidirán con los valores calculados manualmente: γ_3(H)=0.250 y γ_3(L)=0.749.

Complejidad

Forward y Backward requieren cada uno O(T × N²) de tiempo y O(T × N) de espacio, la misma clase de complejidad que Viterbi. Ejecutar ambos solo añade un factor constante. La elegancia de los algoritmos HMM surge de esta estructura compartida.

Correspondencia en informática: combinación de forward y backward

Esta es la versión probabilística de la combinación DP hacia delante + DP hacia atrás tratada en DryBench.

• Forward = probabilidad conjunta desde el inicio hasta t (probabilidad del prefijo)
• Backward = probabilidad condicional posterior a t (probabilidad del sufijo)
• Producto = marginal en el instante t condicionado a todas las observaciones

Esta estructura aparece en muchos lugares: la retropropagación de redes neuronales, el cálculo de marginales en CRF e incluso las iteraciones forward/backward de PageRank para rastreadores web. Calcular las probabilidades de llegada desde ambos sentidos y multiplicarlas para obtener la marginal es un patrón que se redescubre en todos los modelos gráficos probabilísticos.

Obtención de resultados en Rosalind

Abramos el problema HMMFB de Rosalind. Pide calcular la probabilidad Forward y obtener la probabilidad total de la observación P(O). Solo hay que mostrar log_P_O del código anterior.

Posibles continuaciones

• Próximo artículo (M20): Baum-Welch EM — aprende los parámetros del HMM (π, a, b) a partir de α, β y γ de Forward-Backward. Es un bucle de autoaprendizaje que ajusta automáticamente las tablas de probabilidad utilizando solo las observaciones.
• Dos artículos más adelante (M21): Profile HMM + detección de genes — aplicación práctica para detectar estructuras génicas y buscar dominios proteicos (HMMER).
• Extensión relacionada: CRF (Conditional Random Field) — extensión condicional de los HMM y estándar para el etiquetado de secuencias en NLP y bioinformática.

Para profundizar

• Rabiner (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE — las secciones III-A y III-B presentan la derivación canónica de Forward-Backward.
• Harvard STAT115 — serie sobre HMM de la profesora Xiaole Shirley Liu (con subtítulos completos). Desarrolla Forward-Backward paso a paso en la pizarra.
• Durbin et al. Biological Sequence Analysis Capítulo 3 — descripción estándar del contexto bio.
• Michael Collins Lecture Notes (Columbia) — nueva deducción de Forward-Backward desde la perspectiva de CRF. Útil como puente de HMM a CRF.

Antes de pasar al siguiente capítulo, repasemos los cálculos manuales de este capítulo. Debe poder explicar de dónde proviene la diferencia entre α_2(H)=0.064 en Forward y V_2(H)=0.056 en Viterbi; esto le dará una intuición sobre por qué EM funciona bien en M20.