8.4 Care and feeding of your Markov chain

Not every Markov chain behaves well. This section covers diagnostics such as R-hat, effective sample size, and trace plots, and explains how to address common problems.

8.4.1 R-hat (convergence diagnostic)

R-hat ( $\hat{R}$ ) compares the variance within chains to the variance between chains. If the chains have converged, within-chain and between-chain variance should be similar, and $\hat{R} \approx 1$ .

$\hat{R} \approx 1 \quad \text{(converged)}$ $\hat{R} > 1.01 \quad \text{(not converged)}$

# 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 Effective sample size (ESS)

Effective sample size (ESS) tells us how many independent samples our autocorrelated draws are worth. If the chain has high autocorrelation, ESS is much smaller than the actual number of draws.

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

# 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 Trace plots and divergent transitions

A trace plot shows the chain's values over iterations. A well-mixed chain looks like a "hairy caterpillar" — no trends, no stuck periods. A poorly mixed chain shows trends, stuck periods, or one chain exploring a different region than others.

Divergent transitions occur in HMC when the sampler encounters a region of high curvature that it cannot navigate. Divergences are warnings that the posterior may not be fully explored.

# 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 Common problems and solutions

High R-hat: Run more iterations, reparameterize, or use stronger priors.
Low ESS: Use HMC instead of random-walk, reparameterize, or run more iterations.
Divergent transitions: Reparameterize (e.g., center predictors), use non-centered parameterization, or adjust step size.
Stuck chains: Try different starting values or reparameterize.

Rethinking

Trust but verify. MCMC is powerful, but it is not infallible. The diagnostics — R-hat, ESS, trace plots, divergences — are your verification tools. Never report results from a model without checking these diagnostics. A model with R-hat = 1.05 or 500 divergences is not trustworthy, no matter how good the estimates look. The golem does not know when it has failed — you must check.

Overthinking

Non-centered parameterization. A common cause of divergences is a hierarchical model with "funnel" geometry. The non-centered parameterization fixes this by sampling a deviation from the group mean rather than the parameter directly:

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

instead of $\theta_j \sim \text{Normal}(\mu, \sigma)$ . This separates the sampling of $\eta_j$ (easy, standard normal) from the sampling of $\mu$ and $\sigma$ (potentially correlated). This is the single most important trick for getting HMC to work on hierarchical models.

Summary of Section 8.4

R-hat ≈ 1 indicates convergence; R-hat > 1.01 indicates problems.
ESS tells you how many independent samples your chain is worth — aim for ESS > 1000.
Trace plots should look like a "hairy caterpillar" — no trends or stuck periods.
Divergent transitions signal regions the sampler cannot navigate — reparameterize to fix.
Always check diagnostics before reporting results.

Easy HMC: map2stan

Summary