8.4 Cuidado y alimentación de su cadena de Markov

No todas las cadenas de Markov se comportan bien. Esta sección cubre diagnósticos como R-hat, tamaño efectivo de muestra y gráficos de traza, y explica cómo abordar problemas comunes.

8.4.1 R-hat (diagnóstico de convergencia)

R-hat ( $\hat{R}$ ) compara la varianza dentro de las cadenas con la varianza entre cadenas. Si las cadenas han convergido, la varianza dentro y entre cadenas debería ser similar, y $\hat{R} \approx 1$ .

$\hat{R} \approx 1 \quad \text{(convergido)}$ $\hat{R} > 1.01 \quad \text{(no convergido)}$

# Demonstrate R-hat computation
set.seed(42)

# Simulate 4 chains from the same distribution (converged)

chains_good <- matrix(rnorm(4 * 1000, mean = 0, sd = 1), nrow = 1000, ncol = 4)

# Simulate 4 chains from different distributions (not converged)

chains_bad <- cbind(
rnorm(1000, mean = 0, sd = 1),
rnorm(1000, mean = 2, sd = 1),
rnorm(1000, mean = -1, sd = 1),
rnorm(1000, mean = 1, sd = 1.5)
)

# Compute R-hat

compute_rhat <- function(chains) {
m <- ncol(chains) # number of chains
n <- nrow(chains) # chain length

chain_means <- colMeans(chains)
chain_vars <- apply(chains, 2, var)

W <- mean(chain_vars) # within-chain variance
B <- n * var(chain_means) # between-chain variance

var_hat <- (n - 1) / n * W + B / n
rhat <- sqrt(var_hat / W)
return(rhat)
}

cat("R-hat for converged chains:", round(compute_rhat(chains_good), 4), "\n")
cat("R-hat for non-converged chains:", round(compute_rhat(chains_bad), 4), "\n")
cat("\nR-hat near 1 = converged\n")
cat("R-hat > 1.01 = not converged\n")

Click Run to execute this code in your browser with WebR

8.4.2 Tamaño efectivo de muestra (ESS)

El tamaño efectivo de muestra (ESS) nos dice cuántas muestras independientes valen nuestras extracciones autocorrelacionadas. Si la cadena tiene alta autocorrelación, ESS es mucho menor que el número real de extracciones.

$\text{ESS} \gg 1000 \quad \text{(bueno)}$

# Demonstrate effective sample size
set.seed(42)

# Independent samples (ESS = N)

independent <- rnorm(10000)

# Highly autocorrelated samples (AR(1) process)

rho <- 0.9 # autocorrelation
autocorr <- numeric(10000)
autocorr[1] <- rnorm(1)
for (t in 2:10000) {
autocorr[t] <- rho _ autocorr[t-1] + sqrt(1 - rho^2) _ rnorm(1)
}

# Compute autocorrelation

acf_indep <- acf(independent, plot = FALSE, lag.max = 50)
acf_auto <- acf(autocorr, plot = FALSE, lag.max = 50)

# ESS approximation: N / (1 + 2 * sum(autocorrelations))

ess*indep <- 10000 / (1 + 2 * sum(acf*indep$acf[-1]))
ess_auto <- 10000 / (1 + 2 * sum(acf_auto$acf[-1]))

cat("Effective sample size:\n")
cat(" Independent samples: ESS ≈", round(ess_indep, 0), "(out of 10000)\n")
cat(" Autocorrelated (rho=0.9): ESS ≈", round(ess_auto, 0), "(out of 10000)\n")

# Plot autocorrelation

par(mfrow = c(1, 2))
plot(acf_indep, main = "Independent: low autocorrelation")
plot(acf_auto, main = "Autocorrelated: high autocorrelation")
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

8.4.3 Gráficos de traza y transiciones divergentes

Un gráfico de traza muestra los valores de la cadena a lo largo de las iteraciones. Una cadena bien mezclada parece un "oruga peluda" — sin tendencias, sin períodos atascados. Una cadena mal mezclada muestra tendencias, períodos atascados, o una cadena explorando una región diferente que las demás.

Las transiciones divergentes ocurren en HMC cuando el muestreador encuentra una región de alta curvatura que no puede navegar. Las divergencias son advertencias de que la distribución a posteriori puede no estar completamente explorada.

# Demonstrate trace plots
set.seed(42)

# Good chain (well-mixed)

good_chain <- rnorm(2000)

# Bad chain (stuck / slow mixing)

bad_chain <- numeric(2000)
bad_chain[1] <- 0
for (t in 2:2000) {
bad_chain[t] <- bad_chain[t-1] + rnorm(1, sd = 0.05)
}

par(mfrow = c(2, 1))
plot(good_chain, type = "l", col = "steelblue",
main = "Good chain: well-mixed (hairy caterpillar)",
xlab = "Iteration", ylab = "Parameter")
plot(bad_chain, type = "l", col = "firebrick",
main = "Bad chain: slow mixing (random walk)",
xlab = "Iteration", ylab = "Parameter")
par(mfrow = c(1, 1))

cat("Good chain: looks like noise, no trends\n")
cat("Bad chain: shows trends, slow exploration\n")

Click Run to execute this code in your browser with WebR

8.4.4 Problemas comunes y soluciones

R-hat alto: Ejecute más iteraciones, reparametrice, o use distribuciones a priori más fuertes.
ESS bajo: Use HMC en lugar de caminata aleatoria, reparametrice, o ejecute más iteraciones.
Transiciones divergentes: Reparametrice (ej., centrar predictores), use parametrización no centrada, o ajuste el tamaño de paso.
Cadenas atascadas: Pruebe diferentes valores iniciales o reparametrice.

Rethinking

Confíe pero verifique. MCMC es poderoso, pero no es infalible. Los diagnósticos — R-hat, ESS, gráficos de traza, divergencias — son sus herramientas de verificación. Nunca reporte resultados de un modelo sin revisar estos diagnósticos. Un modelo con R-hat = 1.05 o 500 divergencias no es confiable, no importa lo buenos que parezcan los estimados. El golem no sabe cuándo ha fallado — usted debe verificar.

Overthinking

Parametrización no centrada. Una causa común de divergencias es un modelo jerárquico con geometría de "embudo". La parametrización no centrada soluciona esto muestreando una desviación de la media del grupo en lugar del parámetro directamente:

$\theta_j = \mu + \sigma \cdot \eta_j, \quad \eta_j \sim \text{Normal}(0, 1)$

en lugar de $\theta_j \sim \text{Normal}(\mu, \sigma)$ . Esto separa el muestreo de $\eta_j$ (fácil, normal estándar) del muestreo de $\mu$ y $\sigma$ (potencialmente correlacionados). Este es el truco más importante para hacer que HMC funcione en modelos jerárquicos.

Resumen de la Sección 8.4

R-hat ≈ 1 indica convergencia; R-hat > 1.01 indica problemas.
ESS le dice cuántas muestras independientes vale su cadena — busque ESS > 1000.
Los gráficos de traza deberían parecer una "oruga peluda" — sin tendencias ni períodos atascados.
Las transiciones divergentes señalan regiones que el muestreador no puede navegar — reparametrice para solucionar.
Siempre revise los diagnósticos antes de reportar resultados.

HMC fácil: map2stan

Resumen