Gibbs sampling desde cero

Estadística · Lección 23

L3 Derivación Fase F5 Libro Think Bayes 2e §19 Prereqs Metropolis-Hastings · Normal multivariante

Objetivo

Al terminar esta lección se puede:

  1. Escribir las condicionales completas de un modelo y con ellas el muestreador de Gibbs.
  2. Demostrar que Gibbs es Metropolis-Hastings con la condicional como propuesta, y que por eso su razón de aceptación vale exactamente uno.
  3. Demostrar que cada barrido deja invariante el objetivo, sin apoyarse en la lección 22.
  4. Demostrar que sobre una normal bivariante con correlación r la cadena de Gibbs es autorregresiva de coeficiente r^{2}, calcular su factor de ineficiencia (1+r^{2})/(1-r^{2}) y deducir de ahí cuándo Gibbs deja de servir y qué hacer.

De dónde viene

  • Estadística 22: el balance detallado y Metropolis-Hastings. Gibbs se obtiene eligiendo la propuesta de una manera concreta, y todo lo demostrado allí sigue valiendo.
  • Estadística 12: las condicionales de una normal multivariante. Son el ejemplo donde todo sale cerrado, y también el que revela el defecto del método.
  • Estadística 19: la conjugación. Gibbs funciona cuando cada bloque de parámetros es conjugado dado el resto, aunque el modelo entero no lo sea.
  • Estadística 5: la correlación. Aquí es a la vez lo que hace interesante el problema y lo que hace lento al método.

Para qué sirve después

  • Lección 24: la cadena de Gibbs sobre parámetros correlacionados es el caso de libro donde los diagnósticos tienen que dar la alarma.
  • Lección 25: PyMC usa Gibbs para los bloques donde la condicional es conocida y otro muestreador para el resto; saber cuál es cuál explica sus tiempos.
  • Estadística 15: la imputación múltiple de datos faltantes es un Gibbs donde uno de los bloques son los propios datos que faltan.
  • Estadística 13: el modelo de regresión con varianza desconocida es el ejemplo estándar de condicionales completas, y aquí se resuelve entero.

Notación

Símbolo Se lee Significado
\theta=(\theta_1,\dots,\theta_d) theta El vector de parámetros, partido en d bloques
\theta_{-j} theta menos jota Todos los bloques menos el j-ésimo
\pi(\theta_j\mid\theta_{-j},x) — La condicional completa del bloque j
r erre La correlación entre los dos parámetros de la sección 4
\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
\text{IG}(\alpha,\beta) i ge La inversa-gamma: si 1/V\sim\text{Gamma}(\alpha,\beta) entonces V\sim\text{IG}(\alpha,\beta)

1. Condicionales completas

Metropolis-Hastings solo pide evaluar el objetivo salvo constante, y eso lo hace universal. El precio es que hay que elegir una propuesta, y elegirla mal cuesta el factor \tau de la Proposición 22.9. Gibbs cambia ese trato: pide más al modelo y no pide nada al usuario.

Definición 23.1 (condicional completa). Para una partición \theta=(\theta_1,\dots,\theta_d), la condicional completa del bloque j es la distribución de \theta_j dados todos los demás y los datos, \pi(\theta_j\mid\theta_{-j},x)\;\propto\;\pi(\theta\mid x), donde la proporcionalidad es como función de \theta_j con \theta_{-j} fijo.

La proporcionalidad merece un comentario, porque es la que hace practicable el método: para obtener la condicional completa se toma el posterior conjunto y se tacha todo lo que no contenga \theta_j. Lo que queda suele ser el núcleo de una densidad conocida, por la Proposición 19.2. Un modelo puede no ser conjugado en bloque y serlo bloque a bloque, y eso es lo normal.

Definición 23.2 (muestreador de Gibbs). Desde el estado \theta^{(t)}, un barrido consiste en recorrer los bloques en orden y sortear cada uno de su condicional completa, usando siempre los valores más recientes de los demás: \theta_1^{(t+1)}\sim\pi(\cdot\mid\theta_2^{(t)},\dots,\theta_d^{(t)},x),\qquad \theta_2^{(t+1)}\sim\pi(\cdot\mid\theta_1^{(t+1)},\theta_3^{(t)},\dots,x),\ \dots

No hay paso de aceptación ni parámetro que ajustar. La razón es la proposición siguiente.

2. Por qué no hace falta aceptar nada

Proposición 23.3 (Gibbs es Metropolis-Hastings con aceptación uno). Sea la propuesta que cambia solo el bloque j sorteándolo de su condicional completa, q(y\mid x)=\pi(y_j\mid x_{-j},\,\text{datos}),\qquad y_{-j}=x_{-j}. Entonces la razón de aceptación de la Definición 22.5 vale \alpha(x,y)=1 para todo par.

Demostración. Como y_{-j}=x_{-j}, escríbase z=x_{-j}=y_{-j}. Por la definición de condicional, \pi(x)=\pi(x_j\mid z)\,\pi(z) y \pi(y)=\pi(y_j\mid z)\,\pi(z), con la misma marginal \pi(z) en los dos casos. Sustituyendo en la razón, \frac{\pi(y)\,q(x\mid y)}{\pi(x)\,q(y\mid x)} =\frac{\pi(y_j\mid z)\,\pi(z)\cdot\pi(x_j\mid z)}{\pi(x_j\mid z)\,\pi(z)\cdot\pi(y_j\mid z)}=1, porque q(x\mid y)=\pi(x_j\mid y_{-j})=\pi(x_j\mid z) y q(y\mid x)=\pi(y_j\mid z). Luego \alpha=\min(1,1)=1. ∎

La cancelación es total: la marginal \pi(z) se va porque el bloque z no se toca, y las dos condicionales se cruzan. Una propuesta que ya está distribuida como el objetivo condicional no puede mejorarse rechazándola.

Conviene también demostrarlo por la vía directa, sin invocar la lección 22, porque la demostración es otra y aclara qué se conserva en cada paso.

Proposición 23.4 (cada paso deja invariante el objetivo). Si \theta^{(t)}\sim\pi y se actualiza solo el bloque j sorteándolo de \pi(\cdot\mid\theta_{-j}^{(t)}), entonces \theta^{(t+1)}\sim\pi. En consecuencia un barrido completo, que es la composición de d pasos de ese tipo, también deja invariante \pi.

Demostración. Con z=\theta_{-j}, la densidad conjunta del nuevo estado es p(\theta_j^{(t+1)},z)=\underbrace{\pi(\theta_j^{(t+1)}\mid z)}_{\text{del sorteo}}\cdot\underbrace{\pi(z)}_{\text{marginal del estado viejo}}, donde se usó que z no cambia y que, por hipótesis, el estado viejo tenía marginal \pi(z). Ese producto es exactamente \pi(\theta_j^{(t+1)},z) por definición de condicional. Para el barrido completo basta componer: cada paso preserva \pi, luego la composición también. ∎

Obsérvese que la Proposición 23.4 no afirma el balance detallado, y de hecho un barrido completo en orden fijo no lo cumple: recorrer los bloques siempre en el mismo sentido introduce una dirección. Es el caso, anunciado en la lección 22, de una cadena con distribución estacionaria correcta y sin reversibilidad. La estacionariedad se demuestra directamente y basta.

3. El modelo normal con media y varianza desconocidas

Este es el caso que justifica el método, porque el posterior conjunto no es de ninguna familia con nombre y las dos condicionales sí.

Con datos x_1,\dots,x_n normales de media \mu y varianza \sigma^2, ambas desconocidas, y el prior \pi(\mu,\sigma^2)\propto1/\sigma^2, el posterior conjunto es proporcional a

(\sigma^{2})^{-n/2-1}\exp\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_i-\mu)^{2}\right).

Tachando lo que no depende de cada parámetro se leen las dos condicionales:

Proposición 23.5 (condicionales completas del modelo normal). Con el posterior anterior, \mu\mid\sigma^{2},x\;\sim\;\mathcal{N}\!\left(\bar{x},\;\frac{\sigma^{2}}{n}\right), \qquad \sigma^{2}\mid\mu,x\;\sim\;\text{IG}\!\left(\frac{n}{2},\;\frac{1}{2}\sum_{i=1}^{n}(x_i-\mu)^{2}\right).

Demostración. Para \mu, el único factor que la contiene es la exponencial. Usando la identidad \sum_i(x_i-\mu)^2=\sum_i(x_i-\bar{x})^2+n(\bar{x}-\mu)^2 y descartando el primer sumando, que no depende de \mu, \pi(\mu\mid\sigma^{2},x)\;\propto\;\exp\left(-\frac{n(\mu-\bar{x})^{2}}{2\sigma^{2}}\right), núcleo de una normal de media \bar{x} y varianza \sigma^{2}/n. Para \sigma^2, con S(\mu)=\sum_i(x_i-\mu)^2 fijo, \pi(\sigma^{2}\mid\mu,x)\;\propto\;(\sigma^{2})^{-n/2-1}\exp\!\big(-\tfrac{S(\mu)}{2\sigma^{2}}\big), que es el núcleo de una \text{IG}(n/2,\,S(\mu)/2). En ambos casos la Proposición 19.2 identifica la densidad. ∎

Ninguna de las dos actualizaciones es el posterior de nadie: son condicionales. El posterior marginal de \mu sí se conoce en este caso —es una t de Student, como en la lección 9— y por eso sirve de comprobación exacta.

Nunca se calculó el posterior conjunto, ni su constante, ni la t de Student. Solo se alternaron dos sorteos elementales, y el resultado reproduce las marginales exactas.

4. Dónde se rompe

El defecto de Gibbs se ve entero en el ejemplo más simple posible, y se puede demostrar.

Proposición 23.6 (Gibbs sobre una normal bivariante es autorregresivo). Sea el objetivo \mathcal{N}(0,\Sigma) con \Sigma=\begin{pmatrix}1&r\\r&1\end{pmatrix} y |r|<1. La cadena de Gibbs que actualiza primero x y luego y cumple x^{(t+1)}=r^{2}\,x^{(t)}+\varepsilon_t, con \varepsilon_t independiente de x^{(t)}. En consecuencia \rho_k=r^{2k} y \tau=1+2\sum_{k\ge1}r^{2k}=\frac{1+r^{2}}{1-r^{2}}.

Demostración. Por la lección 12, las condicionales de la normal bivariante estandarizada son x\mid y\sim\mathcal{N}(ry,\,1-r^{2}) y y\mid x\sim\mathcal{N}(rx,\,1-r^{2}). El barrido es entonces x^{(t+1)}=r\,y^{(t)}+\sqrt{1-r^{2}}\,\xi_t,\qquad y^{(t+1)}=r\,x^{(t+1)}+\sqrt{1-r^{2}}\,\eta_t, con \xi_t,\eta_t normales estándar independientes. Sustituyendo la segunda igualdad, retrasada un paso, en la primera: x^{(t+1)}=r\big(r\,x^{(t)}+\sqrt{1-r^{2}}\,\eta_{t-1}\big)+\sqrt{1-r^{2}}\,\xi_t =r^{2}x^{(t)}+\varepsilon_t, con \varepsilon_t=r\sqrt{1-r^{2}}\,\eta_{t-1}+\sqrt{1-r^{2}}\,\xi_t, que es independiente de x^{(t)} por serlo \xi_t y \eta_{t-1} de todo lo anterior. Un proceso autorregresivo de coeficiente \phi tiene \rho_k=\phi^{k}; aquí \phi=r^{2}. Finalmente, sumando la serie geométrica, \tau=1+2\sum_{k\ge1}r^{2k}=1+\frac{2r^{2}}{1-r^{2}}=\frac{1+r^{2}}{1-r^{2}}. ∎

La fórmula es el diagnóstico entero. Con r=0 vale \tau=1: las muestras son independientes y Gibbs es óptimo. Con r=0{,}9 vale 9{,}5; con r=0{,}99, casi cien; con r=0{,}999, mil. La causa es geométrica: Gibbs solo se mueve paralelo a los ejes, y un contorno muy alargado en diagonal deja un pasillo estrecho por el que hay que avanzar a pasitos ortogonales.

El visual siguiente enseña de dónde sale la fórmula. El marco está anclado a la población —los ejes van de -4 a 4 siempre, y la elipse dibujada es la de la distribución objetivo, no la de los puntos visitados—, así que subir la correlación estrecha el pasillo sin reescalar nada.

5. Qué hacer cuando se rompe

La Proposición 23.6 también dice la solución, porque el problema no está en el objetivo sino en los ejes. Si se rota el sistema de coordenadas a los ejes principales de \Sigma —los eigenvectores de la lección 9—, la correlación entre las nuevas coordenadas es cero y \tau vuelve a uno. Muestrear ahí y deshacer la rotación da muestras del objetivo original.

Eso es reparametrizar, y es la primera herramienta ante una cadena lenta. Importa dónde se aplica: rotar las columnas que una cadena lenta ya produjo no arregla nada, porque la lentitud está en el tiempo y no entre las columnas. Hay que muestrear en las coordenadas nuevas y deshacer la rotación al final. En un modelo real no se conoce \Sigma, pero se puede estimar de una cadena piloto, o quitar la correlación de raíz: centrar los predictores antes de una regresión elimina casi toda la correlación entre la ordenada y la pendiente, que es el caso de la lección 13.

La segunda herramienta es cambiar de muestreador. Gibbs está atado a los ejes por construcción; los métodos que usan el gradiente del logaritmo del objetivo se mueven en cualquier dirección y no sufren esto de la misma manera. Es lo que hace PyMC por defecto, y la lección 25 lo verá desde fuera.

Ejercicios

  1. Demostrar que si los bloques se recorren en orden aleatorio en cada barrido, la cadena resultante sí cumple el balance detallado. Explicar por qué el orden fijo no lo cumple y por qué da igual.
  2. En la Proposición 23.5, comprobar que la condicional de \sigma^2 usa S(\mu)=\sum_i(x_i-\mu)^2 y no \sum_i(x_i-\bar{x})^2, y decir qué se estaría suponiendo al usar la segunda.
  3. Deducir de la Proposición 23.6 el r a partir del cual hacen falta más de mil barridos para obtener diez muestras efectivas.
  4. Generalizar la Proposición 23.6 a tres bloques con matriz de covarianza equicorrelacionada, o explicar qué se complica al intentarlo.
  5. Demostrar que muestrear (\mu,\sigma^2) con Gibbs y quedarse solo con la columna de \mu da muestras de la marginal posterior de \mu, sin tener que integrar \sigma^2 a mano.

Reto

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

  1. Escribir gibbs_normal(datos, N) que devuelva las cadenas de \mu y \sigma^2, y comparar el histograma de \mu con la densidad t exacta de la lección 9. Medir la distancia máxima entre las dos acumuladas y comprobar que baja como 1/\sqrt{N}.
  2. Implementar Gibbs para la regresión lineal con varianza desconocida: condicionales normal para los coeficientes e inversa-gamma para \sigma^2. Comprobar contra los resultados cerrados de la lección 14 y medir el \tau de cada coeficiente antes y después de centrar los predictores.
  3. Verificar la Proposición 23.6 en un barrido de cien valores de r entre 0 y 0,999: dibujar \tau medido contra la fórmula y localizar dónde la estimación de \tau empieza a fallar por cadenas demasiado cortas.

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 libro no escribe Gibbs: llama a PyMC, que elige el muestreador por su cuenta. Esta página lo deriva para que la elección de PyMC en la lección 25 se pueda leer, y sobre todo para que la Proposición 23.6 explique de antemano por qué unos modelos son rápidos y otros no.

Preguntas para leer §19 con lápiz:

  1. Los modelos del capítulo declaran una distribución por parámetro. ¿Cuáles de esas declaraciones dan condicionales completas conocidas, y cuáles obligarían a un paso de aceptación?
  2. El capítulo repite que hay que mirar las trazas. ¿Qué tendría el aspecto de una traza de la cadena de la Proposición 23.6 con r=0{,}99, y en qué se distinguiría de una cadena sana?
  3. PyMC informa del número efectivo de muestras. ¿Con qué cantidad de esta lección se corresponde, y qué relación guarda con el \tau que la Proposición 23.6 predice?

Para el Cerebro

Nota nueva en 10-Conceptos/gibbs-sampling.md, enlazada a [[metropolis-hastings]], [[condicional-completa]] y [[normal-multivariante]]. Conviene rehacer a mano la Proposición 23.6: sustituir una recursión en la otra son dos líneas, y de ahí sale la única fórmula que hace falta para saber si Gibbs va a servir.

¿Qué es una condicional completa?::La distribución de un bloque de parámetros dados todos los demás y los datos; se obtiene tachando del posterior conjunto lo que no contiene ese bloque
¿Qué es un barrido de Gibbs?::Recorrer los bloques en orden sorteando cada uno de su condicional completa, usando siempre los valores más recientes
¿Cuánto vale la aceptación en Gibbs?::Uno exactamente: es Metropolis-Hastings con la condicional como propuesta, y todo se cancela
¿Por qué se cancela?::La marginal del bloque que no se toca es la misma en ambos lados, y las dos condicionales se cruzan
¿Cumple Gibbs el balance detallado?::Con orden fijo no; la estacionariedad se demuestra directamente y basta
¿Cuándo sirve Gibbs?::Cuando cada bloque es conjugado dado el resto, aunque el modelo entero no lo sea
¿Cuáles son las condicionales del modelo normal?::mu | sigma² ~ N(x̄, sigma²/n) y sigma² | mu ~ IG(n/2, Σ(xᵢ−mu)²/2)
¿Qué estructura tiene la cadena sobre una normal bivariante?::Autorregresiva de coeficiente r²: x⁽ᵗ⁺¹⁾ = r²x⁽ᵗ⁾ + ruido
¿Cuánto vale tau ahí?::(1+r²)/(1−r²): uno si r=0, 9,5 si r=0,9 y casi cien si r=0,99
¿Por qué se vuelve lento?::Porque Gibbs solo se mueve paralelo a los ejes, y un contorno inclinado deja un pasillo estrecho
¿Cuál es la primera solución?::Reparametrizar: rotar a los ejes principales, o centrar los predictores en una regresión
¿Y la segunda?::Cambiar de muestreador: los que usan el gradiente no están atados a los ejes

Fuentes

Lo que esta página demuestra sola. Las Proposiciones 23.3, 23.4, 23.5 y 23.6 se demuestran aquí, a partir de la definición de condicional, de la razón de Metropolis-Hastings de la lección 22 y de las condicionales de la normal bivariante de la lección 12. Se comprueban además numéricamente: la razón de aceptación de Gibbs se evalúa en 20 000 pares de estados y se desvía de uno menos de 10^{-12}; el muestreador del modelo normal reproduce la media exacta de \mu dentro de 0{,}005 y su desviación y la media de \sigma^2 dentro del 2 %; y la Proposición 23.6 se contrasta en cuatro correlaciones, con \rho_1 dentro de 0{,}01 de r^2 en todas. La rotación de la sección 5 baja \tau de más de 50 a menos de 3. El visual se contrastó antes de publicarse contra la fórmula en todo el rango de sus controles.

Lo que se usa de otras lecciones sin repetir. La razón de aceptación y el factor \tau son la lección 22. Las condicionales de la normal multivariante son la lección 12. El reconocimiento de núcleos es la Proposición 19.2. La marginal t de Student del modelo normal es la lección 9. Los ejes principales son la lección 9 de Matemática.

Lo que se enuncia sin demostrar. Que un proceso autorregresivo de coeficiente \phi tiene \rho_k=\phi^k se usa en la Proposición 23.6 sin demostrarlo. Que la cadena converja a \pi desde cualquier arranque exige irreducibilidad, que aquí se cumple pero no se verifica. La afirmación de la sección 5 sobre muestreadores basados en el gradiente se enuncia como contexto para la lección 25.

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 terreno con PyMC y no escribe ningún muestreador; publica sus secciones sin numerar, por lo que se cita a nivel de capítulo.

Lo que es mío, no del libro. La Proposición 23.3 demostrada con la cancelación explícita, y la Proposición 23.4 demostrada por la vía directa para dejar claro que Gibbs con orden fijo no es reversible. La Proposición 23.6 entera —la estructura autorregresiva, \rho_k=r^{2k} y la fórmula de \tau—, que convierte «Gibbs va lento con parámetros correlacionados» en un número que se puede calcular antes de correr nada. Y la comprobación de la sección 5 con la rotación.

Í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.