7  Semana 7. HMC, diagnóstico y parametrización

SP-1653 Modelos Mixtos

7.1 Panorama de la semana

Durante las primeras seis semanas construimos modelos multinivel cada vez más expresivos y examinamos las consecuencias de sus supuestos. La semana anterior el foco estuvo en las distribuciones previas y en la simulación predictiva previa. A partir de esta semana aparece una pregunta adicional: aunque el modelo esté bien formulado, ¿tenemos evidencia de que el algoritmo computacional exploró adecuadamente la distribución posterior?

La pregunta orientadora es:

¿Cómo sabemos si el algoritmo exploró adecuadamente la distribución posterior y qué hacemos cuando los diagnósticos indican dificultades?

En modelos sencillos algunas distribuciones posteriores pueden obtenerse analíticamente. En los modelos mixtos y multinivel que hemos usado, en cambio, la posterior suele ser de dimensión alta y contiene dependencias fuertes entre coeficientes, desviaciones estándar y correlaciones. Por ello trabajamos con simulaciones aproximadas de la posterior. Stan utiliza Hamiltonian Monte Carlo (HMC) y, en particular, el algoritmo No-U-Turn Sampler (NUTS) para explorar espacios continuos de parámetros con ayuda del gradiente de la log densidad (McElreath 2020, chap. 9; Gelman et al. 2013, sec. 12.4).

El hecho de que un ajuste termine sin producir un error de software no implica que la aproximación sea confiable. Bayesian Workflow organiza el flujo computacional alrededor de dos etapas complementarias: obtener simulaciones de la posterior y diagnosticar activamente si esas simulaciones representan adecuadamente la distribución objetivo. Sus capítulos 11 y 12 enfatizan cadenas múltiples, \(\widehat R\), tamaños efectivos de muestra, error Monte Carlo, trazas y diagnósticos específicos de HMC, así como la necesidad de investigar la geometría de la posterior antes de responder mecánicamente a una advertencia (Gelman et al. 2026, chaps. 11-12).

En modelos multinivel esta discusión es especialmente importante. Una desviación estándar grupal cercana a cero puede producir regiones de curvatura muy distinta en la posterior. La parametrización centrada natural,

\[ \alpha_j\sim\mathcal N(\mu_\alpha,\tau_\alpha^2), \]

puede ser estadísticamente equivalente a una parametrización no centrada,

\[ z_j\sim\mathcal N(0,1), \qquad \alpha_j=\mu_\alpha+\tau_\alpha z_j, \]

pero las dos representaciones pueden ser muy diferentes desde el punto de vista de HMC. McElreath desarrolla esta conexión a través del funnel y las transiciones divergentes (McElreath 2020, 420-25). Matsuura dedica su capítulo 9 a separar tres fuentes de dificultades: parámetros no identificados, regularización insuficiente y dependencias posteriores que pueden aliviarse mediante reparametrización.

La secuencia de trabajo de esta semana será:

\[ \boxed{ \text{ajustar} \longrightarrow \text{diagnosticar} \longrightarrow \text{localizar la dificultad} \longrightarrow \text{modificar justificadamente} \longrightarrow \text{volver a ajustar} } \]

La palabra justificadamente es crucial. Aumentar adapt_delta, ampliar max_treedepth o ejecutar más iteraciones puede ser útil en ciertas situaciones, pero ninguna de esas acciones corrige automáticamente un modelo no identificado, una escala mal elegida o una estructura jerárquica con geometría desfavorable. Bayesian Workflow insiste en que una dificultad computacional puede revelar una dificultad de modelación y que modificar el modelo solo hasta conseguir que desaparezcan las advertencias puede alterar la inferencia de manera poco transparente (Gelman et al. 2026, secs. 12.3-12.5).

Objetivos de aprendizaje

Al finalizar esta semana, se espera que la persona estudiante pueda:

1. explicar por qué MCMC permite aproximar esperanzas y otras cantidades posteriores;
2. distinguir el número bruto de simulaciones del número efectivo de simulaciones;
3. interpretar autocorrelación, tamaño efectivo de muestra y error estándar Monte Carlo;
4. explicar por qué se utilizan varias cadenas con valores iniciales diferentes;
5. interpretar \(\widehat R\) como un diagnóstico de mezcla entre y dentro de cadenas, sin tratarlo como una prueba definitiva de convergencia;
6. utilizar gráficos de trazas y de rangos para localizar patrones de mezcla deficiente;
7. describir conceptualmente la diferencia entre Metropolis-Hastings, Gibbs, HMC y NUTS;
8. explicar la función de la energía potencial, la energía cinética y el gradiente en HMC;
9. describir el integrador leapfrog y relacionar su tamaño de paso con el error numérico;
10. explicar por qué NUTS adapta la longitud de las trayectorias en lugar de fijarla manualmente;
11. distinguir la fase de warmup de la fase de muestreo posterior;
12. explicar qué parámetros del algoritmo se adaptan durante warmup y por qué el escalamiento de los parámetros importa;
13. interpretar una transición divergente como evidencia de una región de la posterior difícil de integrar numéricamente;
14. distinguir una divergencia de una saturación de la profundidad máxima del árbol;
15. explicar el efecto de aumentar adapt_delta y por qué no debe usarse como sustituto de revisar el modelo;
16. explicar por qué aumentar max_treedepth puede incrementar el costo sin corregir la geometría que originó el problema;
17. reconocer la geometría tipo funnel que puede aparecer en modelos jerárquicos;
18. demostrar la equivalencia probabilística entre parametrizaciones centradas y no centradas de un coeficiente grupal normal;
19. explicar por qué la parametrización no centrada suele ser útil cuando cada grupo aporta poca información sobre su coeficiente;
20. reconocer situaciones en las que la parametrización centrada puede ser más eficiente;
21. distinguir falta de identificación estadística de dificultad computacional;
22. reconocer que una previa propia puede regularizar un parámetro débilmente informado sin convertir información previa en información proveniente de los datos;
23. utilizar posterior, bayesplot, brms y cmdstanr para resumir \(\widehat R\), ESS, MCSE, trazas y diagnósticos de NUTS;
24. documentar un ajuste computacional de forma reproducible, incluyendo semillas, número de cadenas, iteraciones, warmup y controles de NUTS;
25. construir para el proyecto del curso un informe breve que separe diagnóstico computacional, identificación y comprobación predictiva.

7.2 Problema motivador: el mismo modelo, dos comportamientos computacionales

Considere nuevamente un modelo de interceptos variables:

\[ y_{ij}\mid\alpha_j,\sigma \sim \mathcal N(\alpha_j,\sigma^2), \]

\[ \alpha_j\mid\mu_\alpha,\tau_\alpha \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

Una forma directa de escribir la segunda línea es

\[ \alpha_j = \mu_\alpha+\varepsilon_j, \qquad \varepsilon_j\sim\mathcal N(0,\tau_\alpha^2). \]

Equivalentemente, podemos introducir una variable estándar \(z_j\):

\[ z_j\sim\mathcal N(0,1), \]

\[ \alpha_j = \mu_\alpha+\tau_\alpha z_j. \tag{7.1}\]

Las dos formulaciones inducen exactamente la misma distribución normal para \(\alpha_j\) condicional en \(\mu_\alpha\) y \(\tau_\alpha\). Por tanto, desde el punto de vista del modelo probabilístico sobre \(\alpha_j\), no hemos cambiado la afirmación sustantiva.

Sin embargo, las coordenadas usadas durante el muestreo sí cambiaron:

• en la parametrización centrada, el algoritmo explora directamente \(\alpha_j\) junto con \(\mu_\alpha\) y \(\tau_\alpha\);
• en la parametrización no centrada, el algoritmo explora \(z_j\), \(\mu_\alpha\) y \(\tau_\alpha\), y obtiene \(\alpha_j\) mediante Ecuación 7.1.

Cuando \(\tau_\alpha\) puede ser muy pequeña y los datos de cada grupo aportan poca información, la parametrización centrada puede producir una posterior cuya escala cambia bruscamente con \(\tau_\alpha\). Esa geometría es un caso del llamado funnel. HMC puede encontrar dificultades para usar un único tamaño de paso que funcione bien tanto en la parte ancha como en la parte estrecha de la distribución (McElreath 2020, 420-25; Gelman et al. 2026, sec. 12.3).

La situación plantea una idea central de esta semana:

\[ \boxed{ \text{equivalencia probabilística} \not\Rightarrow \text{equivalencia computacional} } \]

ImportanteLos diagnósticos describen la aproximación computacional

Un buen \(\widehat R\), un ESS grande y ausencia de divergencias son evidencia favorable sobre el muestreo de la posterior definida por el modelo. No demuestran que el modelo sea sustantivamente adecuado ni que sus predicciones reproduzcan los aspectos relevantes de los datos.

Por eso distinguiremos siempre:

1. diagnóstico computacional;
2. identificación y regularización;
3. comprobación predictiva del modelo.

7.3 De la posterior a una aproximación Monte Carlo

Sea \(\theta\) el vector de parámetros y sea

\[ p(\theta\mid y) \]

la distribución posterior. Muchas cantidades de interés pueden escribirse como esperanzas posteriores. Para una función \(h(\theta)\),

\[ \operatorname E[h(\theta)\mid y] = \int h(\theta)p(\theta\mid y)\,d\theta. \]

Si pudiéramos obtener simulaciones independientes

\[ \theta^{(1)},\ldots,\theta^{(S)} \sim p(\theta\mid y), \]

aproximaríamos la esperanza mediante

\[ \widehat{\operatorname E}[h(\theta)\mid y] = \frac{1}{S}\sum_{s=1}^{S}h(\theta^{(s)}). \tag{7.2}\]

El mismo principio permite aproximar cuantiles, probabilidades posteriores y distribuciones predictivas. El problema es que, en modelos complejos, obtener simulaciones independientes de la posterior no suele ser sencillo. MCMC construye una secuencia dependiente que, bajo condiciones apropiadas, visita el espacio de parámetros de acuerdo con la distribución objetivo (Gelman et al. 2013, chaps. 11-12).

7.3.1 Dependencia entre simulaciones

Una cadena de Markov produce valores

\[ \theta^{(1)},\theta^{(2)},\ldots \]

para los cuales el siguiente estado depende del estado actual. Por tanto, dos iteraciones consecutivas pueden contener información redundante. Una cadena con 4000 simulaciones guardadas no necesariamente contiene la misma información Monte Carlo que 4000 simulaciones independientes.

Para una cantidad escalar \(h(\theta)\), una aproximación pedagógica del tamaño efectivo es

\[ S_{\text{eff}} \approx \frac{S} {1+2\sum_{t=1}^{\infty}\rho_t}, \tag{7.3}\]

donde \(\rho_t\) es la autocorrelación en el rezago \(t\). Esta expresión ayuda a entender la idea: autocorrelaciones positivas persistentes reducen la información efectiva. Las implementaciones modernas utilizadas por posterior y Stan emplean estimadores más robustos, incluyendo transformaciones por rangos y versiones diferenciadas para el centro y las colas de la distribución (Gelman et al. 2026, sec. 11.5).

7.3.2 Error estándar Monte Carlo

La posterior tiene incertidumbre porque desconocemos los parámetros. La aproximación numérica añade una segunda fuente: incertidumbre Monte Carlo. Para la media posterior de una cantidad escalar, una relación útil es

\[ \operatorname{MCSE}(\bar h) \approx \frac{\operatorname{sd}(h\mid y)} {\sqrt{S_{\text{eff}}}}. \tag{7.4}\]

La MCSE responde una pregunta computacional:

Si repitiéramos el algoritmo con otra semilla y el mismo modelo, ¿cuánto podría cambiar el resumen numérico únicamente por usar una simulación finita?

No debe confundirse con la desviación estándar posterior. La desviación posterior describe incertidumbre sobre el parámetro; la MCSE describe incertidumbre en nuestra aproximación numérica de un resumen de esa posterior.

NotaMás iteraciones no equivalen a más información sustantiva

Aumentar \(S\) reduce el error Monte Carlo cuando el algoritmo está explorando la distribución correcta. No reduce la incertidumbre posterior debida a pocos datos ni resuelve una falta de identificación.

7.4 Cadenas, valores iniciales y mezcla

Una única trayectoria puede dar una impresión engañosa de estabilidad. Por ello ejecutamos varias cadenas con valores iniciales distintos. Si todas exploran la misma distribución estacionaria, sus distribuciones marginales deberían ser semejantes después de la fase de adaptación.

Bayesian Workflow recomienda utilizar al menos cuatro cadenas independientes como práctica general, precisamente porque múltiples cadenas aumentan la posibilidad de detectar mala adaptación, mezcla deficiente o multimodalidad (Gelman et al. 2026, sec. 11.4).

7.4.1 ¿Qué significa mezclar?

Informalmente, una cadena mezcla bien cuando recorre con suficiente rapidez las regiones relevantes de la posterior, de manera que:

• no queda atrapada durante largos períodos en una región estrecha;
• cadenas iniciadas en puntos diferentes visitan regiones comparables;
• la dependencia serial no destruye la eficiencia de la simulación;
• las cantidades de interés se estiman con error Monte Carlo adecuado al propósito del análisis.

La mezcla es una propiedad de la interacción entre la distribución objetivo, la parametrización y el algoritmo. No existe un único número que la resuma completamente.

7.4.2 Gráficos de trazas

Un gráfico de trazas representa el valor de una cantidad frente a la iteración, separando las cadenas. Puede revelar:

• cadenas ubicadas en regiones distintas;
• desplazamientos lentos;
• períodos largos con poca exploración;
• multimodalidad visible;
• diferencias sistemáticas entre cadenas.

Cuando las cadenas mezclan bien, las trazas suelen superponerse sin patrones persistentes. Cuando existe un problema, las trazas ayudan a formular una hipótesis sobre su causa. Bayesian Workflow señala que estos gráficos son especialmente útiles para investigar problemas, no necesariamente para decorar un informe cuando todos los diagnósticos son favorables (Gelman et al. 2026, sec. 11.4).

7.4.3 Gráficos basados en rangos

Los gráficos de rangos comparan cómo se distribuyen entre cadenas los rangos de las simulaciones combinadas. Si las cadenas exploran la misma distribución, ninguna debería concentrarse sistemáticamente en rangos bajos, medios o altos. Estos gráficos complementan las trazas porque pueden mostrar diferencias distributivas que son difíciles de detectar visualmente en una serie temporal.

7.5 \(\widehat R\): comparar variación entre y dentro de cadenas

El factor potencial de reducción de escala, \(\widehat R\), compara la variación entre cadenas con la variación dentro de las cadenas. La intuición clásica es:

\[ \boxed{ \text{si las cadenas exploran la misma distribución,} \quad \text{variación entre cadenas} \approx \text{variación dentro de cadenas.} } \]

Cuando las cadenas permanecen en regiones diferentes, la variación combinada tiende a exceder la variación observada dentro de cada cadena y \(\widehat R\) se aleja de 1.

Las implementaciones modernas no utilizan simplemente la versión histórica del diagnóstico. Stan y posterior emplean variantes divididas y normalizadas por rangos que son más sensibles a ciertos patrones de no estacionariedad y colas problemáticas. Para un ajuste final, Bayesian Workflow usa como práctica estándar \(\widehat R<1.01\) para todas las cantidades de interés, pero advierte que este umbral no constituye una garantía de que la posterior se haya explorado perfectamente (Gelman et al. 2026, sec. 11.4).

Advertencia\(\widehat R\approx1\) no certifica el ajuste

Es posible que varias cadenas queden atrapadas en la misma región y produzcan \(\widehat R\) aparentemente favorable. También puede existir una región de masa posterior que ninguna cadena encontró. Por ello \(\widehat R\) se interpreta junto con ESS, trazas, diagnósticos de HMC, experimentos con inicialización y conocimiento de la estructura del modelo.

7.6 Tamaño efectivo de muestra: centro y colas

En una distribución posterior no todas las cantidades se estiman con la misma eficiencia. Por eso los resúmenes modernos distinguen al menos:

• bulk-ESS: eficiencia para cantidades asociadas con la región central de la distribución;
• tail-ESS: eficiencia para cuantiles de cola.

Un modelo puede estimar razonablemente bien la media y, sin embargo, tener poca precisión Monte Carlo en un cuantil extremo. Esto importa cuando comunicamos intervalos posteriores, probabilidades de excedencia o predicciones de eventos raros.

Bayesian Workflow utiliza ESS menores que aproximadamente 100 como una señal de preocupación para evaluar convergencia; para el informe final de un proyecto, la cantidad necesaria puede ser mucho mayor y debe determinarse con la MCSE requerida para las cantidades de interés (Gelman et al. 2026, secs. 11.4-11.6, 12.5).

La regla importante no es “ESS debe superar un número mágico”, sino:

\[ \boxed{ \text{la precisión Monte Carlo debe ser suficiente para la decisión inferencial que queremos comunicar.} } \]

7.7 ¿Por qué HMC?

Antes de estudiar los diagnósticos específicos de Stan conviene ubicar HMC dentro de MCMC.

7.7.1 Metropolis-Hastings

Un algoritmo Metropolis-Hastings propone un nuevo estado a partir del estado actual y lo acepta con una probabilidad diseñada para preservar la distribución objetivo. Si la propuesta se mueve localmente sin información sobre la geometría, la cadena puede comportarse como una caminata aleatoria: pasos grandes se rechazan y pasos pequeños avanzan lentamente.

7.7.2 Gibbs

Gibbs actualiza bloques de parámetros utilizando distribuciones condicionales. Puede ser muy eficaz cuando esas condicionales son fáciles de muestrear y la dependencia entre bloques es moderada. Sin embargo, correlaciones posteriores fuertes pueden producir movimientos lentos porque una actualización condicional cambia solo una parte de la posición a la vez.

7.7.3 Hamiltonian Monte Carlo

HMC utiliza el gradiente de la log posterior para construir propuestas que pueden recorrer distancias grandes manteniéndose en regiones de alta probabilidad. Este es el contraste esencial con una caminata aleatoria. McElreath desarrolla esta intuición mediante trayectorias sobre una superficie y BDA presenta el algoritmo en términos de dinámica hamiltoniana (McElreath 2020, 270-79; Gelman et al. 2013, 300-305).

Bürkner señala que brms delega el muestreo a Stan y que HMC/NUTS produce simulaciones típicamente menos autocorrelacionadas que algoritmos de caminata aleatoria, a cambio de requerir gradientes y mayor trabajo por iteración (Bürkner 2017).

NotaGibbs y Metropolis en este curso

Los utilizamos como referencia conceptual e histórica. El flujo computacional del curso se centra en HMC/NUTS mediante Stan, cmdstanr y brms.

7.8 La mecánica básica de HMC

Sea

\[ \pi(\theta) = p(\theta\mid y) \]

la densidad objetivo, conocida posiblemente hasta una constante. Definimos una energía potencial

\[ U(\theta) = -\log\pi(\theta)+C. \]

HMC introduce una variable auxiliar de momento \(r\) y, comúnmente,

\[ r\sim\mathcal N(0,M), \]

donde \(M\) es una matriz de masa definida positiva. La energía cinética es

\[ K(r) = \frac{1}{2}r^{\mathsf T}M^{-1}r. \]

El hamiltoniano es

\[ H(\theta,r) = U(\theta)+K(r). \tag{7.5}\]

La distribución conjunta artificial sobre posición y momento es proporcional a

\[ \exp\{-H(\theta,r)\}. \]

Al marginalizar \(r\), recuperamos la distribución objetivo sobre \(\theta\). El momento se introduce únicamente para construir movimientos eficientes.

7.8.1 Ecuaciones de Hamilton

La dinámica ideal satisface

\[ \frac{d\theta}{dt} = \frac{\partial H}{\partial r} = M^{-1}r, \]

\[ \frac{dr}{dt} = -\frac{\partial H}{\partial\theta} = \nabla_\theta\log\pi(\theta). \]

El gradiente orienta la trayectoria según la geometría local de la posterior. En una dinámica exacta, \(H(\theta,r)\) permanece constante. En el computador necesitamos aproximar estas ecuaciones numéricamente.

7.9 El integrador leapfrog y el tamaño de paso

Una iteración elemental del integrador leapfrog con tamaño de paso \(\epsilon\) puede escribirse como

\[ r_{t+1/2} = r_t -\frac{\epsilon}{2}\nabla U(\theta_t), \]

\[ \theta_{t+1} = \theta_t +\epsilon M^{-1}r_{t+1/2}, \]

\[ r_{t+1} = r_{t+1/2} -\frac{\epsilon}{2}\nabla U(\theta_{t+1}). \tag{7.6}\]

Después de varios pasos obtenemos una propuesta distante del estado inicial. Como la integración es aproximada, la energía final puede diferir de la inicial. HMC corrige ese error con un paso de aceptación de Metropolis basado en

\[ \alpha = \min\left\{1, \exp\left[-H(\theta',r')+H(\theta,r)\right] \right\}. \]

Un \(\epsilon\) muy grande puede producir errores numéricos severos; un \(\epsilon\) demasiado pequeño produce trayectorias más costosas porque se necesitan muchos pasos para recorrer una distancia útil.

7.9.1 Escalas diferentes y matriz de masa

Si un parámetro varía en centésimas y otro en miles, una geometría isotrópica resulta ineficiente. La matriz de masa permite adaptar las escalas relevantes para el movimiento. Durante warmup, Stan estima escalas y, dependiendo de la configuración, puede adaptar información de covarianza para mejorar la eficiencia (Gelman et al. 2026, secs. 11.3, 12.3).

Esto conecta el diagnóstico con decisiones que ya aparecieron en semanas anteriores: centrar predictores, evitar escalas arbitrariamente distintas y especificar previas razonables puede reducir correlaciones posteriores y curvaturas innecesarias.

7.10 NUTS: no fijar manualmente la longitud de cada trayectoria

HMC básico requiere elegir dos cantidades importantes:

1. el tamaño de paso \(\epsilon\);
2. el número de pasos \(L\) por trayectoria.

Si \(L\) es pequeño, la propuesta apenas se aleja del estado actual. Si es demasiado grande, la trayectoria puede dar la vuelta y comenzar a regresar hacia regiones ya recorridas.

NUTS construye adaptativamente una trayectoria y utiliza un criterio de “U-turn” para detener la expansión cuando continuar deja de ser productivo. De esta manera, no necesitamos escoger manualmente un \(L\) único para todas las iteraciones. Stan utiliza NUTS como su estrategia HMC principal; McElreath presenta su función como una extensión que evita trayectorias que comienzan a regresar hacia el punto de partida (McElreath 2020, 274-76).

La longitud de la trayectoria sigue teniendo consecuencias computacionales. NUTS organiza los pasos en un árbol binario, de modo que una profundidad mayor permite un número de pasos leapfrog que crece aproximadamente de forma exponencial con la profundidad. Por eso alcanzar repetidamente la profundidad máxima puede ser muy costoso.

7.11 Warmup y adaptación

Las iteraciones iniciales de Stan cumplen funciones específicas de adaptación. Durante warmup, el algoritmo intenta alcanzar una región típica de la posterior y ajusta parámetros del muestreador, entre ellos el tamaño de paso y la métrica utilizada por HMC. Las iteraciones usadas para esta adaptación no forman parte de la muestra posterior que se utiliza para inferencia (Gelman et al. 2026, secs. 11.2-11.3).

McElreath resalta que este proceso no debe entenderse simplemente como el “burn-in” de algoritmos MCMC antiguos: el objetivo de warmup en Stan incluye explícitamente calibrar el algoritmo antes de guardar simulaciones (McElreath 2020, 274-75).

AdvertenciaReducir warmup no es una optimización inocua

Un warmup demasiado corto puede dejar al algoritmo con una adaptación pobre. Durante la exploración inicial de modelos puede tener sentido usar ajustes breves para detectar errores rápidamente; para resultados finales conviene utilizar una adaptación suficiente y verificar los diagnósticos (Gelman et al. 2026, secs. 11.3-11.4).

7.12 El conjunto típico y por qué el modo no es suficiente

En dimensión alta, la mayor parte de la masa de probabilidad no se concentra necesariamente cerca del punto de máxima densidad. Bayesian Workflow utiliza el concepto de conjunto típico para destacar que un algoritmo de muestreo debe recorrer una región de valores con densidades representativas de la masa posterior, no simplemente dirigirse hacia la moda (Gelman et al. 2026, sec. 11.1).

Esta idea ayuda a distinguir optimización de muestreo:

• un optimizador intenta localizar valores de alta densidad;
• HMC intenta recorrer la distribución de manera que las frecuencias de visita representen su masa.

En modelos jerárquicos, la forma del conjunto típico puede cambiar considerablemente bajo una reparametrización. Esa es una de las razones por las que dos representaciones matemáticamente equivalentes pueden tener eficiencias computacionales muy distintas.

7.13 Diagnósticos específicos de HMC/NUTS

\(\widehat R\), ESS y MCSE son diagnósticos generales de la simulación. HMC ofrece además señales internas relacionadas con la integración numérica y la construcción de trayectorias.

7.13.1 Transiciones divergentes

Una divergencia indica que la integración hamiltoniana acumuló un error de energía suficientemente grande como para que la trayectoria numérica deje de ser confiable. Esto tiende a ocurrir en regiones de gran curvatura o donde la escala local cambia rápidamente (McElreath 2020, 420-23).

Una divergencia no debe tratarse como una observación que podemos eliminar de la muestra. El problema es que la divergencia señala que la región donde ocurrió puede no estar siendo explorada adecuadamente. Por ello, incluso un número pequeño merece investigación cuando afecta una zona relevante de la posterior.

En modelos jerárquicos una causa frecuente es una geometría tipo funnel. También pueden aparecer por restricciones, escalas extremas o parametrizaciones con dependencias difíciles.

7.13.2 adapt_delta

adapt_delta controla el objetivo de aceptación usado durante la adaptación. Aumentarlo normalmente induce un tamaño de paso menor. Un paso más pequeño puede seguir con mayor precisión una región curva y reducir divergencias, pero también incrementa el número de evaluaciones del gradiente y, por tanto, el costo.

El flujo correcto es:

\[ \text{divergencias} \rightarrow \text{inspeccionar geometría} \rightarrow \begin{cases} \text{probar un paso más pequeño},\\ \text{reparametrizar},\\ \text{reescalar},\\ \text{regularizar justificadamente},\\ \text{simplificar si el modelo no está informado}. \end{cases} \]

McElreath muestra un ejemplo donde aumentar adapt_delta reduce divergencias pero la parametrización no centrada produce una mejora mucho mayor en eficiencia (McElreath 2020, 423-25). Bayesian Workflow advierte, además, que cuando una fracción grande de transiciones diverge, subir adapt_delta rara vez debe ser la primera y única respuesta (Gelman et al. 2026, sec. 12.3).

7.13.3 Profundidad del árbol

NUTS construye una trayectoria mediante expansiones binarias. El diagnóstico treedepth__ registra la profundidad utilizada. Si muchas iteraciones alcanzan max_treedepth, el algoritmo está solicitando trayectorias más largas que el límite permitido.

Esto suele indicar ineficiencia: correlaciones fuertes, escalas mal adaptadas o una geometría que obliga a recorrer trayectorias largas. Aumentar max_treedepth permite más trabajo por iteración, pero no elimina la causa de esa ineficiencia. Bayesian Workflow presenta ejemplos donde centrar un predictor reduce drásticamente la correlación posterior y el tiempo de cálculo, sin necesidad de ampliar el límite (Gelman et al. 2026, sec. 12.3).

ImportanteDivergencia y profundidad máxima no significan lo mismo

• Divergencia: posible fallo de integración en una región difícil; puede comprometer la exploración de la posterior.
• Profundidad máxima: la trayectoria se truncó por un límite de costo; suele señalar baja eficiencia y geometría difícil.

Ambas requieren investigación, pero no deben interpretarse como diagnósticos equivalentes.

7.13.4 Número de pasos leapfrog

El diagnóstico n_leapfrog__ ayuda a medir cuánto trabajo realizó NUTS en cada iteración. Dos ajustes con el mismo número de simulaciones pueden tener costos muy diferentes si uno necesita trayectorias mucho más largas.

Esto recuerda una idea importante: el número de simulaciones no es una medida suficiente de eficiencia. Conviene considerar ESS en relación con el costo computacional y con el número de evaluaciones del gradiente.

7.14 Geometría posterior: correlación y curvatura

Suponga una posterior bivariada aproximadamente normal con correlación cercana a 1. Aunque no haya un problema de identificación, una coordenada puede cambiar solo si la otra cambia casi simultáneamente. Una parametrización alineada con los ejes originales obliga al sampler a moverse a través de una región muy alargada.

En regresión, una fuente frecuente es un intercepto definido lejos de la región observada del predictor. Si \(x\) toma valores alrededor de 2000, el intercepto en \(x=0\) y la pendiente pueden quedar fuertemente correlacionados: una pequeña modificación en la pendiente exige una gran compensación en el intercepto. Centrar \(x\) alrededor de una referencia observada puede eliminar gran parte de esa correlación sin cambiar las predicciones del modelo.

Bayesian Workflow presenta precisamente un ejemplo de este tipo y muestra que la corrección conceptual —centrar el predictor— mejora simultáneamente interpretación y cómputo (Gelman et al. 2026, sec. 12.3).

7.14.1 Reescalamiento no es manipular el resultado

Si definimos

\[ x^*=\frac{x-c}{s}, \]

podemos transformar también el intercepto, la pendiente y sus previas para conservar el mismo modelo sobre la escala original. Esta reexpresión cambia las coordenadas, no la pregunta sustantiva. Una buena parametrización intenta que las escalas y dependencias sean manejables para el algoritmo sin alterar las implicaciones probabilísticas que queremos representar.

7.15 El funnel de Neal

Considere el siguiente sistema, sin datos observados:

\[ v\sim\mathcal N(0,3^2), \]

\[ x_j\mid v \sim \mathcal N\left(0,\exp(2v)\right), \qquad j=1,\ldots,J. \tag{7.7}\]

Equivalentemente, la desviación estándar condicional es \(\exp(v)\). Cuando \(v\) es grande, los valores \(x_j\) pueden variar en un intervalo amplio. Cuando \(v\) es muy negativo, todos los \(x_j\) deben quedar en una región muy estrecha alrededor de cero. La densidad conjunta contiene entonces una región ancha conectada con un cuello muy estrecho.

McElreath utiliza una versión bidimensional para mostrar por qué HMC puede generar divergencias en la región estrecha y por qué la transformación

\[ z_j\sim\mathcal N(0,1), \qquad x_j=\exp(v)z_j \]

elimina la escala condicional cambiante de las coordenadas que se muestrean directamente (McElreath 2020, 421-23). Bayesian Workflow desarrolla la misma patología en configuraciones jerárquicas de mayor dimensión y enfatiza que la dificultad puede persistir cuando la verosimilitud es demasiado débil para dominar la geometría de la previa con forma de embudo (Gelman et al. 2026, sec. 12.3).

7.16 Parametrización centrada y no centrada

Regresemos al modelo jerárquico general:

\[ \alpha_j\mid\mu_\alpha,\tau_\alpha \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

7.16.1 Forma centrada

En la forma centrada, \(\alpha_j\) aparece directamente como parámetro:

\[ \alpha_j=\mu_\alpha+\varepsilon_j, \qquad \varepsilon_j\sim\mathcal N(0,\tau_\alpha^2). \]

La posterior puede inducir una dependencia fuerte entre \(\alpha_j\) y \(\tau_\alpha\) cuando los datos individuales de cada grupo son débiles.

7.16.2 Forma no centrada

Definimos

\[ z_j\sim\mathcal N(0,1), \]

\[ \alpha_j=\mu_\alpha+\tau_\alpha z_j. \]

Ahora el parámetro básico \(z_j\) tiene una escala fija que no depende de \(\tau_\alpha\) en la previa. La correlación posterior no desaparece siempre, pero en regímenes de información grupal débil suele reducirse de manera sustancial.

Matsuura desarrolla esta reparametrización para modelos jerárquicos en la sección 9.3.2 y advierte que no constituye una mejora universal: la conveniencia depende de cuánta información aporten los datos sobre los coeficientes específicos de grupo.

7.16.3 El modelo no cambia, las coordenadas sí

Para cualquier \(\mu_\alpha\) y \(\tau_\alpha>0\),

\[ z_j\sim\mathcal N(0,1) \quad\Longrightarrow\quad \mu_\alpha+\tau_\alpha z_j \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

Por tanto, si transformamos correctamente las variables, las dos parametrizaciones implican la misma distribución sobre \(\alpha_j\) y las mismas predicciones. La diferencia es la geometría en el espacio que recorre el sampler.

7.17 ¿Cuándo conviene centrar y cuándo no?

No existe una regla “no centrado siempre”. Una guía cualitativa es comparar la información proveniente de la distribución poblacional con la información que cada grupo aporta mediante su verosimilitud.

Situación	Tendencia computacional frecuente
Pocos datos por grupo o señal grupal débil	No centrada suele funcionar mejor
\(\tau\) cercana a cero y coeficientes grupales poco informados	No centrada suele aliviar el funnel
Muchos datos informativos por grupo	Centrada puede ser más eficiente
Algunos componentes bien informados y otros no	Puede ser útil una parametrización mixta

McElreath subraya que hay situaciones en las que la forma centrada es mejor (McElreath 2020, 425). Bayesian Workflow también muestra que la elección puede depender de cada bloque jerárquico y que en modelos grandes puede ser conveniente mezclar las dos formas (Gelman et al. 2026, sec. 12.3 y estudios de caso).

Notabrms y la parametrización

brms genera código Stan y utiliza estrategias de parametrización diseñadas para modelos multinivel. Para comprender la diferencia entre formas centradas y no centradas utilizaremos Stan explícito en un experimento pequeño. En el trabajo aplicado seguiremos utilizando brms con backend = "cmdstanr".

7.18 Identificación estadística y dificultad computacional

Una dificultad computacional no debe confundirse con falta de identificación.

Considere

\[ y_i \sim \mathcal N\left(a+(b+c)x_i,\sigma^2\right). \tag{7.8}\]

La verosimilitud depende de \(b\) y \(c\) únicamente mediante la suma

\[ \delta=b+c. \]

Si reemplazamos

\[ b^*=b+d, \qquad c^*=c-d, \]

entonces

\[ b^*+c^*=b+c, \]

y el modelo para los datos observados no cambia. Los datos pueden informar \(\delta\), pero no distinguir por separado \(b\) y \(c\) sin estructura adicional. Matsuura utiliza esencialmente esta construcción en la sección 9.1 para ilustrar parámetros no identificados.

7.18.1 Una previa propia puede regularizar, pero no determina qué información aportaron los datos

Supongamos que asignamos

\[ b \sim \mathcal N(0,1), \qquad c \sim \mathcal N(0,1). \]

Estas previas pueden hacer que la distribución posterior sea propia y computacionalmente más estable, porque restringen los valores plausibles de \(b\) y \(c\). Sin embargo, esto no implica que los datos permitan distinguir claramente ambos coeficientes.

Si la verosimilitud depende principalmente de una combinación como \(b+c\), distintas parejas de valores de \(b\) y \(c\) pueden ser igualmente compatibles con los datos siempre que produzcan aproximadamente la misma suma. En ese caso, la información necesaria para separar ambos coeficientes proviene parcial o totalmente de las distribuciones previas. Por ello, sería incorrecto afirmar que los datos observados identifican por separado a \(b\) y \(c\).

Esta distinción es especialmente importante en modelos multinivel. Por ejemplo, una desviación estándar entre grupos puede estar débilmente informada cuando hay pocos grupos o poca información dentro de ellos. Una previa razonable puede restringir valores extremos y facilitar el ajuste computacional, pero no crea información que los datos no contienen. En consecuencia, una posterior bien comportada desde el punto de vista computacional no debe interpretarse automáticamente como evidencia de que todos los parámetros están bien identificados por los datos.

El diagnóstico debe separar entonces dos preguntas:

1. ¿La posterior puede explorarse adecuadamente mediante el algoritmo de muestreo?
2. ¿Qué tanto de la información posterior sobre cada parámetro proviene de los datos y qué tanto de la previa?

Resolver la primera pregunta no resuelve necesariamente la segunda.

7.18.2 Una tabla diagnóstico para separar problemas

Señal	Posible causa	Investigación	Acción defendible
\(\widehat R\) alto y cadenas separadas	multimodalidad o mezcla muy lenta	trazas, rangos, inicialización, simplificar	entender modos; reparametrizar o reformular
ESS bajo, \(\widehat R\) aceptable	autocorrelación alta	trazas, correlaciones posteriores, escalas	reescalar, centrar, reparametrizar, luego más iteraciones si procede
divergencias	curvatura fuerte, funnel, restricciones	pares con divergencias, escalas, estructura jerárquica	paso menor, no centrado, regularización justificada, reformulación
max_treedepth frecuente	trayectorias muy largas	correlaciones, escalas, ESS por costo	centrar, reescalar, reparametrizar; aumentar límite solo si se justifica
coeficientes que se desplazan sin estabilizar	no identificación o posterior impropia	álgebra del predictor lineal, previas, simplificaciones	eliminar redundancia, imponer identificación sustantiva, usar previas propias
posterior sensible a una previa razonable	datos poco informativos	análisis de sensibilidad	comunicar sensibilidad; no ocultarla aumentando iteraciones

7.19 Un flujo de diagnóstico para modelos multinivel

Para los análisis del curso utilizaremos el siguiente orden.

7.19.1 Paso 1. Verificar la formulación y las escalas

Antes de culpar al sampler:

• comprobar unidades y transformaciones;
• revisar que desviaciones estándar estén restringidas a valores positivos;
• verificar índices de grupos;
• revisar el centrado de predictores;
• inspeccionar si hay parámetros redundantes;
• revisar que las previas sean propias y razonables.

7.19.2 Paso 2. Ajustar varias cadenas

Usar varias cadenas, semillas reproducibles y una fase de warmup suficiente. Durante la construcción temprana puede emplearse un ajuste corto para detectar problemas rápidamente; los resultados finales requieren una ejecución con precisión Monte Carlo adecuada.

7.19.3 Paso 3. Revisar diagnósticos generales

Examinar:

• \(\widehat R\);
• bulk-ESS;
• tail-ESS;
• MCSE;
• trazas o rangos cuando sean informativos.

7.19.4 Paso 4. Revisar diagnósticos de HMC

Examinar:

• divergencias;
• saturación de max_treedepth;
• longitud de trayectorias o número de pasos;
• advertencias numéricas.

7.19.5 Paso 5. Localizar la causa

Preguntar si el problema parece provenir de:

• falta de identificación;
• escalas inadecuadas;
• correlaciones posteriores fuertes;
• un funnel jerárquico;
• multimodalidad;
• una previa excesivamente amplia;
• un modelo demasiado flexible para la información disponible.

7.19.6 Paso 6. Modificar una cosa a la vez

Una modificación es útil para diagnóstico solo si podemos atribuirle el cambio observado. Entre las alternativas:

• centrar o reescalar;
• utilizar una parametrización no centrada;
• usar una previa débilmente informativa defendible;
• eliminar una redundancia;
• simplificar temporalmente el modelo para aislar el componente problemático;
• aumentar adapt_delta como experimento controlado;
• ejecutar más iteraciones cuando el problema sea solamente precisión Monte Carlo insuficiente.

7.19.7 Paso 7. Volver a comprobar el modelo

Una vez que el muestreo es confiable, todavía debemos realizar comprobaciones predictivas. Un modelo computacionalmente estable puede describir mal los datos.

7.20 Laboratorio reproducible en R

El laboratorio tendrá tres partes:

1. Monte Carlo y autocorrelación. Mostraremos por qué el número efectivo puede ser mucho menor que el número bruto de simulaciones.
2. El funnel y la parametrización. Visualizaremos la geometría del funnel, mostraremos cómo intentar la forma centrada de manera interactiva y ajustaremos de forma reproducible la representación no centrada con cmdstanr.
3. Aplicación multinivel con sleepstudy. Ajustaremos con brms y construiremos un diagnóstico reproducible de \(\widehat R\), ESS, MCSE, trazas y NUTS.

7.20.1 Experimento 1: autocorrelación y tamaño efectivo

Primero construiremos dos secuencias con la misma distribución marginal normal estándar y el mismo número de elementos. Una será aproximadamente independiente y la otra tendrá dependencia AR(1) fuerte. No son ajustes posteriores reales; son un experimento para aislar el papel de la dependencia serial.

simular_ar1 <- function(n, rho, semilla) {
  set.seed(semilla)

  x <- numeric(n)
  x[1] <- rnorm(1)

  for (s in 2:n) {
    x[s] <- rho * x[s - 1] +
      sqrt(1 - rho^2) * rnorm(1)
  }

  x
}

n_sim <- 4000

sim_ind <- simular_ar1(n_sim, rho = 0.00, semilla = 1653071)
sim_dep <- simular_ar1(n_sim, rho = 0.97, semilla = 1653072)

datos_mc_s7 <- bind_rows(
  tibble(
    iteracion = seq_len(n_sim),
    valor = sim_ind,
    serie = "Dependencia baja"
  ),
  tibble(
    iteracion = seq_len(n_sim),
    valor = sim_dep,
    serie = "Dependencia alta"
  )
)

Ambas series tienen aproximadamente la misma distribución marginal.

datos_mc_s7 |>
  group_by(serie) |>
  summarise(
    media = mean(valor),
    sd = sd(valor),
    .groups = "drop"
  )

# A tibble: 2 × 3
  serie               media    sd
  <chr>               <dbl> <dbl>
1 Dependencia alta -0.170   0.943
2 Dependencia baja -0.00397 1.00 

 datos_mc_s7 |>
  filter(iteracion <= 600) |>
  ggplot(aes(x = iteracion, y = valor)) +
  geom_line(linewidth = 0.4) +
  facet_wrap(~ serie, ncol = 1) +
  labs(
    x = "Iteración",
    y = "Valor simulado"
  ) +
  theme_minimal(base_size = 12)

Figura 7.1: Dos secuencias con distribuciones marginales semejantes pero dependencia serial muy distinta.

Ahora convertimos cada secuencia en un objeto de posterior y calculamos ESS y MCSE.

draws_ind <- posterior::as_draws_matrix(
  matrix(sim_ind, ncol = 1, dimnames = list(NULL, "theta"))
)

draws_dep <- posterior::as_draws_matrix(
  matrix(sim_dep, ncol = 1, dimnames = list(NULL, "theta"))
)

bind_rows(
  posterior::summarise_draws(
    draws_ind,
    "mean", "sd", "ess_bulk", "ess_tail", "mcse_mean"
  ) |>
    mutate(serie = "Dependencia baja"),
  posterior::summarise_draws(
    draws_dep,
    "mean", "sd", "ess_bulk", "ess_tail", "mcse_mean"
  ) |>
    mutate(serie = "Dependencia alta")
) |>
  select(
    serie,
    variable,
    mean,
    sd,
    ess_bulk,
    ess_tail,
    mcse_mean
  )

# A tibble: 2 × 7
  serie            variable     mean    sd ess_bulk ess_tail mcse_mean
  <chr>            <chr>       <dbl> <dbl>    <dbl>    <dbl>     <dbl>
1 Dependencia baja theta    -0.00397 1.00    3835.     3883.    0.0161
2 Dependencia alta theta    -0.170   0.943     90.2     201.    0.0992

El número bruto de valores es el mismo, pero la serie con dependencia fuerte contiene mucha menos información efectiva para estimar la media. Este experimento también muestra por qué no es suficiente reportar “4000 iteraciones”.

7.20.2 Experimento 2: visualizar directamente un funnel

Antes de ejecutar HMC, simulemos directamente desde Ecuación 7.7. Esto nos permite visualizar la distribución objetivo sin confundir sus propiedades con el comportamiento del algoritmo.

set.seed(1653073)

n_funnel <- 8000

funnel_directo <- tibble(
  v = rnorm(n_funnel, 0, 3)
) |>
  mutate(
    x1 = rnorm(n_funnel, 0, exp(v))
  )

funnel_directo |>
  ggplot(aes(x = v, y = x1)) +
  geom_point(alpha = 0.18, size = 0.8) +
  coord_cartesian(ylim = c(-25, 25)) +
  labs(
    x = "v",
    y = "x[1]"
  ) +
  theme_minimal(base_size = 12)

Figura 7.2: Simulación directa de una distribución tipo funnel. Cuando v es negativo, la escala condicional de x se vuelve muy pequeña.

La figura corresponde directamente al caso bidimensional que utilizaremos en el laboratorio: \(v\) y una sola coordenada \(x_1\). En un modelo jerárquico puede haber muchas coordenadas de este tipo; al aumentar su número, la geometría puede volverse aún más difícil. La parte estrecha representa el régimen en el que \(\exp(v)\) es muy pequeño.

7.20.3 Experimento 3: la forma centrada como experimento interactivo

La distribución del experimento anterior puede escribirse directamente en Stan como

\[ v\sim\mathcal N(0,3), \qquad x\sim\mathcal N\!\left(0,\exp(v)\right). \]

Esta es precisamente la clase de geometría para la cual HMC puede presentar dificultades severas. Hay un detalle importante para un capítulo reproducible: el objetivo de este ejemplo es provocar una distribución difícil. Dependiendo de la versión de Stan, de la adaptación y de detalles numéricos, el ajuste centrado puede terminar con divergencias y ESS muy pequeño, pero también puede ocurrir que una o varias cadenas aborten por completo. En este último caso no existen diagnósticos posteriores que puedan extraerse de esas cadenas.

Por esta razón, no ejecutaremos automáticamente el ajuste centrado durante el render del capítulo. El siguiente bloque queda como experimento interactivo. Esta decisión separa dos objetivos distintos:

• el documento debe poder compilar de manera reproducible;
• el estudiantado debe poder observar, de manera controlada, cómo una geometría patológica afecta HMC.

stan_funnel_c <- '
data {
  int<lower=1> J;
}
parameters {
  real v;
  vector[J] x;
}
model {
  v ~ normal(0, 3);
  x ~ normal(0, exp(v));
}
'

archivo_funnel_c <- file.path(
  "_stan",
  "s7_funnel_centered.stan"
)

writeLines(stan_funnel_c, archivo_funnel_c)

mod_funnel_c <- cmdstanr::cmdstan_model(
  archivo_funnel_c
)

fit_funnel_c <- mod_funnel_c$sample(
  data = list(J = 1L),
  seed = 1653074,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000,
  adapt_delta = 0.95,
  max_treedepth = 10,
  refresh = 100
)

fit_funnel_c$summary(
  variables = c("v", "x[1]")
)

fit_funnel_c$diagnostic_summary()

AdvertenciaUn ajuste que no produce cadenas utilizables no tiene diagnósticos posteriores

Si todas las cadenas terminan con error, llamadas como

fit_funnel_c$sampler_diagnostics()

no pueden funcionar: no hay simulaciones terminadas sobre las cuales calcular divergencias, profundidad del árbol o número de pasos leapfrog. En ese caso el problema está antes del resumen diagnóstico. Aumentar iteraciones no resuelve automáticamente esta situación.

El resultado del experimento centrado no debe evaluarse por una cifra específica de divergencias. Tres resultados son informativos:

1. el ajuste termina, pero presenta divergencias;
2. el ajuste termina sin divergencias, pero con baja eficiencia o trayectorias muy costosas;
3. una o más cadenas no consiguen terminar.

Los tres resultados son compatibles con la idea que queremos estudiar: una representación matemáticamente válida puede inducir una geometría muy desfavorable para HMC.

7.20.4 Experimento 4: forma no centrada del mismo funnel

En lugar de muestrear directamente \(x\) condicionado a \(v\), introducimos una variable estandarizada

\[ z\sim\mathcal N(0,1) \]

y reconstruimos

\[ x=\exp(v)z. \]

La distribución conjunta inducida para \((v,x)\) es la misma, pero HMC explora las coordenadas \((v,z)\), cuya geometría es mucho más regular. Como \(x\) no interviene en la densidad objetivo, lo calculamos en generated quantities; así la transformación sirve para interpretar y graficar los draws sin formar parte de la dinámica de HMC. Esta distinción entre el modelo probabilístico y la representación computacional es el punto central de la parametrización no centrada.

stan_funnel_nc <- '
data {
  int<lower=1> J;
}
parameters {
  real v;
  vector[J] z;
}
model {
  v ~ normal(0, 3);
  z ~ std_normal();
}
generated quantities {
  vector[J] x;
  x = exp(v) * z;
}
'

archivo_funnel_nc <- file.path(
  "_stan",
  "s7_funnel_noncentered.stan"
)

lineas_nuevas_nc <- strsplit(
  stan_funnel_nc,
  "\n",
  fixed = TRUE
)[[1]]

if (
  !file.exists(archivo_funnel_nc) ||
    !identical(readLines(archivo_funnel_nc), lineas_nuevas_nc)
) {
  writeLines(lineas_nuevas_nc, archivo_funnel_nc)
}

mod_funnel_nc <- cmdstanr::cmdstan_model(
  archivo_funnel_nc
)

fit_funnel_nc <- mod_funnel_nc$sample(
  data = list(J = 1L),
  seed = 1653076,
  init = 0,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 1000,
  adapt_delta = 0.90,
  max_treedepth = 10,
  refresh = 0
)

Running MCMC with 4 parallel chains...

Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.2 seconds.

Construimos ahora una función para resumir los diagnósticos de NUTS. A diferencia del ajuste centrado opcional, este ajuste es el que utilizaremos para el render reproducible del capítulo.

resumir_nuts_cmdstan <- function(fit, max_treedepth) {
  diag <- fit$sampler_diagnostics(format = "df") |>
    as.data.frame()

  tibble(
    divergencias = sum(diag$divergent__ == 1),
    treedepth_max_observado = max(diag$treedepth__),
    saturaciones_treedepth = sum(
      diag$treedepth__ >= max_treedepth
    ),
    mediana_n_leapfrog = median(diag$n_leapfrog__),
    max_n_leapfrog = max(diag$n_leapfrog__)
  )
}

resumen_funnel_nc <- resumir_nuts_cmdstan(
  fit_funnel_nc,
  max_treedepth = 10
)

resumen_funnel_nc

# A tibble: 1 × 5
  divergencias treedepth_max_observado saturaciones_treedepth mediana_n_leapfrog
         <int>                   <dbl>                  <int>              <dbl>
1            0                       3                      0                  3
# ℹ 1 more variable: max_n_leapfrog <dbl>

También examinamos \(\widehat R\) y los tamaños efectivos de muestra para \(v\), \(z\) y la cantidad transformada \(x\).

fit_funnel_nc$summary(
  variables = c("v", "z[1]", "x[1]")
) |>
  select(
    variable,
    mean,
    sd,
    rhat,
    ess_bulk,
    ess_tail
  )

# A tibble: 3 × 6
  variable    mean     sd  rhat ess_bulk ess_tail
  <chr>      <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 v        -0.0222   3.01  1.00    2832.    2473.
2 z[1]     -0.0399   1.05  1.00    2745.    2301.
3 x[1]     -1.30   758.    1.00    2897.    2236.

La parametrización no centrada no cambia la distribución objetivo de \(v\) y \(x\). Cambia las coordenadas en las que HMC realiza la exploración. En las coordenadas \((v,z)\) la posterior del experimento es simplemente el producto de distribuciones normales independientes, por lo que desaparece la garganta estrecha que complica la dinámica hamiltoniana en las coordenadas \((v,x)\).

Para verificar que la transformación reconstruye la geometría original, graficamos las simulaciones posteriores de \((v,x)\) obtenidas a partir del ajuste no centrado.

draws_funnel_nc <- fit_funnel_nc$draws(
  variables = c("v", "x[1]"),
  format = "df"
) |>
  as.data.frame()

draws_funnel_nc |>
  ggplot(aes(x = v, y = .data[["x[1]"]])) +
  geom_point(alpha = 0.25, size = 0.8) +
  coord_cartesian(ylim = c(-25, 25)) +
  labs(
    x = "v",
    y = "x[1]"
  ) +
  theme_minimal(base_size = 12)

Figura 7.3: Simulaciones de (v,x) obtenidas mediante la parametrización no centrada. La distribución transformada conserva la forma de funnel aunque HMC explore las coordenadas más regulares (v,z).

Tip¿Dónde quedó la comparación de adapt_delta?

El ajuste centrado del funnel es deliberadamente inestable: utilizarlo como dependencia obligatoria del render hace que la compilación del capítulo dependa de que un ejemplo patológico consiga terminar. Por eso la comparación entre adapt_delta = 0.95 y adapt_delta = 0.99 se deja como ejercicio interactivo. Si el ajuste centrado termina, compare divergencias, profundidad del árbol y número de pasos leapfrog. Si no termina, ese resultado también debe documentarse: adapt_delta no sustituye una parametrización apropiada.

7.20.5 Experimento 5: una cresta de no identificación

No necesitamos ejecutar MCMC para reconocer la no identificación de Ecuación 7.8. Construiremos una superficie de log verosimilitud para \(b\) y \(c\) manteniendo fijos los demás parámetros.

set.seed(1653077)

n_id <- 40
x_id <- seq(-2, 2, length.out = n_id)

a_true <- 1
suma_true <- 2
sigma_true <- 0.6

y_id <- rnorm(
  n_id,
  mean = a_true + suma_true * x_id,
  sd = sigma_true
)

grid_id <- tidyr::crossing(
  b = seq(-2, 6, length.out = 180),
  c = seq(-2, 6, length.out = 180)
) |>
  mutate(
    loglik = purrr::map2_dbl(
      b,
      c,
      ~ sum(
        dnorm(
          y_id,
          mean = a_true + (.x + .y) * x_id,
          sd = sigma_true,
          log = TRUE
        )
      )
    )
  )

grid_id |>
  ggplot(aes(x = b, y = c, z = loglik)) +
  geom_contour(bins = 20) +
  geom_abline(
    slope = -1,
    intercept = suma_true,
    linetype = 2
  ) +
  labs(
    x = "b",
    y = "c"
  ) +
  coord_equal() +
  theme_minimal(base_size = 12)

Figura 7.4: Superficie de log verosimilitud para dos coeficientes que entran al predictor solo mediante su suma. La cresta muestra que múltiples pares (b,c) producen el mismo ajuste.

Aumentar adapt_delta no cambia esta superficie. Tampoco lo hace ejecutar un millón de iteraciones. Para separar \(b\) y \(c\) necesitamos eliminar la redundancia o introducir información adicional que distinga sus papeles.

7.21 Aplicación: diagnóstico de un modelo de trayectorias con sleepstudy

Retomaremos sleepstudy, utilizado en la semana 6 para discutir previas. El objetivo aquí no es volver a elegir las previas desde cero, sino completar un informe computacional reproducible.

Los datos contienen mediciones repetidas de tiempo de reacción durante días de privación de sueño. Utilizaremos interceptos y pendientes variables por sujeto. Para que el intercepto represente una condición dentro del rango observado, centraremos Days en su media.

datos_s7 <- lme4::sleepstudy |>
  as_tibble() |>
  mutate(
    Subject = factor(Subject),
    Days_c = Days - mean(Days)
  )

datos_s7 |>
  summarise(
    n = n(),
    n_subject = n_distinct(Subject),
    min_days = min(Days),
    max_days = max(Days),
    media_days = mean(Days),
    faltantes = sum(is.na(Reaction))
  )

# A tibble: 1 × 6
      n n_subject min_days max_days media_days faltantes
  <int>     <int>    <dbl>    <dbl>      <dbl>     <int>
1   180        18        0        9        4.5         0

datos_s7 |>
  ggplot(
    aes(
      x = Days,
      y = Reaction,
      group = Subject
    )
  ) +
  geom_line(alpha = 0.45) +
  geom_point(alpha = 0.75, size = 1.2) +
  labs(
    x = "Día de privación de sueño",
    y = "Tiempo de reacción (ms)"
  ) +
  theme_minimal(base_size = 12)

Figura 7.5: Trayectorias observadas de tiempo de reacción por sujeto.

7.21.1 Especificación del modelo

Escribimos

\[ Reaction_{ij} \sim \mathcal N(\mu_{ij},\sigma^2), \]

\[ \mu_{ij} = \alpha_j+\beta_j Days^c_{ij}, \]

con

\[ \begin{pmatrix} \alpha_j\\ \beta_j \end{pmatrix} \sim \mathcal N_2 \left( \begin{pmatrix} \mu_\alpha\\ \mu_\beta \end{pmatrix}, \Sigma \right). \]

En brms:

formula_s7 <- bf(
  Reaction ~ 1 + Days_c + (1 + Days_c | Subject)
)

formula_s7

Reaction ~ 1 + Days_c + (1 + Days_c | Subject) 

7.21.2 Previas

Usaremos previas explícitas coherentes con la escala trabajada en la semana 6. Debido al centrado, el intercepto corresponde aproximadamente a la mitad del período de observación, por lo que usamos una previa centrada en una magnitud de reacción plausible para esa región.

priors_s7 <- c(
  prior(normal(300, 100), class = "Intercept"),
  prior(normal(0, 20), class = "b", coef = "Days_c"),
  prior(
    exponential(0.02),
    class = "sd",
    group = "Subject",
    coef = "Intercept"
  ),
  prior(
    exponential(0.05),
    class = "sd",
    group = "Subject",
    coef = "Days_c"
  ),
  prior(lkj(2), class = "cor", group = "Subject"),
  prior(exponential(0.02), class = "sigma")
)

priors_s7

             prior     class      coef   group resp dpar nlpar   lb   ub tag
  normal(300, 100) Intercept                                   <NA> <NA>    
     normal(0, 20)         b    Days_c                         <NA> <NA>    
 exponential(0.02)        sd Intercept Subject                 <NA> <NA>    
 exponential(0.05)        sd    Days_c Subject                 <NA> <NA>    
            lkj(2)       cor           Subject                 <NA> <NA>    
 exponential(0.02)     sigma                                   <NA> <NA>    
 source
   user
   user
   user
   user
   user
   user

Estas previas se utilizan aquí como continuación del ejemplo de la semana 6, no como valores universales para cualquier modelo longitudinal.

7.21.3 Ajuste con brms y cmdstanr

Fijaremos explícitamente las decisiones computacionales que queremos documentar.

fit_s7 <- brm(
  formula = formula_s7,
  data = datos_s7,
  family = gaussian(),
  prior = priors_s7,
  backend = "cmdstanr",
  seed = 1653078,
  chains = 4,
  cores = 4,
  iter = 2500,
  warmup = 1000,
  control = list(
    adapt_delta = 0.95,
    max_treedepth = 12
  ),
  refresh = 0,
  file = file.path(
    "_fits",
    "s7_sleepstudy_diagnostico_cmdstanr"
  ),
  file_refit = "on_change"
)

El objeto ajustado contiene simulaciones posteriores y diagnósticos. El orden de revisión importa: primero comprobamos el muestreo; luego interpretamos parámetros y predicciones.

7.22 Diagnóstico general del ajuste sleepstudy

7.22.1 Resumen de brms

summary(fit_s7)

 Family: gaussian 
  Links: mu = identity 
Formula: Reaction ~ 1 + Days_c + (1 + Days_c | Subject) 
   Data: datos_s7 (Number of observations: 180) 
  Draws: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
         total post-warmup draws = 6000

Multilevel Hyperparameters:
~Subject (Number of levels: 18) 
                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept)            38.51      6.86    27.68    54.15 1.00     1416
sd(Days_c)                6.14      1.32     3.98     9.09 1.00     1971
cor(Intercept,Days_c)     0.61      0.18     0.19     0.88 1.00     2810
                      Tail_ESS
sd(Intercept)             2815
sd(Days_c)                3229
cor(Intercept,Days_c)     4080

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   298.27      9.37   279.84   316.66 1.00     1030     1771
Days_c       10.42      1.59     7.23    13.61 1.00     1607     2782

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    25.79      1.53    22.97    29.00 1.00     5594     4395

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Este resumen muestra \(\widehat R\) y ESS, pero conviene construir una tabla explícita para los parámetros principales.

variables_clave_s7 <- c(
  "b_Intercept",
  "b_Days_c",
  "sd_Subject__Intercept",
  "sd_Subject__Days_c",
  "cor_Subject__Intercept__Days_c",
  "sigma"
)

draws_s7 <- posterior::as_draws_df(fit_s7)

posterior::subset_draws(
  draws_s7,
  variable = variables_clave_s7
) |>
  posterior::summarise_draws(
    "mean", "sd", "rhat", "ess_bulk", "ess_tail", "mcse_mean", "mcse_sd"
  ) |>
  select(
    variable,
    mean,
    sd,
    rhat,
    ess_bulk,
    ess_tail,
    mcse_mean,
    mcse_sd
  )

# A tibble: 6 × 8
  variable                  mean    sd  rhat ess_bulk ess_tail mcse_mean mcse_sd
  <chr>                    <dbl> <dbl> <dbl>    <dbl>    <dbl>     <dbl>   <dbl>
1 b_Intercept            298.    9.37   1.00    1030.    1771.   0.291   0.175  
2 b_Days_c                10.4   1.59   1.00    1607.    2782.   0.0397  0.0241 
3 sd_Subject__Intercept   38.5   6.86   1.00    1416.    2815.   0.179   0.119  
4 sd_Subject__Days_c       6.14  1.32   1.00    1971.    3229.   0.0301  0.0204 
5 cor_Subject__Intercep…   0.614 0.175  1.00    2810.    4080.   0.00327 0.00243
6 sigma                   25.8   1.53   1.00    5594.    4395.   0.0208  0.0199 

La lectura sugerida es:

1. ¿hay alguna cantidad con \(\widehat R\) claramente por encima de 1.01?;
2. ¿algún bulk-ESS o tail-ESS es pequeño respecto al resto?;
3. ¿la MCSE es suficientemente pequeña para los dígitos que pretendemos reportar?;
4. ¿los problemas se concentran en parámetros de escala o correlación, como suele ocurrir en modelos jerárquicos?

7.22.2 Gráficos de trazas

bayesplot::mcmc_trace(
  as.array(fit_s7),
  pars = variables_clave_s7
)

Figura 7.6: Trazas MCMC para parámetros poblacionales y de covarianza seleccionados.

No buscamos una textura visual específica. Buscamos cadenas que recorran regiones semejantes sin desplazamientos persistentes ni separación entre ellas.

7.22.3 Rangos por cadena

bayesplot::mcmc_rank_overlay(
  as.array(fit_s7),
  pars = variables_clave_s7
)

Figura 7.7: Distribución de rangos por cadena para parámetros seleccionados.

Una concentración sistemática de una cadena en ciertos rangos sugeriría que las cadenas no están explorando la misma distribución marginal.

7.23 Diagnósticos de NUTS en brms

No accederemos directamente a fit_s7$fit$..., porque el tipo interno de ese objeto depende del backend. brms::nuts_params() ofrece una interfaz común y evita depender de que el objeto subyacente sea S4 o CmdStanMCMC.

nuts_s7 <- brms::nuts_params(fit_s7)

nuts_s7 |>
  count(Parameter)

      Parameter    n
1 accept_stat__ 6000
2   treedepth__ 6000
3    stepsize__ 6000
4   divergent__ 6000
5  n_leapfrog__ 6000
6      energy__ 6000

7.23.1 Divergencias

n_divergencias_s7 <- sum(
  nuts_s7$Parameter == "divergent__" &
    nuts_s7$Value == 1
)

tibble(
  diagnostico = "Divergencias post-warmup",
  valor = n_divergencias_s7
)

# A tibble: 1 × 2
  diagnostico              valor
  <chr>                    <int>
1 Divergencias post-warmup     0

Una ejecución final deseable no presenta divergencias. Si aparecen, la respuesta no es eliminarlas del archivo de draws. Debemos investigar los parámetros asociados y revisar la geometría.

7.23.2 Profundidad del árbol

Como fijamos max_treedepth = 12, podemos contar cuántas iteraciones alcanzaron ese límite.

profundidades_s7 <- nuts_s7 |>
  filter(Parameter == "treedepth__") |>
  pull(Value)

tibble(
  profundidad_maxima_observada = max(profundidades_s7),
  iteraciones_en_limite = sum(profundidades_s7 >= 12)
)

# A tibble: 1 × 2
  profundidad_maxima_observada iteraciones_en_limite
                         <dbl>                 <int>
1                            6                     0

Si hubiera muchas saturaciones pero ninguna divergencia, investigaríamos eficiencia y correlación antes de aumentar el límite.

7.23.3 Número de pasos leapfrog

leapfrog_s7 <- nuts_s7 |>
  filter(Parameter == "n_leapfrog__") |>
  pull(Value)

tibble(
  mediana = median(leapfrog_s7),
  q90 = quantile(leapfrog_s7, 0.90),
  q99 = quantile(leapfrog_s7, 0.99),
  maximo = max(leapfrog_s7)
)

# A tibble: 1 × 4
  mediana   q90   q99 maximo
    <dbl> <dbl> <dbl>  <dbl>
1      31    63    63    127

Estos valores no tienen un umbral universal. Sirven para comparar versiones del mismo modelo o detectar trayectorias inusualmente costosas.

7.24 Un resumen computacional compacto

Construimos una tabla que podría incorporarse al informe del proyecto.

resumen_param_s7 <- posterior::subset_draws(
  draws_s7,
  variable = variables_clave_s7
) |>
  posterior::summarise_draws(
    "rhat", "ess_bulk", "ess_tail", "mcse_mean"
  ) |>
  summarise(
    max_rhat = max(rhat, na.rm = TRUE),
    min_ess_bulk = min(ess_bulk, na.rm = TRUE),
    min_ess_tail = min(ess_tail, na.rm = TRUE),
    max_mcse_media = max(mcse_mean, na.rm = TRUE)
  )

resumen_nuts_s7 <- tibble(
  divergencias = n_divergencias_s7,
  saturaciones_treedepth = sum(profundidades_s7 >= 12),
  profundidad_maxima = max(profundidades_s7),
  mediana_n_leapfrog = median(leapfrog_s7)
)

bind_cols(
  tibble(
    cadenas = 4,
    iteraciones_totales_por_cadena = 2500,
    warmup_por_cadena = 1000,
    adapt_delta = 0.95,
    max_treedepth_configurado = 12
  ),
  resumen_param_s7,
  resumen_nuts_s7
)

# A tibble: 1 × 13
  cadenas iteraciones_totales_por_cadena warmup_por_cadena adapt_delta
    <dbl>                          <dbl>             <dbl>       <dbl>
1       4                           2500              1000        0.95
# ℹ 9 more variables: max_treedepth_configurado <dbl>, max_rhat <dbl>,
#   min_ess_bulk <dbl>, min_ess_tail <dbl>, max_mcse_media <dbl>,
#   divergencias <int>, saturaciones_treedepth <int>, profundidad_maxima <dbl>,
#   mediana_n_leapfrog <dbl>

Una tabla así debe acompañarse de una interpretación. No basta con escribir “\(\widehat R=1\)”. El informe debería explicar si los diagnósticos son coherentes entre sí y si hubo decisiones correctivas.

7.25 ¿Qué haríamos si el ajuste presentara divergencias?

Supongamos que n_divergencias_s7 > 0. Un flujo razonable sería:

1. localizar si las divergencias se concentran en combinaciones de parámetros de escala y coeficientes grupales;
2. revisar si el predictor está centrado y razonablemente escalado;
3. revisar la previa para desviaciones estándar y correlaciones;
4. comprobar si la estructura de pendientes variables está informada por suficientes observaciones dentro de cada grupo;
5. probar una parametrización alternativa si trabajamos directamente en Stan;
6. aumentar adapt_delta de manera controlada y volver a diagnosticar;
7. documentar qué cambio eliminó o redujo el problema y si alteró las cantidades sustantivas de interés.

La conclusión no sería “subimos adapt_delta hasta que Stan dejó de quejarse”. Queremos una explicación de la causa probable y del efecto de la modificación.

7.26 ¿Qué haríamos si solamente el ESS fuera bajo?

Si \(\widehat R\) está cerca de 1, no hay divergencias ni saturaciones importantes y la geometría parece estable, pero el ESS de una cantidad es insuficiente para la precisión requerida, ejecutar más iteraciones puede ser la solución correcta.

Por ejemplo, si estimamos una media posterior con desviación 10 y ESS 100, Ecuación 7.4 sugiere aproximadamente

\[ MCSE \approx \frac{10}{\sqrt{100}} =1. \]

Si necesitamos una precisión Monte Carlo de aproximadamente 0.25, y la eficiencia por iteración se mantiene, necesitaríamos alrededor de

\[ S_{\text{eff}} \approx \left(\frac{10}{0.25}\right)^2 =1600 \]

simulaciones efectivas para esa cantidad.

Esta es una justificación concreta para ejecutar más iteraciones. Es distinta de prolongar indefinidamente una cadena que no mezcla entre modos.

7.27 Comprobación predictiva posterior: un diagnóstico diferente

Una vez que la aproximación computacional es confiable, realizamos una comprobación predictiva mínima.

pp_check(
  fit_s7,
  type = "dens_overlay",
  ndraws = 80
)

Figura 7.8: Comprobación predictiva posterior global para el modelo sleepstudy.

Esta figura responde a una pregunta diferente de \(\widehat R\):

• \(\widehat R\), ESS y divergencias preguntan si aproximamos adecuadamente la posterior del modelo;
• la comprobación predictiva pregunta si ese modelo ajustado reproduce características relevantes de los datos.

Podemos tener un ajuste sin ninguna advertencia de HMC y una comprobación predictiva deficiente. También podemos tener un modelo sustantivamente razonable cuya parametrización inicial produzca divergencias. Ambos problemas deben resolverse, pero no son el mismo problema.

7.28 Predicción y parametrización

Si dos parametrizaciones son verdaderamente equivalentes y ambas se muestrean correctamente, deben producir la misma distribución posterior sobre cantidades observables y, por tanto, las mismas predicciones dentro del error Monte Carlo.

Esta propiedad ofrece una comprobación útil cuando escribimos Stan directamente:

1. ajustar forma centrada y no centrada;
2. transformar ambas a una cantidad común, por ejemplo \(\alpha_j\);
3. comparar resúmenes y predicciones;
4. atribuir diferencias sustanciales persistentes a error Monte Carlo, problemas de muestreo o a una implementación que no es realmente equivalente.

La eficiencia puede cambiar mucho aun cuando las predicciones sean las mismas.

7.29 El papel de las previas en el diagnóstico computacional

La semana 6 mostró que una previa se justifica por la escala y el conocimiento disponible. Esta semana añadimos una consecuencia: previas extremadamente amplias pueden permitir regiones de parámetros con geometría difícil y predicciones absurdas.

Esto no significa que debamos escoger una previa porque produce cero divergencias. Una previa más restrictiva añade información al modelo. Bayesian Workflow advierte específicamente contra usar restricciones únicamente para hacer más rápido el algoritmo; si una previa más fuerte se utiliza para mejorar estabilidad, debe ser defendible sustantivamente y conviene evaluar sensibilidad (Gelman et al. 2026, sec. 12.4).

La secuencia correcta es:

\[ \text{conocimiento y escala} \rightarrow \text{previa defendible} \rightarrow \text{comprobación previa} \rightarrow \text{ajuste} \rightarrow \text{diagnóstico computacional}. \]

No debemos invertirla y ajustar la previa hasta que desaparezca una advertencia.

7.30 Problemas computacionales frente a problemas de identificación

Conviene resumir la distinción.

7.30.1 Problema computacional

La posterior está bien definida, pero el algoritmo la recorre de manera ineficiente o poco confiable debido a su geometría. Ejemplos:

• correlación posterior extrema;
• funnel;
• escalas muy diferentes;
• tamaño de paso inadecuado;
• trayectorias demasiado largas.

Una reparametrización puede cambiar radicalmente el rendimiento sin cambiar el modelo sobre las cantidades observables.

7.30.2 Problema de identificación

Diferentes valores de parámetros producen el mismo, o casi el mismo, comportamiento para los datos. Ejemplos:

• parámetros redundantes como \(b\) y \(c\) en Ecuación 7.8;
• intercepto y todos los indicadores de categoría incluidos sin una restricción identificadora;
• componentes de un modelo cuya separación no está sustentada por el diseño o los datos.

Más iteraciones no crean información ausente. Una previa puede regularizar o identificar la posterior en un sentido bayesiano, pero entonces la separación de parámetros depende de esa información previa.

7.30.3 Identificación débil

Entre ambos extremos existe la identificación débil: los datos distinguen formalmente los parámetros, pero solo débilmente. Este régimen es común en modelos multinivel con pocos grupos o pocas observaciones por grupo. Aquí regularización, reparametrización y diseño del estudio interactúan estrechamente.

7.31 Errores frecuentes

7.31.1 1. Declarar “convergencia” porque \(\widehat R=1.00\)

\(\widehat R\) es un diagnóstico necesario pero no suficiente. Debe leerse junto con ESS, trazas, diagnósticos de HMC y estructura del modelo.

7.31.2 2. Reportar únicamente el número bruto de iteraciones

“Usamos 4000 draws” no informa cuántas simulaciones efectivas se obtuvieron ni cuál fue la MCSE.

7.31.3 3. Confundir incertidumbre posterior con MCSE

La primera es parte del resultado inferencial; la segunda es error numérico debido a una simulación finita.

7.31.4 4. Eliminar draws divergentes

Una divergencia no es un dato atípico. Señala que una región de la posterior puede no haberse explorado correctamente.

7.31.5 5. Aumentar adapt_delta automáticamente

Puede ayudar cuando el problema es un tamaño de paso demasiado grande, pero puede ocultar que la geometría necesita reparametrización o que el modelo está débilmente identificado.

7.31.6 6. Aumentar max_treedepth como primera respuesta

Permite trayectorias más largas y costosas. No reduce una correlación posterior causada por un predictor mal centrado.

7.31.7 7. Ejecutar más iteraciones ante multimodalidad severa

Si las cadenas permanecen en modos distintos, duplicar una longitud moderada puede no lograr mezcla entre modos. Primero hay que entender la estructura posterior.

7.31.8 8. Suponer que no centrado siempre es mejor

La eficiencia depende de cuánta información aporta cada grupo. Con datos muy informativos por unidad, la forma centrada puede ser superior.

7.31.9 9. Confundir reparametrización con cambiar el modelo

Una reparametrización correcta conserva la distribución sobre las cantidades sustantivas. Cambiar una previa o eliminar un componente sí cambia el modelo.

7.31.10 10. Usar una previa más fuerte solo para eliminar advertencias

Si una previa modifica la posterior, esa información debe justificarse y comunicarse. La computación no es una justificación sustantiva suficiente.

7.31.11 11. Interpretar falta de identificación como “Stan necesita más tiempo”

Cuando dos parámetros entran solo mediante una combinación redundante, la dificultad es estructural.

7.31.12 12. Diagnosticar solo coeficientes poblacionales

En modelos multinivel los problemas pueden concentrarse en desviaciones estándar, correlaciones o coeficientes de grupos pequeños.

7.31.13 13. Revisar únicamente trazas

Las trazas son útiles, pero no sustituyen \(\widehat R\), ESS y los diagnósticos específicos de HMC.

7.31.14 14. Ignorar la unidad de las variables

Predictores en escalas muy distintas pueden perjudicar la adaptación y producir correlaciones difíciles. Reescalar también mejora la especificación de previas.

7.31.15 15. Confundir ausencia de advertencias con buen modelo

El sampler puede representar perfectamente la posterior de un modelo que falla en las comprobaciones predictivas.

7.31.16 16. Usar acceso directo a objetos internos de brms

Expresiones como fit$fit$diagnostic_summary() dependen del backend y de la clase del objeto interno. Para diagnósticos NUTS en un brmsfit, utilice brms::nuts_params() u otras interfaces documentadas; para un objeto creado directamente con cmdstanr, use sus métodos de CmdStanMCMC.

7.32 Hito del proyecto de curso

Al finalizar esta semana, cada grupo debe poder producir un diagnóstico computacional básico de la primera versión de su modelo. El informe puede ser breve, pero debe permitir reproducir y evaluar las decisiones.

Incluya al menos:

1. fórmula del modelo y estructura de agrupamiento;
2. versiones de R, brms, cmdstanr y CmdStan;
3. semilla;
4. número de cadenas;
5. iteraciones totales y de warmup;
6. adapt_delta y max_treedepth si se modificaron;
7. máximo \(\widehat R\) entre cantidades relevantes;
8. mínimo bulk-ESS y tail-ESS entre cantidades relevantes;
9. MCSE para al menos dos cantidades sustantivamente importantes;
10. número de divergencias;
11. número de iteraciones que alcanzaron la profundidad máxima configurada;
12. una figura de trazas o rangos cuando ayude a interpretar un problema;
13. descripción de cualquier modificación realizada para mejorar el ajuste;
14. evidencia de que la modificación fue una respuesta al problema identificado y no una búsqueda mecánica de advertencias cero;
15. una comprobación predictiva posterior mínima, claramente separada del diagnóstico computacional.

Un formato posible para la conclusión es:

Se ejecutaron cuatro cadenas con 1500 simulaciones post-warmup por cadena. Las cantidades principales presentaron \(\widehat R\) cercanos a 1, ESS suficientes para la precisión requerida y MCSE pequeñas frente a la incertidumbre posterior. No se observaron divergencias ni saturaciones de profundidad máxima. Las trazas no mostraron separación persistente entre cadenas. Por separado, la comprobación predictiva indicó […].

El texto debe completarse con los resultados reales del proyecto y no utilizarse como plantilla para afirmar diagnósticos que no se comprobaron.

7.33 Síntesis

Los modelos bayesianos complejos dependen de una capa computacional que debe ser diagnosticada. Las ideas centrales de la semana son:

1. MCMC aproxima integrales posteriores mediante simulaciones dependientes.
2. El número bruto de iteraciones no mide por sí solo la precisión; ESS y MCSE cuantifican la eficiencia y el error Monte Carlo.
3. Varias cadenas permiten comparar exploraciones iniciadas desde regiones diferentes.
4. \(\widehat R\) evalúa discrepancias entre y dentro de cadenas; valores cercanos a 1 son evidencia favorable, no una certificación absoluta.
5. HMC utiliza gradientes para construir trayectorias que evitan el comportamiento de paseo aleatorio de propuestas locales.
6. NUTS adapta la longitud de las trayectorias y Stan adapta durante warmup el tamaño de paso y características de la métrica.
7. Las divergencias señalan regiones donde la integración hamiltoniana puede ser poco confiable.
8. Alcanzar max_treedepth repetidamente señala trayectorias costosas y suele motivar una investigación de geometría y eficiencia.
9. Un funnel aparece cuando la escala de unos parámetros depende fuertemente de otros, situación frecuente en modelos jerárquicos débilmente informados.
10. La parametrización no centrada conserva el modelo probabilístico pero cambia las coordenadas del muestreo.
11. La forma no centrada suele ser útil cuando hay poca información por grupo; la centrada puede ser mejor cuando los coeficientes grupales están fuertemente informados por los datos.
12. Un problema de identificación no se corrige ejecutando más iteraciones.
13. Reescalamiento, regularización y reparametrización pueden mejorar la computación, pero cada modificación debe tener una interpretación estadística clara.
14. Un ajuste computacionalmente fiable todavía debe someterse a comprobación predictiva.

El principio operativo puede resumirse así:

\[ \boxed{ \text{una advertencia computacional es información sobre el ajuste, no ruido que debemos silenciar.} } \]

7.34 Ejercicios

7.34.1 Ejercicios conceptuales

7.34.1.1 Ejercicio 1. ESS y MCSE

Dos ajustes producen 4000 draws para un parámetro \(\theta\).

• Ajuste A: \(\operatorname{sd}(\theta\mid y)=4\), ESS = 1600.
• Ajuste B: \(\operatorname{sd}(\theta\mid y)=4\), ESS = 100.

1. Aproximar la MCSE de la media posterior en cada caso.
2. Explicar por qué el número bruto de draws no permite comparar la precisión Monte Carlo.
3. Si el objetivo es reportar la media con error Monte Carlo menor que 0.2, ¿cuál ajuste satisface el objetivo aproximadamente?

7.34.1.2 Ejercicio 2. \(\widehat R\) no es una prueba de validez

Considere cuatro cadenas que producen \(\widehat R=1.001\) para todos los parámetros.

1. ¿Qué evidencia aporta este resultado?
2. ¿Qué problemas podrían existir todavía?
3. Mencione al menos tres diagnósticos adicionales.
4. ¿Podemos concluir a partir de \(\widehat R\) que el modelo reproduce adecuadamente los datos? Explique.

7.34.1.3 Ejercicio 3. Divergencias y adapt_delta

Un modelo presenta 180 divergencias entre 4000 iteraciones post-warmup. Una persona propone cambiar únicamente adapt_delta de 0.8 a 0.999.

1. ¿Por qué ese cambio puede reducir divergencias?
2. ¿Por qué 180 divergencias sugieren investigar el modelo antes de confiar en esa modificación?
3. Enumere tres causas estructurales posibles.
4. Describa qué resultados compararía después de la modificación.

7.34.1.4 Ejercicio 4. Profundidad del árbol

Dos modelos no presentan divergencias.

• Modelo A: 0 de 4000 iteraciones alcanzan max_treedepth.
• Modelo B: 3500 de 4000 iteraciones alcanzan max_treedepth.

1. ¿Por qué no son computacionalmente equivalentes?
2. ¿Qué investigaría en el modelo B?
3. ¿Por qué aumentar el límite puede no ser la mejor primera acción?
4. ¿Qué papel podrían tener el centrado y el reescalamiento?

7.34.1.5 Ejercicio 5. Centrada versus no centrada

Para

\[ \alpha_j\sim\mathcal N(\mu,\tau^2), \]

1. escribir la parametrización no centrada;
2. demostrar que induce la misma distribución marginal condicional para \(\alpha_j\);
3. explicar por qué la geometría puede ser distinta;
4. indicar qué forma espera que funcione mejor con dos observaciones muy ruidosas por grupo;
5. indicar por qué su respuesta al punto anterior no constituye una regla universal.

7.34.1.6 Ejercicio 6. Identificación

Considere

\[ y_i\sim\mathcal N(a+b x_i+c x_i,\sigma^2). \]

1. Demuestre que la verosimilitud depende de \(b\) y \(c\) solo por su suma.
2. Proponga una reparametrización identificada usando \(\delta=b+c\).
3. Explique por qué aumentar el número de iteraciones no identifica \(b\) y \(c\) por separado.
4. Explique qué cambia si se asignan previas muy informativas y distintas a \(b\) y \(c\).

7.34.1.7 Ejercicio 7. Computación versus comprobación predictiva

Para cada afirmación, indique si corresponde principalmente a diagnóstico computacional, identificación o comprobación predictiva.

1. Las cadenas tienen \(\widehat R=1.3\).
2. El modelo reproduce medias globales pero no la variación entre grupos.
3. Dos coeficientes solo aparecen en el predictor como \(\beta_1+\beta_2\).
4. Hay divergencias concentradas cuando \(\tau\) se aproxima a cero.
5. El tail-ESS de un cuantil de interés es 40.
6. Las réplicas posteriores nunca producen la asimetría observada.

7.34.2 Ejercicios computacionales

7.34.2.1 Ejercicio 8. ESS bajo autocorrelación

Simule procesos AR(1) estacionarios con

\[ \rho\in\{0,0.5,0.9,0.99\} \]

y \(S=5000\).

1. Calcule bulk-ESS para cada serie usando posterior.
2. Grafique ESS contra \(\rho\).
3. Compare MCSE de la media.
4. Explique por qué la distribución marginal puede ser casi idéntica mientras la precisión Monte Carlo cambia drásticamente.

7.34.2.2 Ejercicio 9. adapt_delta en el funnel

Use el programa centrado del laboratorio.

1. Ajuste con adapt_delta = 0.8, 0.9, 0.95, 0.99.
2. Registre divergencias, mediana de n_leapfrog__, tiempo de ejecución y ESS de \(v\).
3. Grafique costo frente a divergencias.
4. Determine si aumentar adapt_delta elimina el problema o solo lo hace más costoso.
5. Compare con la parametrización no centrada.

7.34.2.3 Ejercicio 10. Número de grupos y parametrización

Simule datos de

\[ y_{ij}\sim\mathcal N(\alpha_j,1), \qquad \alpha_j\sim\mathcal N(0,\tau^2) \]

bajo combinaciones de:

• \(J\in\{5,20,80\}\);
• \(n_j\in\{2,20\}\);
• \(\tau\in\{0.1,1\}\).

Implemente versiones centrada y no centrada en Stan.

1. Compare divergencias y ESS.
2. Determine en qué regímenes cada parametrización es más eficiente.
3. Relacione el patrón con la cantidad de información por grupo.
4. Explique por qué este experimento es más informativo que memorizar una regla fija.

7.34.2.4 Ejercicio 11. Centrado de un predictor y correlación posterior

Simule una regresión donde \(x\) se encuentre alrededor de 2000:

\[ x_i\in[1990,2010]. \]

Ajuste dos modelos equivalentes:

1. con \(x\) original;
2. con \(x_c=x-2000\).

Compare:

• correlación posterior entre intercepto y pendiente;
• ESS;
• saturaciones de profundidad;
• tiempo de ejecución;
• predicciones en la escala original.

Explique por qué las predicciones pueden coincidir aunque el costo computacional sea distinto.

7.34.2.5 Ejercicio 12. Informe reproducible con brms

Seleccione uno de los modelos de las semanas 2–6 y produzca automáticamente una tabla con:

1. versión de brms;
2. versión de cmdstanr;
3. versión de CmdStan;
4. número de cadenas;
5. \(\widehat R\) máximo para parámetros seleccionados;
6. bulk-ESS mínimo;
7. tail-ESS mínimo;
8. divergencias;
9. profundidad máxima observada;
10. número de saturaciones del límite configurado.

Añada una interpretación de no más de 150 palabras.

7.34.2.6 Ejercicio 13. Sensibilidad entre prior y computación

Tome un modelo jerárquico con pocos grupos y compare dos previas sustantivamente defendibles para \(\tau\).

1. Realice comprobación predictiva previa para ambas.
2. Ajuste ambos modelos.
3. Compare divergencias y ESS.
4. Compare la posterior de \(\tau\) y las predicciones para un grupo nuevo.
5. Si una previa mejora el cómputo pero cambia las predicciones, explique por qué no es correcto elegirla únicamente por el diagnóstico computacional.

7.34.2.7 Ejercicio 14. Criticar una respuesta de software o IA

Se proporciona la siguiente recomendación:

“Stan reportó divergencias. Aumente adapt_delta a 0.999 y max_treedepth a 20. Si las advertencias desaparecen, el modelo convergió.”

Redacte una crítica técnica que:

1. identifique al menos cuatro problemas con la recomendación;
2. proponga un flujo de diagnóstico mejor;
3. distinga divergencias, profundidad del árbol e identificación;
4. explique qué evidencia sería necesaria antes de interpretar la posterior.

7.35 Referencias de la semana

Las lecturas centrales para esta unidad son:

• McElreath, Statistical Rethinking, capítulo 9 para MCMC, HMC, NUTS y diagnóstico de cadenas, y sección 13.4 para divergencias, el funnel y parametrización no centrada (McElreath 2020, chap. 9, pp. 420-425).
• Gelman et al., Bayesian Workflow, capítulo 11 para conjunto típico, inicialización, warmup, cadenas, \(\widehat R\), ESS y MCSE; capítulo 12 para modos de fallo, divergencias, profundidad del árbol, reparametrización y estrategias de corrección (Gelman et al. 2026, chaps. 11-12).
• Gelman et al., Bayesian Data Analysis, capítulo 11 para evaluación de simulaciones iterativas y capítulo 12, especialmente secciones 12.4–12.5, para HMC y su aplicación a modelos jerárquicos (Gelman et al. 2013, chaps. 11-12).
• Matsuura, Bayesian Statistical Modeling with Stan, R, and Python, capítulo 9, especialmente secciones 9.1–9.3 para no identificación, previas débilmente informativas y reparametrización de modelos jerárquicos.
• Bürkner para la conexión entre brms, Stan y HMC/NUTS, y para la implementación práctica de modelos multinivel bayesianos (Bürkner 2017).

Bürkner, Paul-Christian. 2017. «brms: An R Package for Bayesian Multilevel Models Using Stan». Journal of Statistical Software 80 (1): 1-28. https://doi.org/10.18637/jss.v080.i01.

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, y Donald B. Rubin. 2013. Bayesian Data Analysis. 3.ª ed. Chapman & Hall/CRC.

Gelman, Andrew, Aki Vehtari, Richard McElreath, et al. 2026. Bayesian Workflow. Chapman & Hall/CRC.

McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2.ª ed. Chapman & Hall/CRC.