Metropolis-Hastings desde cero

Estadística · Lección 22

L3 Derivación Fase F5 Libro Think Bayes 2e §19 Prereqs Priors conjugados · Bootstrap

Objetivo

Al terminar esta lección se puede:

1. Definir cadena de Markov, distribución estacionaria y balance detallado, y demostrar que el balance detallado implica la estacionariedad.
2. Demostrar que el núcleo de Metropolis-Hastings cumple el balance detallado respecto del objetivo, y deducir que la constante de normalización se cancela: basta conocer el posterior salvo un factor.
3. Escribir el algoritmo en diez líneas y verificarlo contra un posterior conjugado de respuesta conocida.
4. Demostrar cómo la autocorrelación de la cadena infla la varianza de sus promedios, calcular el tamaño de muestra efectivo y usarlo para elegir el paso.

De dónde viene

• Estadística 21: la verosimilitud marginal, la integral que no se sabe calcular. Este método existe porque no hace falta calcularla para muestrear el posterior.
• Estadística 19: los casos conjugados. Aquí sirven de banco de pruebas: se aplica el método a un problema cuya respuesta ya se sabe y se comprueba que la reproduce.
• Estadística 11: el bootstrap ya resolvía por simulación lo que no salía en forma cerrada. La idea es la misma; lo que cambia es que aquí las muestras están correlacionadas, y eso cuesta.
• Estadística 5: la correlación entre observaciones sucesivas, que en la sección 4 es lo que decide cuántas muestras hacen falta.

Para qué sirve después

• Lección 23: Gibbs es un caso particular en el que la propuesta es la condicional exacta y la aceptación vale siempre uno.
• Lección 24: los diagnósticos de cadena responden a la pregunta que esta lección deja abierta, que es cuándo parar.
• Lección 25: PyMC hace esto por dentro con propuestas mejores; saber qué hace evita tratarlo como una caja negra.
• Estadística 12: en más de una dimensión la forma del posterior —su matriz de covarianza— decide cuánto cuesta recorrerlo, y ahí la lección 12 vuelve a hacer falta.

Notación

Símbolo	Se lee	Significado
\pi(\theta)	pi de theta	La distribución objetivo, aquí el posterior; se conoce salvo constante
\tilde\pi(\theta)	pi tilde de theta	Una función con \pi=\tilde\pi/Z: el posterior sin normalizar
Z	zeta	La constante desconocida, \int\tilde\pi
q(y\mid x)	cu de ye dado equis	La densidad de propuesta: desde dónde se salta y adónde
\alpha(x,y)	alfa	La probabilidad de aceptar el salto de x a y
P(x,y)	pe de equis, ye	El núcleo de transición de la cadena
\rho_k	rho sub ka	La autocorrelación de la cadena a retardo k
\tau	tau	El factor de ineficiencia, 1+2\sum_k\rho_k

1. El problema

Las lecciones 19 y 20 vivían de la conjugación. Fuera de ella, el posterior se conoce solo así:

\pi(\theta\mid x)=\frac{p(x\mid\theta)\,\pi(\theta)}{Z},\qquad Z=\int p(x\mid\theta)\,\pi(\theta)\,d\theta,

con Z imposible de calcular en cuanto \theta tiene más de dos o tres componentes: una rejilla de cien puntos por dimensión son 10^{10} evaluaciones en cinco dimensiones y 10^{20} en diez. El numerador, en cambio, se evalúa en un instante para cualquier \theta.

La pregunta que resuelve esta lección es si eso basta. Y basta: se puede muestrear de \pi conociendo solo \tilde\pi=p(x\mid\theta)\pi(\theta), sin tocar Z. Con muestras, cualquier esperanza se aproxima por un promedio, y el posterior queda descrito.

2. Cadenas y balance detallado

Definición 22.1 (cadena de Markov). Una sucesión \theta_0,\theta_1,\dots es una cadena de Markov con núcleo P si la distribución de \theta_{t+1} depende del pasado solo a través de \theta_t: P(x,A)=\Pr(\theta_{t+1}\in A\mid \theta_t=x).

Definición 22.2 (distribución estacionaria). Una distribución \pi es estacionaria para P si al arrancar en \pi se sigue en \pi: en el caso discreto, \sum_x \pi(x)P(x,y)=\pi(y) para todo y.

Definición 22.3 (balance detallado). El núcleo P cumple el balance detallado respecto de \pi si \pi(x)\,P(x,y)=\pi(y)\,P(y,x)\qquad\text{para todo par }x,y. Es decir, en equilibrio pasa tanta masa de x a y como de y a x.

Proposición 22.4 (el balance detallado basta). Si P cumple el balance detallado respecto de \pi, entonces \pi es estacionaria para P.

Demostración. Sumando la igualdad del balance detallado sobre x, \sum_x \pi(x)P(x,y)=\sum_x \pi(y)P(y,x)=\pi(y)\sum_x P(y,x)=\pi(y), porque las probabilidades de transición desde y suman uno. Eso es exactamente la Definición 22.2. ∎

El recíproco es falso, y conviene saberlo: hay cadenas con distribución estacionaria que no cumplen el balance detallado. Lo que la proposición ofrece es una condición suficiente y, sobre todo, local: se comprueba mirando pares de estados, sin sumar sobre nada. Eso es lo que hace posible diseñar un núcleo a medida.

3. El algoritmo

Definición 22.5 (Metropolis-Hastings). Dada una función \tilde\pi>0 proporcional al objetivo y una densidad de propuesta q(\cdot\mid\cdot), desde el estado x:

1. Proponer y\sim q(\cdot\mid x).
2. Calcular la razón de aceptación \alpha(x,y)=\min\left(1,\;\frac{\tilde\pi(y)\,q(x\mid y)}{\tilde\pi(x)\,q(y\mid x)}\right).
3. Con probabilidad \alpha(x,y) moverse a y; en otro caso quedarse en x y anotar x otra vez.

El tercer punto es donde se equivoca casi todo el que lo implementa por primera vez: un rechazo no se descarta, se repite. Si se descartara, la cadena no tendría la distribución estacionaria que la proposición siguiente le atribuye.

Proposición 22.6 (Metropolis-Hastings cumple el balance detallado). El núcleo de la Definición 22.5 cumple el balance detallado respecto de \pi=\tilde\pi/Z, cualquiera que sea Z>0.

Demostración. Para x\ne y el núcleo es P(x,y)=q(y\mid x)\,\alpha(x,y), porque llegar a y exige proponerlo y aceptarlo. Sustituyendo esa expresión en la Definición 22.3, lo que hay que probar es \pi(x)\,q(y\mid x)\,\alpha(x,y)=\pi(y)\,q(x\mid y)\,\alpha(y,x). Escríbase r=\dfrac{\pi(y)\,q(x\mid y)}{\pi(x)\,q(y\mid x)}, de modo que \alpha(x,y)=\min(1,r) y \alpha(y,x)=\min(1,1/r). Se distinguen dos casos.

Si r\le1: entonces \alpha(x,y)=r y \alpha(y,x)=1, y el lado izquierdo vale \pi(x)\,q(y\mid x)\cdot\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}=\pi(y)\,q(x\mid y), que es el lado derecho.

Si r>1: entonces \alpha(x,y)=1 y \alpha(y,x)=1/r, y ahora es el lado derecho el que vale \pi(y)\,q(x\mid y)\cdot\frac{\pi(x)q(y\mid x)}{\pi(y)q(x\mid y)}=\pi(x)\,q(y\mid x), que es el lado izquierdo. Para x=y la igualdad es trivial. Por la Proposición 22.4, \pi es estacionaria. ∎

Corolario 22.7 (la constante se cancela). La razón \alpha de la Definición 22.5 no cambia si se sustituye \tilde\pi por c\,\tilde\pi con c>0. En particular el algoritmo se puede correr con \tilde\pi=p(x\mid\theta)\pi(\theta) sin conocer Z.

Demostración. \tilde\pi aparece en \alpha solo como el cociente \tilde\pi(y)/\tilde\pi(x), y ese cociente es invariante al multiplicar por c. ∎

Ahí está todo. La integral de la lección 21, la que no se sabe hacer, entra en la fórmula dividida por sí misma. El método no la calcula: la evita.

Proposición 22.8 (caso simétrico). Si la propuesta es simétrica, q(y\mid x)=q(x\mid y) —por ejemplo y=x+\varepsilon con \varepsilon de una densidad simétrica alrededor de cero—, entonces \alpha(x,y)=\min\left(1,\;\frac{\tilde\pi(y)}{\tilde\pi(x)}\right).

Demostración. Los factores q del numerador y del denominador de la Definición 22.5 son iguales y se cancelan. ∎

Con eso la regla se lee sola: si la propuesta sube, se acepta siempre; si baja, se acepta con probabilidad igual a cuánto baja. La cadena sube sin condiciones y baja con reticencia, y ese desequilibrio es exactamente el que reproduce \pi.

Comprobarlo con la matriz entera

En un espacio de estados pequeño no hace falta simular: el núcleo se puede escribir como matriz y verificar la Proposición 22.6 con álgebra lineal.

La cadena olvida dónde empezó —todas las filas de P^{400} son iguales— y lo que recuerda es \pi. Y en ningún momento se usó la suma 26 que normaliza el objetivo: solo cocientes.

El caso continuo, contra una respuesta conocida

El visual siguiente deja mover el paso de la propuesta y enseña las dos cosas que decide: el recorrido de la cadena y el histograma que produce, contra la densidad exacta. El marco no depende de la cadena: el eje del histograma es [0,1] porque \theta es una probabilidad, y la altura tope es constante.

Paso de la propuesta 0,30 Pasos de cadena 1200

4. Lo que la Proposición 22.6 no promete

La proposición dice que \pi es estacionaria. No dice que la cadena llegue a \pi desde donde se arranque —para eso hace falta que sea irreducible y aperiódica—, ni dice cuándo. Y, sobre todo, no dice que las muestras sean independientes. Lo son las del bootstrap de la lección 11; estas, no, y eso tiene un precio exacto.

Proposición 22.9 (el precio de la correlación). Sea \theta_1,\dots,\theta_n una cadena estacionaria con varianza \sigma^2 y autocorrelaciones \rho_k=\text{Corr}(\theta_t,\theta_{t+k}). Entonces \text{Var}(\bar\theta_n)=\frac{\sigma^{2}}{n}\left(1+2\sum_{k=1}^{n-1}\Big(1-\frac{k}{n}\Big)\rho_k\right)\;\xrightarrow[n\to\infty]{}\;\frac{\sigma^{2}}{n}\,\tau,\qquad \tau=1+2\sum_{k=1}^{\infty}\rho_k, suponiendo que la serie converge. El tamaño de muestra efectivo es n_{\text{ef}}=n/\tau: la cadena da la misma precisión que n_{\text{ef}} muestras independientes.

Demostración. Por bilinealidad de la covarianza, \text{Var}\Big(\frac1n\sum_{t=1}^n\theta_t\Big)=\frac{1}{n^{2}}\sum_{t=1}^{n}\sum_{s=1}^{n}\text{Cov}(\theta_t,\theta_s) =\frac{\sigma^{2}}{n^{2}}\sum_{t,s}\rho_{|t-s|}. En la matriz de índices, el retardo k aparece exactamente n-|k| veces para k=-(n-1),\dots,n-1. Agrupando por retardo y usando \rho_{-k}=\rho_k, \frac{\sigma^{2}}{n^{2}}\Big[n\rho_0+2\sum_{k=1}^{n-1}(n-k)\rho_k\Big]=\frac{\sigma^{2}}{n}\Big[1+2\sum_{k=1}^{n-1}\Big(1-\frac{k}{n}\Big)\rho_k\Big], que es el enunciado. El límite es convergencia dominada sobre la serie, y la definición de n_{\text{ef}} sale de igualar \sigma^2\tau/n con \sigma^2/n_{\text{ef}}. ∎

Con \rho_k\ge0 —lo habitual en una caminata aleatoria— se tiene \tau\ge1, y el tamaño efectivo nunca supera al nominal. Un \tau de 150 significa que de cada 150 iteraciones se saca la información de una.

El paso de la propuesta controla \tau, y lo hace por dos caminos opuestos. Un paso pequeño acepta casi siempre pero avanza poco, así que \rho_1 es alta; un paso grande propone lejos pero se rechaza casi siempre, y quedarse quieto también es correlación. El óptimo está en medio y se puede medir.

El mejor de los cinco pasos es 0{,}30, y su tasa de aceptación cae en 0{,}44. No es casualidad: para una caminata aleatoria en una dimensión, la tasa de aceptación que maximiza el tamaño efectivo es conocida y vale aproximadamente 0{,}44. Los dos extremos pierden por motivos opuestos: con paso 0{,}02 se acepta el 95 % y \tau pasa de cien; con paso 3{,}00 se acepta el 5 % y \tau vuelve a dispararse.

El visual siguiente barre el paso entero en lugar de mirar cinco valores sueltos. Corre una cadena por cada paso del barrido y dibuja las dos curvas que importan: la fracción de muestra efectiva y la tasa de aceptación.

Pasos por cadena 1500 Forma del objetivo Beta(9,5) Beta(2,2), ancha Beta(40,20), estrecha

5. Qué mirar antes de creerse una cadena

Tres cosas, y la lección 24 las convierte en números.

La primera es el arranque. La cadena empieza donde se la ponga y tarda en olvidarlo; los primeros estados no son muestras de \pi y se descartan. En las celdas de arriba se tiraron cinco mil de sesenta mil, por exceso de prudencia y sin justificarlo.

La segunda es la tasa de aceptación, que es gratis de calcular y detecta los dos desastres: cerca de uno significa paso demasiado corto, cerca de cero paso demasiado largo.

La tercera es el tamaño efectivo. Una cadena de un millón de iteraciones con \tau=10^4 tiene cien muestras útiles, y reportar un intervalo como si tuviera un millón es la forma más común de mentir con MCMC sin querer.

Ejercicios

1. Demostrar que si q(y\mid x)>0 siempre que \pi(y)>0, entonces la cadena de la Definición 22.5 es irreducible. ¿Qué pasa si la propuesta solo puede moverse dentro de un intervalo que no cubre todo el soporte?
2. Comprobar que el algoritmo de Metropolis-Hastings con q(y\mid x)=\pi(y) —proponer directamente del objetivo— acepta siempre. ¿Por qué no sirve eso de nada en la práctica?
3. Demostrar que el balance detallado no es necesario: construir una cadena de tres estados con distribución estacionaria uniforme que lo viole. (Pista: un ciclo que gira siempre en el mismo sentido.)
4. En la Proposición 22.9, suponer \rho_k=\phi^{k} con 0<\phi<1 y demostrar que \tau=(1+\phi)/(1-\phi). Calcular el \phi que corresponde a \tau=150.
5. Demostrar que si se descartan los rechazos en vez de repetir el estado, la distribución estacionaria deja de ser \pi. Exhibir el sesgo en el ejemplo de seis estados de la sección 3.

Reto

En proyectos/notebooks/F1-retos.ipynb, sección Est 22:

1. Escribir metropolis(log_objetivo, x0, paso, n, semilla) que devuelva la cadena y la tasa de aceptación, y tau_ess(cadena). Verificar contra tres posteriores conjugados distintos —Beta, Gamma y normal— que la media y la desviación de la cadena caen dentro de tres errores estándar del valor exacto.
2. Reproducir la tabla de la sección 4 con un barrido fino de pasos y localizar numéricamente el que maximiza el tamaño efectivo. Comprobar que la tasa de aceptación en ese punto se acerca a 0{,}44, y repetir en dos dimensiones para ver hacia dónde se mueve.
3. Implementar la comprobación del ejercicio 5: correr la versión incorrecta —descartar los rechazos— sobre el espacio de seis estados, estimar su distribución estacionaria y medir la distancia a \pi. Explicar con la Proposición 22.6 dónde se rompe la demostración.

Del libro

Think Bayes se publica bajo licencia CC BY-NC-SA 4.0, que permite citarlo textualmente. El capítulo de MCMC publica sus secciones sin numerar en la versión en línea, así que aquí se cita a nivel de capítulo.

MCMC, which stands for Markov chain Monte Carlo.

Downey, Think Bayes 2e, §19 MCMC

El capítulo 19 usa PyMC desde la primera página y no escribe el algoritmo. Esa es la diferencia de camino: esta lección lo deriva del balance detallado y lo implementa en diez líneas, para que la lección 25 pueda usar PyMC sabiendo qué hay debajo.

Preguntas para leer §19 con lápiz:

1. El capítulo presenta MCMC como alternativa a la rejilla. ¿A partir de cuántas dimensiones deja de caber una rejilla de cien puntos por eje, y cómo se compara ese número con el coste de una cadena?
2. PyMC pide especificar el modelo pero no la propuesta. ¿Qué parte de la Definición 22.5 está eligiendo la librería por su cuenta, y qué riesgo se corre al no verla?
3. El capítulo muestra trazas de las cadenas. ¿Qué de lo que se ve en esas trazas corresponde al \tau de la Proposición 22.9, y qué habría que medir para no depender del ojo?

Para el Cerebro

Nota nueva en 10-Conceptos/metropolis-hastings.md, enlazada a [[comparacion-de-modelos]], [[cadena-de-markov]] y [[bootstrap]]. Conviene rehacer a mano la Proposición 22.6: son dos casos y cuatro líneas, y de ella sale el Corolario 22.7, que es la razón entera de que el método exista.

¿Qué problema resuelve MCMC?::Muestrear del posterior cuando la constante de normalización no se puede calcular
¿Qué es el balance detallado?::π(x)P(x,y) = π(y)P(y,x): en equilibrio pasa tanta masa en un sentido como en el otro
¿Qué garantiza el balance detallado?::Que π es estacionaria; es condición suficiente y local, no necesaria
¿Cuál es la razón de aceptación de Metropolis-Hastings?::min(1, π̃(y)q(x|y) / π̃(x)q(y|x))
¿Por qué no hace falta la constante Z?::Porque π̃ solo aparece como cociente π̃(y)/π̃(x), y la constante se cancela
¿Qué se hace al rechazar una propuesta?::Quedarse en el estado actual y anotarlo otra vez; descartarlo rompe la estacionariedad
¿Qué pasa con propuesta simétrica?::La razón se reduce a min(1, π̃(y)/π̃(x)): se sube siempre y se baja con reticencia
¿Qué no garantiza la demostración?::Que la cadena converja desde cualquier arranque, ni cuándo, ni que las muestras sean independientes
¿Cuánto vale la varianza de la media de una cadena?::σ²τ/n con τ = 1 + 2Σρₖ, en vez de σ²/n
¿Qué es el tamaño de muestra efectivo?::n/τ: cuántas muestras independientes darían la misma precisión
¿Qué tasa de aceptación conviene en una dimensión?::Alrededor de 0,44; en el ejemplo el mejor paso la dio en 0,44
¿Cuáles son los dos desastres del paso?::Muy corto acepta casi todo y avanza poco; muy largo se rechaza casi todo y la cadena se queda quieta

Fuentes

Lo que esta página demuestra sola. Las Proposiciones 22.4, 22.6, 22.8 y 22.9 y el Corolario 22.7 se demuestran aquí, a partir de las definiciones. Se comprueban además numéricamente: sobre un espacio de seis estados se construye la matriz de transición entera y se verifica el balance detallado y la estacionariedad por debajo de 10^{-12}, que el autovector de autovalor uno es \pi, y que todas las filas de P^{400} coinciden con \pi a 10^{-10}; en el caso continuo, la cadena reproduce la media y la desviación de una \text{Beta}(9,5) dentro de 0{,}005; y la Proposición 22.9 se usa para medir \tau en cinco pasos de propuesta. El visual se contrastó antes de publicarse contra la densidad exacta en todo el rango de sus controles.

Lo que se usa de otras lecciones sin repetir. La conjugación Beta–binomial, que aquí solo sirve de patrón de comparación, es la lección 19. La integral que se evita es la de la lección 21. La idea de responder por simulación lo que no sale en forma cerrada es la lección 11. La autocorrelación es la lección 5.

Lo que se enuncia sin demostrar. Que irreducibilidad y aperiodicidad garantizan la convergencia desde cualquier arranque se enuncia y no se demuestra: es el teorema ergódico para cadenas de Markov, que excede el alcance de la página y queda señalado como tal en la sección 4. La convergencia dominada del límite de la Proposición 22.9 se invoca sin justificar las hipótesis. El valor 0{,}44 como tasa de aceptación óptima en una dimensión se cita como resultado conocido y se comprueba numéricamente, pero no se deriva.

Lo que viene de los libros.

• Think Bayes 2e (Downey, CC BY-NC-SA 4.0), §19 MCMC — citada textualmente una frase. El capítulo cubre el mismo terreno usando PyMC y sin escribir el algoritmo; publica sus secciones sin numerar, por lo que se cita a nivel de capítulo.

Lo que es mío, no del libro. La derivación entera desde el balance detallado, y la Proposición 22.6 demostrada por los dos casos de la razón. La comprobación con la matriz de transición completa sobre seis estados, que permite verificar la teoría sin simular nada. La Proposición 22.9 enunciada con el término (1-k/n) explícito y demostrada por conteo de retardos. Y la tabla de la sección 4, construida para que la tasa óptima aparezca sola en lugar de citarla.

Índice de Think Bayes 2e verificado el 12-09-2026: 193 secciones de los capítulos 1 a 18 y 20; el 19 publica sus secciones sin numerar.