Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Conceptos clave de inferencia bayesiana

Un modelo frecuentista responde «qué tan raros serían estos datos si el parámetro valiera X». Un modelo bayesiano responde la pregunta que suelen hacerte en el trabajo: «dado lo que observé, ¿qué valores del parámetro son plausibles y con qué probabilidad?». Para eso necesitas prior, verosimilitud, posterior, predictiva posterior y un método de cálculo.

Bayes con densidades

En las videolecciones viste el Teorema de Bayes con eventos. Con parámetros continuos $\theta$ y datos $y$ se escribe con densidades:

$$p(\theta \mid y) = \frac{p(y \mid \theta)\, p(\theta)}{p(y)} \propto p(y \mid \theta)\, p(\theta)$$

• Prior $p(\theta)$: lo que se sabe del parámetro antes de ver los datos; débilmente informativo (descarta valores absurdos) o informativo (resume estudios previos).
• Verosimilitud $p(y \mid \theta)$: cómo se generan los datos. Como función de $\theta$ no integra 1.
• Posterior $p(\theta \mid y)$: la distribución actualizada; el resultado completo, no un solo número.
• Evidencia $p(y) = \int p(y \mid \theta)\, p(\theta)\, d\theta$: la constante que hace que el posterior integre 1. Es la parte difícil de calcular.

Aplicación: aproximación por rejilla

Toma el estudio de sueño del curso: 18 sujetos y su tiempo de reacción en milisegundos. En el día 0 quieres estimar la media poblacional $\mu$. Supón un prior $\mu \sim \mathcal{N}(300, 50)$ y una verosimilitud normal con $\sigma = 35$ conocida. Con un solo parámetro basta evaluar prior por verosimilitud en una rejilla y normalizar:

import numpy as np
from scipy import stats
rng = np.random.default_rng(42)            # mismos datos que make_sueno()
u0 = rng.normal(0, 25, 18); rng.normal(0, 6, 18)
y = (250 + u0 + rng.normal(0, 25, (18, 10))[:, 0]).round(1)  # día 0
mu = np.linspace(150, 400, 2001)            # rejilla de valores de mu
log_post = stats.norm.logpdf(mu, 300, 50)   # prior N(300, 50)
log_post += stats.norm.logpdf(y[:, None], mu, 35).sum(axis=0)  # verosimilitud
post = np.exp(log_post - log_post.max())
post /= post.sum()
cdf = np.cumsum(post)
print("media:", (mu * post).sum().round(2))
print("IC 94%:", mu[np.searchsorted(cdf, [0.03, 0.97])].round(1))
print("P(mu > 260):", post[mu > 260].sum().round(3))

La media posterior es 250.64 ms, el intervalo de credibilidad al 94% va de 235.4 a 266.0 y $P(\mu > 260 \mid y) = 0.124$. El prior pedía 300, pero 18 observaciones lo dominan: la media muestral es 249.3. Se trabaja en escala logarítmica para evitar que el producto de densidades se redondee a cero.

El caso conjugado y la predictiva posterior

Con prior normal y verosimilitud normal de $\sigma$ conocida, el posterior también es normal y las precisiones (inversos de varianzas) se suman: $1/s_n^2 = 1/s_0^2 + n/\sigma^2$. La predictiva posterior responde otra pregunta: qué valor tendrá una observación nueva. Integra la incertidumbre del parámetro y el ruido de los datos, $p(\tilde y \mid y) = \int p(\tilde y \mid \mu)\, p(\mu \mid y)\, d\mu$, y se aproxima por simulación:

import numpy as np
rng = np.random.default_rng(42)
u0 = rng.normal(0, 25, 18); rng.normal(0, 6, 18)
y = (250 + u0 + rng.normal(0, 25, (18, 10))[:, 0]).round(1)
n, sigma, m0, s0 = len(y), 35, 300, 50
prec = 1 / s0**2 + n / sigma**2                 # las precisiones se suman
m_post = (m0 / s0**2 + y.sum() / sigma**2) / prec
s_post = prec ** -0.5
print("posterior: N(%.2f, %.2f)" % (m_post, s_post))
mu_draws = rng.normal(m_post, s_post, 4000)     # muestras del posterior
y_nuevo = rng.normal(mu_draws, sigma)           # predictiva posterior
print("y nuevo, 94%:", np.percentile(y_nuevo, [3, 97]).round(1))
print("P(y nuevo > 300):", (y_nuevo > 300).mean().round(3))

El resultado exacto es $\mathcal{N}(250.64, 8.14)$, igual que la rejilla. El intervalo predictivo va de 183.1 a 318.4 ms, mucho más ancho que el de $\mu$: un sujeto nuevo varía más que la media. Es la base de los chequeos con az.plot_ppc.

Credibilidad contra confianza

Un intervalo de credibilidad al 94% dice: dado el modelo y los datos, $\mu$ está en el intervalo con probabilidad 0.94. Un intervalo de confianza al 95% dice otra cosa: si repitieras el experimento muchas veces, el 95% de los intervalos construidos así contendría el valor verdadero; sobre un intervalo concreto no asigna probabilidad. ArviZ reporta por defecto el HDI al 94% (columnas hdi_3% y hdi_97%), el intervalo más corto que contiene esa masa.

Por qué hacen falta muestreadores

La rejilla funciona con uno o dos parámetros. El modelo de pendientes aleatorias que ajustarás con Bambi tiene 41: intercepto, pendiente, $\sigma$, dos desviaciones de grupo y 36 desviaciones por sujeto. Con 100 puntos por eje serían $100^{41}$ evaluaciones. Los métodos MCMC (Markov chain Monte Carlo) construyen una cadena cuyas muestras se distribuyen, a la larga, según el posterior, sin calcular $p(y)$. NUTS (No-U-Turn Sampler), el muestreador por defecto de PyMC, usa el gradiente del log-posterior para dar pasos largos en direcciones útiles y ajusta solo la longitud de cada trayectoria. Bambi traduce una fórmula como tiempo_reaccion ~ dias + (dias|sujeto) a un modelo de PyMC con priors débilmente informativos y lanza el muestreo; en el curso usarás nutpie, una implementación de NUTS más rápida.

Trampas comunes

• Leer el posterior como un punto. Reporta siempre un intervalo.
• Confundir el intervalo de $\mu$ con el de una observación nueva. Para predecir individuos usa la predictiva posterior.
• Priors planos «para no influir». En modelos jerárquicos dificultan el muestreo; uno débilmente informativo es más estable.
• Confiar en MCMC sin diagnosticar. Revisa $\hat R$ y el tamaño efectivo de muestra.

Cierre

Antes de ajustar un modelo, escribe el prior de cada parámetro, la verosimilitud y la pregunta que responderá el posterior. Bambi automatiza el cálculo; esas decisiones siguen siendo tuyas.