# introduccion-a-la-estadistica-bayesiana-con-pymc

Introducción a la estadística bayesiana con PyMC

Volver a la portada del curso
Lectura (Markdown)

Algoritmos Avanzados de Inferencia

Sección 4 · Inferencia Bayesiana Avanzada

Por qué Metropolis no alcanza

Metropolis-Hastings propone un salto aleatorio y lo acepta o lo rechaza. Con dos parámetros funciona; con cincuenta, casi todas las propuestas caen en zonas de baja probabilidad y se rechazan, o el paso se vuelve tan pequeño que la cadena apenas se mueve. El resultado son muestras muy autocorrelacionadas y un tamaño efectivo de muestra diminuto. Los algoritmos modernos resuelven esto usando información que Metropolis ignora: el gradiente de la posterior.

HMC: moverse siguiendo la pendiente

Hamiltonian Monte Carlo (HMC) trata los parámetros $\theta$ como la posición de una partícula y agrega un momento $p$ aleatorio. La energía total combina la energía potencial (la log-posterior negativa) y la cinética:

$$H(\theta, p) = -\log p(\theta \mid y) + \tfrac{1}{2}\, p^\top M^{-1} p$$

La partícula se desplaza sobre esa superficie siguiendo las ecuaciones de Hamilton, simuladas con el integrador leapfrog: $L$ pasos de tamaño $\varepsilon$. Como el gradiente empuja hacia las zonas de alta densidad, cada propuesta llega lejos y aun así tiene alta probabilidad de aceptarse. Por eso PyMC necesita PyTensor: sin gradientes automáticos no hay HMC.

El problema práctico es elegir $\varepsilon$ y $L$. Un paso grande acumula error numérico; uno pequeño desperdicia cómputo. Pocos pasos se parecen a una caminata aleatoria; demasiados hacen que la trayectoria regrese al punto de partida.

NUTS: HMC que se ajusta solo

El No-U-Turn Sampler (NUTS) elimina la elección de $L$: extiende la trayectoria hasta que empieza a dar la vuelta sobre sí misma. El tamaño de paso $\varepsilon$ se ajusta durante la fase de tune con dual averaging, buscando una tasa de aceptación media igual a target_accept, que en PyMC vale 0.8 por defecto. NUTS es lo que pm.sample() usa automáticamente cuando todas las variables son continuas.

Divergencias

Una transición divergente ocurre cuando el integrador falla: la energía $H$ se dispara porque la curvatura de la posterior cambia más rápido de lo que el paso $\varepsilon$ puede seguir. El caso típico es el embudo de un modelo jerárquico. Las divergencias indican que hay regiones de la posterior que el muestreador no explora bien, así que los resultados pueden estar sesgados. Subir target_accept a 0.9 o 0.95 reduce el paso; si no basta, hay que reparametrizar el modelo (por ejemplo, la forma no centrada).

Aplicación: elegir muestreador y diagnosticar

import pymc as pm
import arviz as az

with modelo:
    # NUTS de PyMC, con paso más conservador
    idata = pm.sample(target_accept=0.95, random_seed=42)
    # Alternativas: "nutpie" (Rust + Numba), "numpyro" o "blackjax" (JAX)
    idata_np = pm.sample(nuts_sampler="nutpie", random_seed=42)

n_div = int(idata.sample_stats["diverging"].sum())
print("Divergencias:", n_div)
print(az.summary(idata, kind="diagnostics"))  # r_hat, ess_bulk, ess_tail

Los muestreadores alternativos implementan el mismo algoritmo NUTS sobre otro motor. nutpie compila el modelo con Numba y suele ser el más rápido en CPU para modelos pequeños y medianos; numpyro y blackjax requieren JAX y aprovechan GPU en modelos grandes. Todos devuelven un InferenceData, así que el resto del análisis no cambia.

Inferencia variacional: velocidad a cambio de exactitud

La inferencia variacional convierte el muestreo en optimización: busca, dentro de una familia simple $q_\phi(\theta)$, la distribución más cercana a la posterior maximizando la ELBO:

$$\text{ELBO}(\phi) = \mathbb{E}_{q_\phi}[\log p(y, \theta)] - \mathbb{E}_{q_\phi}[\log q_\phi(\theta)]$$

ADVI (Automatic Differentiation Variational Inference) usa una normal con covarianza diagonal en el espacio no restringido de los parámetros. En PyMC:

with modelo:
    aprox = pm.fit(n=20000, method="advi", random_seed=42)
    idata_vi = aprox.sample(draws=1000)

print(az.summary(idata_vi, kind="stats"))

ADVI es útil con muchos datos, para prototipos rápidos o para inicializar MCMC. Su limitación es estructural: al suponer independencia entre parámetros subestima la incertidumbre y no captura correlaciones ni multimodalidad. method="fullrank_advi" recupera correlaciones con más costo. Para conclusiones finales, NUTS sigue siendo la referencia.

Diagnósticos mínimos

  • R-hat: compara la varianza dentro y entre cadenas; el criterio de Vehtari et al. (2021) exige $\hat{R} < 1.01$.
  • ESS bulk y tail: muestras efectivas para el centro y las colas de la distribución; como referencia, al menos 100 por cadena (400 con 4 cadenas).
  • Divergencias: idealmente cero después del tune; unas pocas ya justifican revisar el modelo.

Trampas comunes

  • Silenciar las divergencias: subir target_accept a 0.99 sin reparametrizar suele ocultar el problema y volver el muestreo muy lento.
  • Tratar ADVI como MCMC: sus intervalos son más estrechos de lo real; no reportes incertidumbre final con ADVI sin comparar.
  • Usar Metropolis por costumbre: con variables continuas, NUTS casi siempre gana en ESS por segundo.
  • Revisar solo R-hat: un R-hat de 1.00 con ESS tail de 50 todavía no permite confiar en los intervalos.

Cierre

El orden práctico es NUTS por defecto, nutpie cuando importa la velocidad en CPU, JAX para modelos grandes y ADVI solo cuando el costo de MCMC es prohibitivo. En cualquier caso, ningún resultado se interpreta antes de revisar R-hat, ESS y divergencias.

Recursos