Diagnóstico de cadenas

Estadística · Lección 24

L4 Derivación Fase F5 Libro autocontenida Prereqs Gibbs sampling · Metropolis-Hastings

Objetivo

Al terminar esta lección se puede:

  1. Definir \hat{R} a partir de la varianza dentro y entre cadenas, y demostrar la identidad exacta que lo expresa como cociente de esas dos dispersiones.
  2. Demostrar que con cadenas ya estacionarias \hat{R}^{2}\approx1+(\tau-1)/n, y convertir con eso el umbral habitual \hat{R}<1{,}01 en una afirmación sobre el número de muestras efectivas por cadena.
  3. Calcular el error de Monte Carlo de cualquier resumen y decidir con él cuántas cifras se pueden escribir.
  4. Construir un caso donde \hat{R} da el visto bueno y la respuesta está mal, y explicar qué clase de fallo ningún diagnóstico de este tipo puede ver.

De dónde viene

  • Estadística 22: el factor de ineficiencia \tau y el tamaño efectivo. Aquí se conectan con \hat{R} por una fórmula, de modo que los dos diagnósticos dejan de ser independientes.
  • Estadística 23: la cadena de Gibbs sobre parámetros correlacionados, que es el banco de pruebas de toda la lección porque su \tau se conoce exactamente.
  • Estadística 4: la descomposición de una varianza en la parte de dentro de los grupos y la parte de entre grupos. \hat{R} es esa descomposición con las cadenas como grupos.
  • Estadística 7: el error estándar de una media. El error de Monte Carlo es el mismo objeto con n_{\text{ef}} en lugar de n.

Para qué sirve después

  • Lección 25: PyMC reporta \hat{R} y el tamaño efectivo en cada ajuste; esta lección es la que permite leerlos en vez de obedecerlos.
  • Estadística 16: un posterior mal muestreado da pronósticos mal calibrados, y el defecto no se distingue de un modelo malo si no se diagnostica antes la cadena.
  • Estadística 21: estimar una verosimilitud marginal desde una cadena exige mucha más precisión que estimar una media, y el error de Monte Carlo de la sección 4 dice cuánta.

Notación

Símbolo Se lee Significado
m eme Número de cadenas independientes
n ene Longitud de cada cadena tras descartar el arranque
\bar\theta_j theta barra sub jota La media de la cadena j
W doble ve La varianza media dentro de las cadenas
S_\mu^{2} ese dos sub mu La varianza muestral de las m medias de cadena
B be La varianza entre cadenas, n\,S_\mu^{2}
\hat{V} ve gorro El estimador combinado de la varianza del objetivo
\hat{R} erre gorro \sqrt{\hat{V}/W}, el diagnóstico de Gelman y Rubin
\tau tau El factor de ineficiencia de la Proposición 22.9

1. Dos preguntas distintas

Una cadena puede fallar de dos maneras y conviene no mezclarlas.

La primera es no haber llegado: la cadena sigue arrastrando el punto donde arrancó y todavía no explora el objetivo. Eso se detecta comparando cadenas que arrancaron en sitios distintos, y es lo que mide \hat{R}.

La segunda es haber llegado pero con poca información: la cadena está en el objetivo pero sus muestras se parecen tanto entre sí que mil iteraciones valen por diez. Eso se mide con el \tau de la lección 22 y no necesita más de una cadena.

Las dos preguntas no son independientes, y la sección 3 lo demuestra: con cadenas ya estacionarias, \hat{R} queda determinado por \tau y por n.

2. La definición y lo que mide

Definición 24.1 (\hat{R} de Gelman-Rubin). Con m cadenas de longitud n, sean W=\frac{1}{m}\sum_{j=1}^{m}s_j^{2},\qquad S_\mu^{2}=\frac{1}{m-1}\sum_{j=1}^{m}(\bar\theta_j-\bar\theta)^{2},\qquad B=n\,S_\mu^{2}, con s_j^2 la varianza muestral de la cadena j. Se define \hat{V}=\frac{n-1}{n}\,W+\frac{B}{n},\qquad \hat{R}=\sqrt{\frac{\hat{V}}{W}}.

La idea de fondo es la de un análisis de varianza con las cadenas como grupos. W mide cuánto se mueve cada cadena por dentro; S_\mu^2, cuánto discrepan entre sí. Si todas exploran el mismo objetivo, lo segundo tiene que ser pequeño frente a lo primero.

Proposición 24.2 (identidad exacta). Para todo m,n\ge2, \hat{R}^{2}=1-\frac{1}{n}+\frac{S_\mu^{2}}{W}.

Demostración. Sustituyendo B=n\,S_\mu^2 en la definición de \hat{V}, \hat{V}=\frac{n-1}{n}W+\frac{n\,S_\mu^{2}}{n}=\Big(1-\frac1n\Big)W+S_\mu^{2}. Dividiendo entre W queda el enunciado. ∎

La identidad quita todo el misterio: \hat{R}^2-1 es, salvo el término 1/n que se desvanece, la razón entre la dispersión de las medias de cadena y la dispersión dentro de las cadenas. Un \hat{R}=1{,}01 significa S_\mu^2/W\approx0{,}02: las medias de las cadenas discrepan entre sí un dos por ciento de lo que se mueve cada una por dentro.

3. Qué valor tiene \hat{R} cuando todo va bien

Que las cadenas hayan convergido no hace \hat{R} igual a uno: lo hace igual a algo que depende de \tau, y esa es la parte que suele faltar.

Proposición 24.3 (\hat{R} en el régimen estacionario). Supóngase que las m cadenas son independientes entre sí y estacionarias con el objetivo, de varianza \sigma^2 y factor de ineficiencia \tau en el sentido de la Proposición 22.9. Entonces, para n grande, E[S_\mu^{2}]=\frac{\sigma^{2}\tau}{n},\qquad E[W]\to\sigma^{2},\qquad\text{luego}\qquad E[\hat{R}^{2}]\;\approx\;1+\frac{\tau-1}{n}.

Demostración. Cada \bar\theta_j es la media de una cadena estacionaria, así que por la Proposición 22.9 su varianza es \sigma^2\tau/n. Las m medias son independientes entre sí y tienen todas la misma esperanza, luego su varianza muestral S_\mu^2 es insesgada para esa varianza: E[S_\mu^2]=\sigma^2\tau/n. Por otro lado s_j^2 estima la varianza marginal de la cadena, que es \sigma^2, con sesgo que tiende a cero; luego E[W]\to\sigma^2. Sustituyendo las dos en la Proposición 24.2, E[\hat{R}^{2}]\approx1-\frac1n+\frac{\sigma^{2}\tau/n}{\sigma^{2}}=1+\frac{\tau-1}{n}. ∎

Corolario 24.4 (qué dice el umbral). Si se exige \hat{R}\le1+\epsilon con \epsilon pequeño, entonces aproximadamente \frac{n}{\tau}\;\gtrsim\;\frac{1}{2\epsilon}. Con el umbral habitual \epsilon=0{,}01 eso son unas cincuenta muestras efectivas por cadena.

Demostración. De la Proposición 24.3, \hat{R}\le1+\epsilon equivale a 1+(\tau-1)/n\le(1+\epsilon)^2=1+2\epsilon+\epsilon^2. Despreciando \epsilon^2 y aproximando \tau-1 por \tau, queda \tau/n\le2\epsilon, es decir n/\tau\ge1/(2\epsilon). Con \epsilon=0{,}01, n/\tau\ge50. ∎

El corolario explica por qué \hat{R} y el tamaño efectivo dicen casi lo mismo cuando las cadenas ya convergieron, y por qué aun así hay que mirar los dos: cincuenta muestras efectivas por cadena es un mínimo muy bajo para reportar un intervalo, y \hat{R} se queda tranquilo mucho antes de que la precisión sea aceptable.

Las tres filas siguen la fórmula de cerca. Y la última columna dice lo que importa: con n=8000 y \tau\approx100 hay ochenta muestras efectivas por cadena, que es poco, y \hat{R} ya vale 1{,}007.

El visual siguiente deja ver las dos cosas a la vez. El marco está anclado a la población: los ejes no dependen de las cadenas que se dibujan.

4. Cuántas cifras se pueden escribir

Definición 24.5 (error de Monte Carlo). Para un resumen que sea una media sobre la cadena, con desviación posterior \sigma y tamaño efectivo total n_{\text{ef}}, el error de Monte Carlo es \text{EMC}=\frac{\sigma}{\sqrt{n_{\text{ef}}}}.

Proposición 24.6 (cifras defendibles). Escribir k decimales de una media posterior exige \text{EMC}\lesssim\tfrac12\cdot10^{-k}, es decir n_{\text{ef}}\;\gtrsim\;4\,\sigma^{2}\,10^{2k}.

Demostración. El último decimal escrito es defendible si el error de muestreo no llega a cambiarlo, lo que pide \text{EMC}\le\frac12 10^{-k}. Elevando al cuadrado la Definición 24.5, \sigma^2/n_{\text{ef}}\le\frac14 10^{-2k}, y despejando queda el enunciado. ∎

La consecuencia práctica es incómoda y conviene tenerla presente: con \sigma=0{,}124, escribir tres decimales de la media posterior exige unas 61\,200 muestras efectivas, que con \tau=100 son más de seis millones de iteraciones. La mayoría de los resultados publicados con MCMC llevan más cifras de las que su cadena sostiene.

5. Lo que ningún \hat{R} puede ver

Proposición 24.7 (el punto ciego). Supóngase que el objetivo tiene dos modas separadas por una región de probabilidad despreciable, y que las m cadenas arrancan todas dentro de la misma moda y no la abandonan en n pasos. Entonces cada cadena es aproximadamente estacionaria respecto del objetivo restringido a esa moda, de modo que \hat{R}\to1 y todos los diagnósticos basados en comparar cadenas dan el visto bueno, mientras el estimador de cualquier esperanza converge al valor de la moda visitada y no al del objetivo.

Demostración. Por hipótesis ninguna cadena cruza, así que las m trayectorias son realizaciones de la misma cadena restringida al mismo subconjunto A. La Proposición 24.3 aplicada a ese objetivo restringido da E[\hat{R}^2]\approx1+(\tau_A-1)/n\to1. Por otro lado, para cualquier f, \frac1{mn}\sum_{j,i}f(\theta_{ji})\;\longrightarrow\;E[f\mid A]\neq E[f] salvo coincidencia, porque el objetivo pone masa P(A)<1 en A. El diagnóstico compara cadenas entre sí y todas coinciden; nada en la comparación mira fuera de A. ∎

Ejemplo 24.8. Objetivo 0{,}7\,\mathcal{N}(4,1)+0{,}3\,\mathcal{N}(-4,1), de media verdadera 1{,}6. Cuatro cadenas de Metropolis con paso 0{,}5 arrancadas todas en +4: ninguna cruza en 20 000 pasos, \hat{R}=1{,}0003 y la media reportada es 4{,}00. El diagnóstico es perfecto y la respuesta está mal por 2{,}4. Repartiendo los arranques entre las dos modas, \hat{R} sube a 1{,}80 y la alarma suena, aunque la respuesta sigue siendo mala: lo que el diagnóstico detecta es que las cadenas discrepan, no que acierten.

El visual siguiente deja mover la separación entre las modas y el reparto de los arranques, para ver dónde empieza y dónde acaba el punto ciego.

De la Proposición 24.7 salen las dos únicas defensas que hay, y ninguna es un número. La primera es arrancar disperso: tantas cadenas como se pueda, desde puntos alejados entre sí, porque solo así la discrepancia tiene oportunidad de aparecer. La segunda es saber qué forma puede tener el posterior antes de muestrearlo: una mezcla, una verosimilitud con simetrías, un modelo no identificable. El diagnóstico confirma lo que se sospecha; no descubre lo que no se miró.

6. La lista

Antes de reportar nada salido de una cadena, cuatro cosas:

  1. \hat{R} por parámetro, no solo para el que interesa: un parámetro auxiliar mal muestreado contamina al resto. Umbral 1{,}01, con el Corolario 24.4 en mente sobre lo poco que ese umbral exige.
  2. Tamaño efectivo por parámetro, y compararlo con lo que la Proposición 24.6 pide para las cifras que se piensan escribir.
  3. Arranques dispersos y varias cadenas, por la Proposición 24.7. Cuatro es el mínimo habitual y no es generoso.
  4. La traza, mirada. Los diagnósticos resumen; una traza con escalones planos o con un salto en medio dice cosas que ningún resumen recoge.

Ejercicios

  1. Demostrar a partir de la Proposición 24.2 que \hat{R}\ge\sqrt{1-1/n} siempre, y que \hat{R}<1 es posible. ¿Qué significaría encontrarlo?
  2. Comprobar que \hat{V} de la Definición 24.1 es un promedio ponderado de W y de B, e identificar qué estima cada uno cuando las cadenas han convergido y cuando no.
  3. Usar el Corolario 24.4 para decir qué umbral de \hat{R} haría falta para garantizar mil muestras efectivas por cadena.
  4. En la Proposición 24.6, calcular cuántas muestras efectivas hacen falta para escribir tres decimales de un cuantil extremo, sabiendo que la varianza del estimador de un cuantil al 1 % es mucho mayor que la de una media. ¿Por qué los cuantiles de las colas son lo último que converge?
  5. Construir un objetivo con dos modas de igual peso y demostrar que, arrancando todas las cadenas en una de ellas, la media reportada es correcta por simetría aunque el muestreo esté roto. ¿Qué resumen sí delataría el fallo?

Reto

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

  1. Escribir diagnostico(cadenas) que devuelva \hat{R}, el tamaño efectivo y el error de Monte Carlo de la media, y verificarlo contra la Proposición 24.3 sobre la cadena de Gibbs con \tau conocido, para al menos cinco valores de r y tres longitudes.
  2. Reproducir el Ejemplo 24.8 y barrer la separación de las modas de 1 a 6. Encontrar la separación a partir de la cual \hat{R} deja de detectar el problema con arranques juntos, y relacionarla con la probabilidad de cruzar el valle en un paso.
  3. Implementar \hat{R} partido: cortar cada cadena por la mitad y tratar las mitades como cadenas distintas. Comprobar que eso detecta una cadena con deriva lenta que el \hat{R} ordinario deja pasar, y explicar por qué.

Del libro

Ninguno de los libros con índice verificado en verificar/indices.json trata el diagnóstico de cadenas. Think Bayes usa PyMC y menciona de pasada las cantidades que la librería reporta, sin definirlas; los demás no entran en MCMC. Por eso esta lección es autocontenida y no cita texto de nadie: la regla 3 pide preferir una página que se demuestra sola a una referencia que no se puede verificar.

Las cantidades tienen origen conocido —\hat{R} es de Gelman y Rubin, y la variante partida y las correcciones posteriores son de la literatura de Stan—, pero como no hay una fuente abierta y verificada en el índice, aquí se enuncian y se demuestran desde cero en lugar de atribuirse a un texto que no se ha abierto.

Para el Cerebro

Nota nueva en 10-Conceptos/diagnostico-de-cadenas.md, enlazada a [[metropolis-hastings]], [[gibbs-sampling]] y [[tamano-efectivo]]. Conviene rehacer a mano la Proposición 24.3: es la que convierte \hat{R} de regla de pulgar en una consecuencia de \tau, y de ella sale el Corolario 24.4, que dice lo poco que el umbral habitual exige.

¿Qué dos fallos distintos puede tener una cadena?::No haber llegado al objetivo, y haber llegado con poca información; R̂ mide el primero y τ el segundo
¿Qué compara R̂?::La dispersión entre las medias de varias cadenas contra la dispersión dentro de cada una
¿Cuál es la identidad exacta de R̂?::R̂² = 1 − 1/n + S²_μ/W, sin aproximaciones
¿Qué significa R̂ = 1,01?::Que las medias de las cadenas discrepan un 2 % de lo que se mueve cada cadena por dentro
¿Cuánto vale R̂ con cadenas ya estacionarias?::E[R̂²] ≈ 1 + (τ−1)/n: queda determinado por τ y por n, no vale uno
¿Qué exige realmente el umbral R̂ < 1,01?::Unas 50 muestras efectivas por cadena, que es muy poco para reportar un intervalo
¿Qué es el error de Monte Carlo?::σ/√n_ef: el error estándar de la 7 con el tamaño efectivo en lugar de n
¿Cuántas muestras hacen falta para k decimales?::n_ef ≳ 4σ²·10^(2k); con σ=0,124 y tres decimales, unas 61 200
¿Cuál es el punto ciego de R̂?::Si todas las cadenas se quedan en la misma moda y no cruzan, R̂ tiende a uno con la respuesta mal
¿Por qué no puede verlo?::Porque compara las cadenas entre sí, y todas coinciden; nada en la comparación mira fuera de la moda visitada
¿Cuáles son las dos defensas?::Arrancar disperso con varias cadenas, y saber de antemano qué forma puede tener el posterior
¿Qué mirar antes de reportar?::R̂ por parámetro, tamaño efectivo por parámetro, arranques dispersos, y la traza con los ojos

Fuentes

Lo que esta página demuestra sola. Todo. Las Proposiciones 24.2, 24.3, 24.6 y 24.7 y el Corolario 24.4 se demuestran aquí a partir de la Definición 24.1 y de la Proposición 22.9. Se comprueban además numéricamente: la Proposición 24.3 se contrasta sobre la cadena de Gibbs de la lección 23, cuyo \tau se conoce exactamente, en tres longitudes con 40 cadenas cada una y coincidencia en la tercera cifra; la Proposición 24.7 se realiza en el Ejemplo 24.8 con cuatro cadenas de 20 000 pasos, cero cruces, \hat{R}=1{,}0003 y un error de 2{,}4 en la media. Los dos visuales se contrastaron antes de publicarse contra las fórmulas en todo el rango de sus controles.

Lo que se usa de otras lecciones sin repetir. El factor de ineficiencia \tau y el tamaño efectivo son la Proposición 22.9. La fórmula \tau=(1+r^2)/(1-r^2) de la cadena de Gibbs es la Proposición 23.6, y es lo que permite comparar contra un valor exacto en lugar de contra otra estimación. El error estándar de una media es la lección 7. La descomposición dentro/entre grupos es la lección 4.

Lo que se enuncia sin demostrar. Que E[W]\to\sigma^2 en la Proposición 24.3 se argumenta pero no se acota: el sesgo de la varianza muestral de una serie correlacionada es de orden \tau/n y se desprecia. La Proposición 24.7 supone que ninguna cadena cruza en n pasos, lo cual es una hipótesis del enunciado y no un hecho demostrado. El \hat{R} partido del reto 3 se propone sin desarrollarlo.

Lo que viene de los libros. Nada. Ningún libro del índice verificado trata este material, y por la regla 3 no se cita lo que no se ha podido abrir. El origen de \hat{R} es conocido y se menciona en la sección Del libro sin atribuirle texto.

Lo que es mío, no del libro. La Proposición 24.2 escrita como identidad exacta, que quita el misterio a la definición. La Proposición 24.3, que ata \hat{R} a \tau, y el Corolario 24.4, que traduce el umbral 1{,}01 en «cincuenta muestras efectivas por cadena». La Proposición 24.6 sobre cifras defendibles. Y el Ejemplo 24.8, construido con pesos desiguales para que el fallo no se salve por simetría.

Lección autocontenida; ninguna cita de libro. Índices verificados el 12-09-2026.