Inferencia bayesiana y MCMC

1. Motivación: del mejor ajuste a la distribución de los parámetros

El periodograma nos dio un periodo candidato; ahora queremos los parámetros orbitales con sus incertidumbres y correlaciones: ¿cuánto vale $K$ y con qué error? ¿Está $e$ realmente distinta de cero? ¿Hay un segundo planeta? Mínimos cuadrados responde a la primera pregunta solo si el modelo es casi lineal y los errores gaussianos; la inferencia bayesiana responde a todas y lo hace con un único objeto, la distribución posterior de los parámetros, que explora con cadenas de Markov. En exoplanetas es la metodología estándar (Gregory 2005; Foreman-Mackey et al. 2013), y se aplica igual a curvas de luz de tránsito, a cosmología o a ondas gravitacionales.

2. Teorema de Bayes y verosimilitud gaussiana

Para parámetros $\boldsymbol\theta$ y datos $D$,

con la verosimilitud $\mathcal L(\boldsymbol\theta) = p(D|\boldsymbol\theta)$, el prior $p(\boldsymbol\theta)$ y la evidencia $p(D)$ (una constante para estimar parámetros; esencial para comparar modelos). Si cada dato tiene error gaussiano independiente de desviación $\sigma_i$ conocida y el modelo predice $m_i(\boldsymbol\theta)$,

Con prior uniforme, el máximo a posteriori coincide con el mínimo de $\chi^2$: mínimos cuadrados es el caso particular bayesiano de priors planos y errores gaussianos (problema 1). Si sospechamos que los errores están subestimados (ruido estelar), se añade un término de jitter $s$: $\sigma_i^2 \to \sigma_i^2 + s^2$, con $s$ un parámetro más (y entonces $\sum\ln\sigma_i$ ya no es constante).

Priors

El prior codifica lo que sabemos antes de ver los datos y debe declararse. Para un parámetro de posición (como $\gamma$ o $t_p$) el prior uniforme es natural; para un parámetro de escala positivo que abarca órdenes de magnitud ($P$, $K$) es más honesto el prior de Jeffreys $p(\theta) \propto 1/\theta$ (uniforme en $\ln\theta$), que no favorece ninguna década (problema 6). Para $e$ y $\omega$ se suele muestrear en $(\sqrt e\cos\omega, \sqrt e\sin\omega)$, que da prior uniforme en $e$ y evita la singularidad de $\omega$ en $e = 0$. Si el prior es informativo (p. ej. $M_\star$ de una isocrona) se incluye como gaussiana.

3. Metropolis–Hastings

Queremos muestras $\boldsymbol\theta^{(1)}, \boldsymbol\theta^{(2)}, \ldots$ distribuidas según $\pi(\boldsymbol\theta) \propto \mathcal L(\boldsymbol\theta)p(\boldsymbol\theta)$ sin conocer la constante de normalización. El algoritmo de Metropolis–Hastings (Metropolis et al. 1953; Hastings 1970):

1. Desde el estado actual $\mathbf x$, proponer $\mathbf x'$ con una densidad de propuesta $q(\mathbf x'|\mathbf x)$.
2. Calcular $\alpha = \min\!\left\{1,\ \dfrac{\pi(\mathbf x')\,q(\mathbf x|\mathbf x')}{\pi(\mathbf x)\,q(\mathbf x'|\mathbf x)}\right\}$.
3. Con probabilidad $\alpha$ aceptar ($\mathbf x^{(k+1)} = \mathbf x'$); si no, repetir el estado ($\mathbf x^{(k+1)} = \mathbf x$).

Con propuesta gaussiana simétrica (paseo aleatorio), $q(\mathbf x'|\mathbf x) = q(\mathbf x|\mathbf x')$ y $\alpha = \min\{1, \pi(\mathbf x')/\pi(\mathbf x)\}$: siempre se sube, y se baja con probabilidad igual al cociente. En logaritmos: aceptar si $\ln u < \ln\pi(\mathbf x') - \ln\pi(\mathbf x)$ con $u \sim U(0,1)$. Es lo que implementa metropolisHastings en src/lib/math/mcmc.ts.

Por qué funciona: balance detallado

Teorema. La cadena de Metropolis–Hastings tiene a $\pi$ como distribución estacionaria.

Demostración. El núcleo de transición para $\mathbf x' \ne \mathbf x$ es $T(\mathbf x'|\mathbf x) = q(\mathbf x'|\mathbf x)\,\alpha(\mathbf x'|\mathbf x)$. Supongamos $\pi(\mathbf x')q(\mathbf x|\mathbf x') \le \pi(\mathbf x)q(\mathbf x'|\mathbf x)$ (el otro caso es simétrico). Entonces $\alpha(\mathbf x'|\mathbf x) = \pi(\mathbf x')q(\mathbf x|\mathbf x')/[\pi(\mathbf x)q(\mathbf x'|\mathbf x)]$ y $\alpha(\mathbf x|\mathbf x') = 1$, de modo que

Esta igualdad es el balance detallado. Integrando en $\mathbf x$: $\int\pi(\mathbf x)T(\mathbf x'|\mathbf x)\,d\mathbf x = \pi(\mathbf x')\int T(\mathbf x|\mathbf x')\,d\mathbf x = \pi(\mathbf x')$, es decir, si el estado actual se distribuye según $\pi$, el siguiente también: $\pi$ es estacionaria. (Los rechazos, que dejan $\mathbf x' = \mathbf x$, satisfacen el balance trivialmente.) $\blacksquare$

Que la cadena converja a $\pi$ desde cualquier punto exige además irreducibilidad y aperiodicidad, garantizadas por una propuesta gaussiana con soporte en todo el espacio. Lo que el teorema no dice es cuánto tarda: de eso tratan los diagnósticos.

La escala de la propuesta

Con pasos muy pequeños casi todo se acepta pero la cadena avanza como un paseo aleatorio lentísimo; con pasos muy grandes casi todo se rechaza. Para blancos gaussianos de dimensión $d$ la tasa de aceptación óptima es $\approx 0.23$ (0.44 en $d = 1$), con escala $\approx 2.4\,\sigma_{\rm post}/\sqrt d$ (Roberts, Gelman & Gilks 1997). En nuestro ejemplo (§5), multiplicar la escala por 0.2, 1 y 5 da aceptaciones de 0.90, 0.56 y 0.045, y tiempos de autocorrelación de $K$ de 203, 23 y 36 pasos: el centro gana.

4. Diagnósticos de convergencia

Gelman–Rubin

Se corren $m$ cadenas desde puntos iniciales dispersos, de longitud $n$ tras descartar el calentamiento (burn-in). Para cada parámetro escalar, sean $W$ la media de las varianzas intracadena y $B = n\,\mathrm{Var}(\text{medias de las cadenas})$ la varianza entre cadenas. El estimador de la varianza de la posterior

sobrestima $\mathrm{Var}(\theta)$ si las cadenas no se han mezclado (cada una explora una región distinta: $B$ grande), mientras que $W$ la subestima; ambos convergen y $\hat R \to 1$. Criterio práctico: $\hat R < 1.1$ (1.01 en trabajos exigentes). Un $\hat R \approx 1$ no demuestra convergencia (todas las cadenas pueden estar atrapadas en el mismo modo), pero $\hat R \gg 1$ sí demuestra lo contrario. Implementación: gelmanRubin.

Autocorrelación y tamaño efectivo

Las muestras consecutivas están correlacionadas. Para una cadena $x_k$ con función de autocorrelación normalizada $\rho(\ell)$, el tiempo de autocorrelación integrado es

y la varianza de la media muestral es $\mathrm{Var}(\bar x) = \tau\,\sigma^2/n$: una cadena de $n$ pasos vale como $n_{\rm eff} = n/\tau$ muestras independientes. La suma se trunca con la ventana automática de Sokal (parar en el primer $M \ge c\tau(M)$, $c = 5$; es lo que hace emcee) porque los $\rho(\ell)$ de retardo grande son puro ruido. Para un proceso AR(1) con coeficiente $\phi$, $\tau = (1+\phi)/(1-\phi)$ (test con $\phi = 0.9$: $\tau = 19$). Regla práctica: correr al menos $50\tau$ pasos por cadena, y descartar unos pocos $\tau$ como calentamiento. Implementación: tiempoAutocorrelacion, tamanoEfectivo.

emcee: el muestreador de conjunto

El paseo aleatorio requiere elegir una escala por parámetro, y falla con correlaciones fuertes. Goodman & Weare (2010) propusieron evolucionar un conjunto de $L$ caminantes donde cada uno se mueve a lo largo de la recta que lo une con otro caminante elegido al azar («stretch move», $\mathbf x' = \mathbf x_j + Z(\mathbf x - \mathbf x_j)$ con $Z$ de densidad $\propto 1/\sqrt Z$ en $[1/a, a]$, $a = 2$), aceptando con probabilidad $\min\{1, Z^{d-1}\pi(\mathbf x')/\pi(\mathbf x)\}$. El algoritmo es invariante afín: su rendimiento es el mismo para una gaussiana muy alargada y para una esférica, sin ajuste manual. emcee (Foreman-Mackey et al. 2013) es la implementación de referencia en Python, con estimación de $\tau$ incorporada; en esta lección usamos Metropolis–Hastings por transparencia, pero en la práctica profesional se usa emcee o muestreo anidado (dynesty, MultiNest) cuando además se quiere la evidencia.

5. Ejemplo: la curva de velocidad radial de 51 Peg

Datos: las mismas 40 medidas sintéticas de la lección anterior ($P = 4.2293$ d, $K = 59$ m s⁻¹, $\gamma = -33\,250$ m s⁻¹, $\sigma = 13$ m s⁻¹, 61 noches). Modelo circular de cuatro parámetros $\boldsymbol\theta = (P, K, \gamma, t_0)$:

con priors uniformes $P \in [4.0, 4.5]$ d, $K \in [0, 200]$, $\gamma \in [-34\,000, -32\,000]$, $t_0 \in [-5, 5]$ d. Cuatro cadenas de 40 000 pasos desde puntos dispersos, escalas de propuesta $(0.002, 2, 2, 0.01)$, calentamiento 5000. Resultados (mcmc.ts, semillas fijas):

Aceptación 0.56 en las cuatro cadenas. Comentarios:

• Comparación con la teoría lineal. Para una sinusoide bien muestreada, $\sigma_K = \sigma\sqrt{2/N} = 2.91$ m s⁻¹ y $\sigma_\gamma = \sigma/\sqrt N = 2.06$: la posterior reproduce ambos (3.0 y 2.1). $K$ cae a 0.9 σ del valor verdadero: una realización de ruido normal.
• Degeneración $P$–$t_0$. Cambiar $P$ en $\delta P$ desplaza la fase en el centro de la campaña en $2\pi(30\ \mathrm d)\delta P/P^2$, que se compensa con $t_0$: la posterior es una «banana» alargada y la cadena de paseo aleatorio la recorre despacio ($\tau \approx 260$ frente a 10–20 en los otros). Con 35 000 pasos útiles, $n_{\rm eff} \approx 135$ por cadena para $P$: suficiente para la media (error de Monte Carlo $\sigma_P/\sqrt{4\times135} = 0.0003$ d) pero escaso para colas. Es el caso típico en que emcee (invariante afín) o reparametrizar ($t_0$ referido al centro de la campaña) gana mucho.
• Precisión del periodo. $\sigma_P = 0.008$ d es 25 veces menor que la resolución del periodograma ($0.03$ d): el periodograma localiza, el ajuste refina.

Comparación de modelos

¿Hay planeta? Comparamos el modelo constante ($k = 1$ parámetro, $\chi^2 = 400.5$) con el circular ($k = 4$, $\chi^2 = 48.2$ en el máximo). El criterio de información bayesiano $\mathrm{BIC} = \chi^2 + k\ln N$ da $\Delta\mathrm{BIC} = -352.3 + 3\ln 40 = -341$: evidencia abrumadora (umbral habitual $|\Delta\mathrm{BIC}| > 10$). El BIC es una aproximación a $-2\ln p(D)$ válida para $N$ grande y posteriores gaussianas; la comparación rigurosa usa la evidencia, que penaliza automáticamente la complejidad (factor de Occam) y que calculan el muestreo anidado o la integración térmica.

6. Análisis de error

• Error de Monte Carlo. La incertidumbre de cualquier promedio de la cadena es $\sigma_{\rm post}\sqrt{\tau/n}$; para percentiles extremos (p. ej. el límite al 99 % de $e$) hace falta $n_{\rm eff} \gtrsim 10^3$–$10^4$.
• Calentamiento. Si se incluye, sesga las medias hacia el punto inicial; se detecta en la traza y con $\hat R$ sobre tramos.
• Multimodalidad. Si el periodograma tiene alias comparables (lección anterior), la posterior tiene varios modos y el paseo aleatorio no salta entre ellos: hay que lanzar cadenas desde cada alias o usar muestreo anidado.
• Modelo erróneo. Ningún muestreador corrige un modelo incompleto (segundo planeta, actividad estelar): se ve en los residuos y en un jitter ajustado mucho mayor que $\sigma_i$.

Problemas

1. (Demostración) Mínimos cuadrados como caso bayesiano. Demuestra que con errores gaussianos independientes y prior uniforme el máximo a posteriori minimiza $\chi^2$, y que si el modelo es lineal, $m = \mathbf X\boldsymbol\theta$, la posterior es gaussiana con media $(\mathbf X^T\mathbf C^{-1}\mathbf X)^{-1}\mathbf X^T\mathbf C^{-1}\mathbf y$ y covarianza $(\mathbf X^T\mathbf C^{-1}\mathbf X)^{-1}$, con $\mathbf C = \mathrm{diag}(\sigma_i^2)$.

2. (Demostración) Hastings con propuesta asimétrica. Demuestra que si se usa una propuesta no simétrica (p. ej. lognormal para un parámetro positivo) pero se acepta con $\min\{1, \pi(\mathbf x')/\pi(\mathbf x)\}$ sin el cociente $q(\mathbf x|\mathbf x')/q(\mathbf x'|\mathbf x)$, la distribución estacionaria no es $\pi$. Calcula cuál es para $q(x'|x) = \mathcal N(x + \delta, s^2)$ con deriva $\delta$ y blanco uniforme en una caja grande.

3. Gelman–Rubin a mano. Tres cadenas de $n = 1000$ pasos para $K$ tienen medias 55.1, 57.8 y 56.4 y varianzas 8.9, 9.4 y 9.1. Calcula $W$, $B$, $\hat V$ y $\hat R$. ¿Han convergido? ¿Y si las medias fueran 50, 57 y 63?

4. (Computacional) Escala de la propuesta. Con el ejemplo de §5 (o con una gaussiana 2D correlacionada, $\rho = 0.95$), corre Metropolis con escalas $0.1\sigma$, $\sigma$, $2.4\sigma/\sqrt d$ y $10\sigma$. Tabula aceptación, $\tau$ y $n_{\rm eff}$ y verifica que el máximo de $n_{\rm eff}$ cae cerca de aceptación 0.2–0.4.

import numpy as np
rng = np.random.default_rng(0)
C = np.array([[1, 0.95], [0.95, 1]]); Ci = np.linalg.inv(C)
logp = lambda x: -0.5 * x @ Ci @ x
def tau_int(x, c=5):
    x = x - x.mean(); n = len(x); acf = np.correlate(x, x, 'full')[n-1:] / (x @ x)
    t = 1.0
    for M in range(1, n // 2):
        t += 2 * acf[M]
        if M >= c * t: return t
    return t
for s in (0.1, 1.0, 2.4 / np.sqrt(2), 10.0):
    x = np.zeros(2); lp = logp(x); X = []; acc = 0
    for _ in range(50000):
        xp = x + s * rng.standard_normal(2); lpp = logp(xp)
        if np.log(rng.random()) < lpp - lp: x, lp, acc = xp, lpp, acc + 1
        X.append(x.copy())
    X = np.array(X)[5000:]; t = tau_int(X[:, 0])
    print(s, acc / 50000, t, len(X) / t)

5. (Computacional) Ajuste completo. Reproduce §5 con tu propio Metropolis (o con emcee): cuatro cadenas, trazas, $\hat R$, $\tau$, histogramas marginales y el diagrama $P$–$t_0$. Reparametriza con $t_0' = t_0 - 30$ d (centro de la campaña) y compara $\tau_P$. Añade después el parámetro de jitter $s$ y comprueba que su posterior es compatible con cero.

6. (Demostración) Prior de Jeffreys. Demuestra que $p(\theta) \propto 1/\theta$ es el único prior en $\theta > 0$ invariante ante cambios de unidades $\theta \to c\theta$, y que equivale a un prior uniforme en $\ln\theta$. Para $K \in [0.1, 1000]$ m s⁻¹, ¿qué fracción de la masa del prior cae por debajo de 1 m s⁻¹ con el prior uniforme y con el de Jeffreys?

7. Occam. Para los datos de §5, un modelo con excentricidad libre ($k = 6$) alcanza $\chi^2 = 46.9$. Calcula $\Delta\mathrm{BIC}$ frente al circular y decide. ¿Qué dice el teorema de Wilks sobre $\Delta\chi^2 = 1.3$ con 2 grados de libertad extra?

Lecturas

• Gregory, P. C. (2005). Bayesian Logical Data Analysis for the Physical Sciences, cap. 3 (fundamentos), cap. 11 (ajuste no lineal) y cap. 12 (MCMC, con el ejemplo de velocidad radial).
• Hogg & Foreman-Mackey (2018). Data analysis recipes: Using MCMC, ApJS 236, 11: guía práctica, breve y honesta sobre qué funciona.
• Foreman-Mackey et al. (2013), PASP 125, 306 (emcee) y Goodman & Weare (2010): el muestreador invariante afín.
• Gelman & Rubin (1992), Stat. Sci. 7, 457; Sokal (1997): diagnósticos de convergencia y autocorrelación.
• Perryman (2018). The Exoplanet Handbook, §2.3 y §6.4: inferencia de órbitas por velocidad radial y tránsitos.