Chapter 5
Markov chain Monte Carlo

Lectures 9–10 · how the posterior is computed: Metropolis, Hamiltonian Monte Carlo and how to tell when a sampler can be trusted

Every model so far has ended with the same line: brm(...), then a summary. Behind that line, software drew thousands of samples from a posterior distribution that no one could write down in closed form. For the first four chapters that was fine: the models were small, the samplers behaved, and the diagnostics were all 1.00 and thousands. But samplers do not always behave, and a researcher who cannot tell when they have failed will eventually publish a result the computer never actually computed.

This chapter opens the black box. You will write two samplers yourself, in about twenty lines of R each, and run them on a real posterior: the probability of re-arrest by age in the COMPAS data of Chapter 4. The first, the Metropolis algorithm, is the original Markov chain Monte Carlo method, simple enough to explain to a child and slow enough to fail in instructive ways. The second, Hamiltonian Monte Carlo, is the engine inside Stan and brms; it borrows the physics of a skateboard rolling in a bowl to move through high-dimensional posteriors efficiently. Then you will learn to read the diagnostics that tell you whether a sampler has done its job, and see a real dataset break a sampler in a way that is entirely the model's fault.

• Lecture 9 is about how sampling works. Why grids fail in high dimensions; Metropolis by hand; why the way you write a model changes how hard it is to sample; and Hamiltonian Monte Carlo by hand.
• Lecture 10 is about trusting the result. Multiple chains, $\widehat{R}$, effective sample size and Monte Carlo error; what divergences mean; the "folk theorem" that computational problems are usually model problems; and how to check, with fake data, that a model can recover the truth.

Lecture 9 — How sampling works

5.1 Why not a grid?

In Section 3.2 we computed a posterior for two parameters by brute force: a grid of $150\times 150=22,500$ points, the posterior evaluated at each. For ten parameters, the same resolution would need ${150}^{10}\approx 6\times {10}^{21}$ points, more than the number of grains of sand on Earth. Grids are hopeless beyond three or four parameters, and real models have dozens or thousands.

There is a deeper problem, and it explains why every successful algorithm works the way it does. In high dimensions, almost none of the posterior's probability is near its peak. Consider a posterior that is a standard normal distribution in $d$ dimensions. Its density is highest at the centre, but there is very little volume near the centre, and a great deal of volume a little further out. Probability is density times volume, so the probability mass ends up concentrated in a thin shell at distance about $\sqrt{d}$ from the peak: the typical set (Betancourt, 2017). A sample from the posterior will almost never be near the mode.

In one dimension, the draws cluster near the peak. In ten, not one draw in a thousand is within distance 1 of it, and in a hundred the middle 98% of draws lie between distances of about 8.4 and 11.7: a thin shell around $\sqrt{100}=10$. This has two consequences. First, finding the mode (as optimisation does) tells you little about where the probability is; that is one reason point estimates can mislead in complex models. Second, an algorithm that explores the posterior must spend its time in the typical set, a thin, curved region that is hard to stay in by taking random steps. The rest of this lecture is about two ways of doing it.

5.2 The Metropolis algorithm

Here is the idea, as a story. A politician must visit a chain of islands, spending time on each in proportion to its population, but she does not know the populations and has no map. She can only ask the population of the island she is on and of one neighbour. Her rule: each day, toss a coin to pick a neighbouring island. If it has more people than the current one, go there. If it has fewer, go there with probability equal to the ratio of the populations (say, a neighbour with 60% of the current population is visited 60% of the time); otherwise stay. Remarkably, in the long run the share of days she spends on each island is exactly proportional to its population. Replace islands with parameter values and population with posterior density, and you have the Metropolis algorithm.

That last point is why Metropolis was revolutionary. Bayes' theorem divides by $p(y)$, an integral over all parameter values that is usually impossible to compute. Metropolis never needs it: it only compares two points at a time.

Now write it. Our target is the posterior for a two-parameter logistic regression from Chapter 4: the probability of re-arrest within two years as a function of age alone, for the 6,172 COMPAS defendants. Age is centred at 35 and measured in decades, so $a$ is the log-odds of re-arrest at age 35 and $b$ the change per decade.

library(tidyverse)
library(brms)
library(posterior)   # summaries and diagnostics of posterior draws
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

url <- "https://raw.githubusercontent.com/propublica/compas-analysis/master/compas-scores-two-years.csv"
compas <- read_csv(url, show_col_types = FALSE, name_repair = make.unique) |>
  filter(between(days_b_screening_arrest, -30, 30), is_recid != -1,
         c_charge_degree != "O", score_text != "N/A")
y   <- compas$two_year_recid             # 1 = re-arrested within two years
age <- compas$age
x_c <- (age - 35) / 10                   # age centred at 35, in decades

log_post <- function(theta, x, y) {      # log posterior, up to a constant
  eta <- theta[1] + theta[2] * x
  sum(y * eta - log1p(exp(eta))) +       # Bernoulli log-likelihood on the logit scale
    dnorm(theta[1], 0, 1.5, log = TRUE) + dnorm(theta[2], 0, 1, log = TRUE)   # priors
}
log_post(c(0, 0), x_c, y)

[1] -4280

metropolis <- function(n_iter, start, step, x, y) {
  draws <- matrix(NA, n_iter, 2, dimnames = list(NULL, c("a", "b")))
  current <- start; lp_current <- log_post(current, x, y); accepted <- 0
  for (t in 1:n_iter) {
    proposal <- current + rnorm(2, 0, step)                  # 1. propose a nearby point
    lp_proposal <- log_post(proposal, x, y)
    if (log(runif(1)) < lp_proposal - lp_current) {          # 2. accept with prob min(1, ratio)
      current <- proposal; lp_current <- lp_proposal; accepted <- accepted + 1
    }
    draws[t, ] <- current                                    # 3. record where we are
  }
  attr(draws, "accept") <- accepted / n_iter
  draws
}
set.seed(1)
chain_c <- metropolis(5000, start = c(1, 1), step = 0.03, x = x_c, y = y)
attr(chain_c, "accept")                                      # acceptance rate
kept <- chain_c[-(1:500), ]                                  # discard the first 500 (warmup)
rbind(mean = colMeans(kept), sd = apply(kept, 2, sd)) |> round(3)

[1] 0.475
          a      b
mean -0.207 -0.342
sd    0.027  0.023

The chain started at $(1,1)$, far from the posterior, and walked downhill into it within a few hundred steps; we discard those first 500 draws as warmup. The remaining 4,500 draws give a posterior mean of about −0.21 for $a$ (a 35-year-old's probability of re-arrest is ${\text{logit}}^{-1}(-0.21)\approx 45\%$) and about −0.34 for $b$: each decade of age lowers the log-odds by about a third. These agree with what brms gives for the same model (Section 5.5).

But 4,500 Metropolis draws are not worth 4,500 independent draws. Each draw is either a copy of the previous one (a rejection) or a small step away from it (an acceptance), so successive draws are strongly correlated. The effective sample size (ESS) measures how many independent draws the chain is worth for estimating a posterior mean. Play with the step size below and watch what it does to the chain.

Three regimes appear. With a tiny step size, nearly every proposal is accepted, but the chain crawls: it takes hundreds of steps to cross the posterior, and the trace looks like a slow random walk. With a huge step size, nearly every proposal lands far out in the tails and is rejected, so the chain sits still for long stretches and the trace looks like a staircase. In between (for this posterior, around $s=0.03$, with about half of proposals accepted) the chain moves well. Theory for well-behaved posteriors says the most efficient random-walk Metropolis accepts about 23% of proposals in many dimensions (Roberts, Gelman and Gilks, 1997) and about 44% in one (Gelman, Roberts and Gilks, 1996).

5.3 How you write the model changes how hard it is to sample

Now switch the widget to the raw parameterisation, with age in years and not centred. The model is mathematically the same (the predictions for every defendant are identical), but the posterior looks completely different: a long, thin, tilted ridge. The intercept is now the log-odds of re-arrest at age zero, a meaningless extrapolation that depends sharply on the slope: make the slope a little steeper and the intercept must rise a lot to keep the line through the data. In this parameterisation the posterior correlation between $a$ and $b$ is −0.95, and the slope's posterior standard deviation (0.0023 per year) is 35 times smaller than the intercept's.

set.seed(2)
chain_r <- metropolis(5000, start = c(0.99, -0.034), step = 0.03, x = age, y = y)
c(accept_centred = attr(chain_c, "accept"), accept_raw = attr(chain_r, "accept"))
c(ess_centred = ess_bulk(chain_c[-(1:500), "b"]), ess_raw = ess_bulk(chain_r[-(1:500), "b"]))

tibble(iteration = rep(1:5000, 2),
       slope_per_decade = c(chain_c[, "b"], 10 * chain_r[, "b"]),
       model = rep(c("centred: age in decades from 35", "raw: age in years"), each = 5000)) |>
  filter(iteration > 500) |>
  ggplot(aes(iteration, slope_per_decade)) +
  geom_line(colour = harbour[["sky"]], linewidth = 0.3) +
  facet_wrap(~ model, ncol = 1) +
  labs(x = "iteration", y = "slope per decade")

accept_centred     accept_raw 
        0.4750         0.0328 
ess_centred     ess_raw 
    522.649       1.457 

Even started at the right answer, the raw chain is nearly useless. With a step size that suits the intercept, nearly every proposal changes the slope by far too much and is rejected; with a step size that suits the slope, the intercept barely moves. It accepts only 3% of its proposals, and its effective sample size for the slope falls from about 520 to about 1. Nothing is wrong with the data or the model; the trouble is entirely in the coordinates. Centring and scaling predictors, which Section 3.3 recommended for interpretability, is also one of the most effective things you can do for computation. (brms centres predictors internally for exactly this reason, before translating the results back.)

5.4 Hamiltonian Monte Carlo

Metropolis explores by guessing: it proposes a random nearby point and checks it. In two dimensions that is fine. In a hundred, a random direction is almost always a bad one (it points out of the thin typical set), so proposals must be tiny and the chain diffuses at a snail's pace. A sampler that knew which way the posterior curves could do much better. Hamiltonian Monte Carlo (HMC) uses the gradient of the log posterior, which tells it the slope of the landscape at every point, to propose distant points that are still likely to be accepted.

The idea comes from physics. Turn the posterior upside down, so that high probability becomes low ground: a bowl whose shape is $-\log p(\theta \mid y)$. Put a frictionless skateboard at the current point and give it a random push (a random momentum). It rolls downhill, speeds up, climbs the far side, slows, turns and keeps rolling. After a while, stop it and record where it is. Because the board follows the shape of the bowl, it travels along the typical set rather than leaving it, and because its energy is conserved, it ends at a point of about the same total energy as it started, which means the proposal is almost always accepted.

A computer simulates the rolling in small discrete steps with the leapfrog integrator: half a step of momentum change using the gradient, a full step of position change using the momentum, and so on. With a small enough step size the simulated energy is almost exactly conserved; with too large a step the simulation becomes unstable and the energy explodes. The widget shows both.

With a moderate step size, the trajectory swings along the ridge and lands far from where it started, and the energy error is tiny, so the move is accepted: one HMC step can cross the whole posterior. Increase $\epsilon$ toward 0.5 and the trajectory starts to zigzag across the narrow direction of the ridge, the energy error grows, and acceptance collapses. Push it further and the trajectory flies off entirely: a divergence, which Section 5.7 explains. Notice also that the trajectory eventually turns back on itself. The No-U-Turn Sampler (NUTS) of Hoffman and Gelman (2014), used by Stan, stops the simulation automatically when it starts to turn back, so the user never has to choose $L$, and it tunes $\epsilon$ during warmup.

Here is HMC in R. It needs the gradient of the log posterior, which for logistic regression is simple: for the intercept, ${\sum }_i(y_i-p_i)$ minus the prior's pull; for the slope, ${\sum }_i(y_i-p_i)x_i$ minus the prior's pull.

grad_log_post <- function(theta, x, y) {
  p <- plogis(theta[1] + theta[2] * x)
  c(sum(y - p) - theta[1] / 1.5^2, sum((y - p) * x) - theta[2] / 1^2)
}
hmc <- function(n_iter, start, eps, L, x, y) {
  draws <- matrix(NA, n_iter, 2, dimnames = list(NULL, c("a", "b")))
  current <- start; accepted <- 0
  for (t in 1:n_iter) {
    q <- current; p <- rnorm(2)                          # a random push (momentum)
    H0 <- -log_post(q, x, y) + sum(p^2) / 2             # total energy at the start
    p <- p + eps / 2 * grad_log_post(q, x, y)            # leapfrog: half step of momentum,
    for (l in 1:L) {
      q <- q + eps * p                                   #   full step of position,
      if (l < L) p <- p + eps * grad_log_post(q, x, y)   #   full step of momentum,
    }
    p <- p + eps / 2 * grad_log_post(q, x, y)            #   final half step of momentum
    H1 <- -log_post(q, x, y) + sum(p^2) / 2             # total energy at the end
    if (log(runif(1)) < H0 - H1) { current <- q; accepted <- accepted + 1 }
    draws[t, ] <- current
  }
  attr(draws, "accept") <- accepted / n_iter
  draws
}
set.seed(3)
t_hmc  <- system.time(chain_h <- hmc(1000, start = c(0, 0), eps = 0.01, L = 20, x = x_c, y = y))
t_metr <- system.time(chain_m <- metropolis(5000, start = c(1, 1), step = 0.03, x = x_c, y = y))
tibble(sampler = c("Metropolis, 5,000 steps", "HMC, 1,000 steps"),
       accept  = c(attr(chain_m, "accept"), attr(chain_h, "accept")),
       ess_b   = c(ess_bulk(chain_m[-(1:500), "b"]), ess_bulk(chain_h[-(1:100), "b"])),
       seconds = c(t_metr[["elapsed"]], t_hmc[["elapsed"]])) |>
  mutate(ess_per_second = ess_b / seconds)

# A tibble: 2 x 5
  sampler                 accept ess_b seconds ess_per_second
  <chr>                    <dbl> <dbl>   <dbl>          <dbl>
1 Metropolis, 5,000 steps  0.486  467.   0.807           579.
2 HMC, 1,000 steps         0.989 1850.   3.93            471.

One thousand HMC iterations, each with 20 leapfrog steps, produce an effective sample size about twice the number of draws kept: successive HMC draws are essentially independent (an ESS larger than the number of draws means successive draws are slightly negatively correlated, which makes averages even more precise). One warning from experience: started at $(1,1)$ instead of $(0,0)$, this same sampler crawls for over 800 iterations before reaching the posterior, because its fixed step size and unit "mass" suit the posterior's centre, not the steep landscape far from it. Stan avoids this by adapting the step size and a mass matrix (the scale of each parameter) during warmup, which is one reason its warmup draws must be discarded. Metropolis needed five times as many iterations to produce about a quarter of the effective draws. In two dimensions, per second of computing, the two samplers come out roughly level (the timings vary from run to run and machine to machine), because each HMC iteration costs 21 gradient evaluations. In hundreds of dimensions the gap becomes enormous, because Metropolis's efficiency falls roughly in proportion to the number of parameters while HMC's falls much more slowly (Neal, 2011). That is why every modern general-purpose Bayesian engine uses gradients.

Lecture 10 — Trusting the result

5.5 Many chains: R-hat and effective sample size

A single chain can look perfectly healthy while exploring only part of the posterior. The standard defence is to run several chains (Stan's default is four) from different random starting points and check that they agree. If the chains have all converged to the same distribution, the variation between chains should be no larger than the variation within each chain. The potential scale reduction factor, $\widehat{R}$ ("R-hat"; Gelman and Rubin, 1992), compares the two.

With the centred parameterisation and a step size near 0.03, the four chains find each other within a few hundred steps and $\widehat{R}$ drops below 1.01 within a few thousand. Now switch to raw. The chains creep along the ridge from their separate starting points and, whatever step size you choose, $\widehat{R}$ stays far above 1: in Exercise 5, four raw chains still have $\widehat{R}\approx 1.7$ after 20,000 steps each. Any one of those chains, examined alone, might have looked like a slowly drifting but plausible trace. Only the comparison between chains reveals that none of them has explored the posterior.

Now fit the centred model with brms and look at its diagnostics the way you should for every model you fit.

d_age <- tibble(recid = y, x_c = x_c)
m_age <- brm(recid ~ 1 + x_c, data = d_age, family = bernoulli(),
             prior = c(prior(normal(0, 1.5), class = Intercept), prior(normal(0, 1), class = b)),
             seed = 5, refresh = 0, file = "fits/ch5_age")
draws_age <- as_draws_array(m_age, variable = c("b_Intercept", "b_x_c"))
summarise_draws(draws_age, "mean", "sd", "rhat", "ess_bulk", "ess_tail", "mcse_mean")
nuts_params(m_age) |> filter(Parameter == "divergent__") |> summarise(divergences = sum(Value))

# A tibble: 2 x 7
  variable      mean     sd  rhat ess_bulk ess_tail mcse_mean
  <chr>        <dbl>  <dbl> <dbl>    <dbl>    <dbl>     <dbl>
1 b_Intercept -0.208 0.0263  1.00    3290.    2575.  0.000458
2 b_x_c       -0.343 0.0236  1.00    2770.    2441.  0.000450
  divergences
1           0

Read the line for the slope, b_x_c. The posterior mean is −0.34, matching our hand-written samplers. $\widehat{R}$ is 1.00, the bulk and tail ESS are in the thousands, there are no divergences, and the Monte Carlo standard error of the mean (about 0.00045) is fifty times smaller than the posterior standard deviation. We can report the slope to two decimal places with complete confidence that the computation is not the limiting factor. Pictures help too: trace plots should look like overlapping "hairy caterpillars", and rank plots (histograms of the rank of each draw among all draws, separately by chain) should be flat for every chain.

dd <- as_draws_df(m_age) |> select(b_x_c, .chain, .iteration) |>
  mutate(chain = factor(.chain), rank = rank(b_x_c))
chain_cols <- c(harbour[["sky"]], harbour[["amber"]], harbour[["rose"]], harbour[["moss"]])
p1 <- ggplot(dd, aes(.iteration, b_x_c, colour = chain)) + geom_line(linewidth = 0.25, alpha = 0.8) +
  labs(x = "iteration (after warmup)", y = "slope b", colour = "chain")
p2 <- ggplot(dd, aes(rank, fill = chain)) + geom_histogram(bins = 20, colour = NA) +
  facet_wrap(~ chain, nrow = 2) + labs(x = "rank of draw among all chains", y = NULL) +
  theme(legend.position = "none", axis.text = element_blank())
patchwork::wrap_plots(p1, p2, widths = c(1.4, 1)) &
  scale_colour_manual(values = chain_cols) & scale_fill_manual(values = chain_cols)

5.6 When the sampler struggles, look at the model

The COMPAS file that ProPublica published contains several columns twice: the decile score appears as decile_score and again, identically, as decile_score.1 (the same is true of the count of prior arrests). This is common in real data assembled by joining tables. Suppose an analyst, not noticing, includes both versions in a model. To keep things fast, we aggregate the data by score, which gives a binomial model with ten rows (Exercise 3 of Chapter 4 showed that aggregation does not change the likelihood).

dup <- compas |>
  group_by(d1 = decile_score, d2 = decile_score.1) |>
  summarise(k = sum(two_year_recid), n = n(), .groups = "drop")
all(dup$d1 == dup$d2)                                       # the two columns are identical
m_dup_flat <- brm(k | trials(n) ~ 1 + d1 + d2, data = dup, family = binomial(),
                  seed = 5, refresh = 0, file = "fits/ch5_dup_flat")   # brms default: flat priors on b
summarise_draws(as_draws_array(m_dup_flat, variable = c("b_d1", "b_d2")),
                "mean", "sd", "rhat", "ess_bulk")
nuts_params(m_dup_flat) |> filter(Parameter == "treedepth__") |>
  summarise(iterations = n(), hit_max_treedepth = sum(Value >= 10))

[1] TRUE
# A tibble: 2 x 5
  variable  mean    sd  rhat ess_bulk
  <chr>    <dbl> <dbl> <dbl>    <dbl>
1 b_d1     -48.2  104.  2.37     4.96
2 b_d2      48.5  104.  2.37     4.96
  iterations hit_max_treedepth
1       4000              3233

The sampler complains loudly. $\widehat{R}$ is far above 1, the effective sample sizes are tiny, and most iterations hit the maximum tree depth (Stan prints a warning saying so): its trajectories kept running without ever turning back. The coefficients' posterior standard deviations are absurd. Nothing is wrong with Stan. The model has a problem: with two identical columns, only the sum $b_1+b_2$ affects the predictions, so any pair with the same sum fits the data equally well. With brms's default flat priors on the coefficients, nothing pins down the individual values, and the posterior is an infinitely long, flat ridge. The chains wander along it without end.

m_dup_prior <- update(m_dup_flat, prior = prior(normal(0, 1), class = b),
                      seed = 5, refresh = 0, file = "fits/ch5_dup_prior")
summarise_draws(as_draws_array(m_dup_prior, variable = c("b_d1", "b_d2")), "mean", "sd", "rhat", "ess_bulk")
sums <- as_draws_df(m_dup_prior) |> mutate(total = b_d1 + b_d2)
c(mean_sum = mean(sums$total), sd_sum = sd(sums$total))

bind_rows(`flat priors` = as_draws_df(m_dup_flat), `Normal(0, 1) priors` = as_draws_df(m_dup_prior),
          .id = "model") |>
  slice_head(n = 250, by = c(model, .chain)) |>                   # plot 250 draws per chain
  mutate(model = factor(model, levels = c("flat priors", "Normal(0, 1) priors"))) |>
  ggplot(aes(b_d1, b_d2, colour = factor(.chain))) +
  geom_point(alpha = 0.3, size = 0.6) +
  facet_wrap(~ model, scales = "free") +
  scale_colour_manual(values = c(harbour[["sky"]], harbour[["amber"]], harbour[["rose"]], harbour[["moss"]])) +
  labs(x = expression(b[1]), y = expression(b[2]), colour = "chain") + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1)))

# A tibble: 2 x 5
  variable  mean    sd  rhat ess_bulk
  <chr>    <dbl> <dbl> <dbl>    <dbl>
1 b_d1     0.159 0.727  1.00     874.
2 b_d2     0.122 0.727  1.00     874.
mean_sum   sd_sum 
 0.28034  0.01027 

With weakly informative Normal(0, 1) priors, the ridge becomes bounded: the prior gently prefers coefficients near zero, so the posterior is a bounded diagonal segment and the sampler converges. But look at what we learn: $b_1$ and $b_2$ each have huge uncertainty and are almost perfectly negatively correlated, while their sum is estimated precisely (it is the same as the coefficient on a single copy of the column). The right fix is not a prior; it is to remove the duplicate.

This example is extreme, but the pattern is general, and it has a name. Andrew Gelman's folk theorem of statistical computing says: when you have computational problems, often there's a problem with your model. Slow mixing, high $\widehat{R}$, maximum tree depth and divergences are frequently symptoms of predictors that are nearly collinear, parameters that the data cannot separate, priors that are too vague for the data, or a likelihood that does not suit the outcome. Before reaching for more iterations, ask what the sampler is telling you about the model.

5.7 Divergences

HMC has a built-in alarm that Metropolis lacks. As you saw in the trajectory widget, when the leapfrog step size is too large for the curvature of the posterior, the simulated energy explodes. Stan watches the energy along every trajectory (Stan Development Team, 2024), and when the error passes a threshold it stops the trajectory and records a divergent transition. Divergences matter because they reveal regions of the posterior that the sampler cannot enter: typically narrow, highly curved regions where the step size tuned for the rest of the posterior is far too large. The draws then systematically under-represent those regions, and the posterior summaries are biased, not just noisy.

A few divergences out of thousands of iterations, scattered randomly, can sometimes be removed by asking Stan for a smaller step size (control = list(adapt_delta = 0.95) or 0.99). Divergences that persist, or that cluster in one region of a pairs plot, need a better parameterisation. The classic case is the "funnel" in hierarchical models, where a group-level standard deviation near zero squeezes all the group effects into a narrow neck. Chapter 6 meets it on real data and shows the fix, the non-centred parameterisation, which is the same idea as centring a predictor here: change the coordinates so the posterior becomes a round hill.

5.8 Can the model recover the truth?

All the diagnostics so far check whether the sampler explored the posterior of the model you wrote. They cannot tell you whether that posterior would lead you to the right answer if the model were true. For that there is a simple, powerful habit: fake-data simulation. Choose parameter values, simulate a dataset from the model with those values (keeping the real predictors), fit the model, and check that the posterior covers the values you chose. Repeat a few times.

truth <- c(a = -0.2, b = -0.35)
recovery <- map_dfr(1:20, function(r) {
  set.seed(100 + r)
  fake <- tibble(x_c = x_c, recid = rbinom(length(x_c), 1, plogis(truth[["a"]] + truth[["b"]] * x_c)))
  fit <- update(m_age, newdata = fake, chains = 2, iter = 1000, seed = r, refresh = 0)
  as_tibble(fixef(fit), rownames = "parameter") |> mutate(rep = r)
}) |>
  mutate(parameter = recode(parameter, Intercept = "a (intercept)", x_c = "b (slope)"),
         true = if_else(parameter == "a (intercept)", truth[["a"]], truth[["b"]]),
         covered = Q2.5 <= true & true <= Q97.5)
recovery |> group_by(parameter) |> summarise(coverage = mean(covered), mean_error = mean(Estimate - true))
ggplot(recovery, aes(Estimate, rep)) +
  geom_vline(aes(xintercept = true), linetype = "dashed", colour = harbour[["muted"]]) +
  geom_errorbarh(aes(xmin = Q2.5, xmax = Q97.5, colour = covered), height = 0) +
  geom_point(aes(colour = covered), size = 1.2) +
  facet_wrap(~ parameter, scales = "free_x") +
  scale_x_continuous(n.breaks = 4) + theme(panel.spacing.x = unit(3, "lines")) +
  scale_colour_manual(values = c(`TRUE` = harbour[["sky"]], `FALSE` = harbour[["rose"]])) +
  labs(x = "posterior mean and 95% interval", y = "simulated dataset", colour = "covers truth")

# A tibble: 2 x 3
  parameter     coverage mean_error
  <chr>            <dbl>      <dbl>
1 a (intercept)     0.95    0.00321
2 b (slope)         0.9     0.00605

The intervals scatter around the true values, and the share that contains them is close to 95%, as it should be. Twenty repetitions is a small check (with 20 intervals, anywhere from 17 to 20 hits is consistent with 95% coverage), but it is enough to catch the commonest errors: a coding mistake, a variable on the wrong scale, a prior that silently dominates. For more complex models, run it before touching the real data, and especially for any parameter whose estimate your argument depends on.

5.9 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv:1701.02434.
2. Duane, S., Kennedy, A. D., Pendleton, B. J. and Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2), 216–222.
3. Gelfand, A. E. and Smith, A. F. M. (1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association, 85(410), 398–409.
4. Gelman, A., Roberts, G. O. and Gilks, W. R. (1996). Efficient Metropolis jumping rules. In J. M. Bernardo, J. O. Berger, A. P. Dawid and A. F. M. Smith (eds), Bayesian Statistics 5, 599–607. Oxford University Press.
5. Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472.
6. Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C. and Modrák, M. (2020). Bayesian workflow. arXiv:2011.01808.
7. Hoffman, M. D. and Gelman, A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, 1593–1623.
8. Kruschke, J. K. (2021). Bayesian analysis reporting guidelines. Nature Human Behaviour, 5, 1282–1291.
9. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087–1092.
10. Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K., Gelman, A. and Vehtari, A. (2025). Simulation-based calibration checking for Bayesian computation: the choice of test quantities shapes sensitivity. Bayesian Analysis, 20(2), 461–488.
11. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In S. Brooks et al. (eds), Handbook of Markov Chain Monte Carlo, 113–162. Chapman and Hall/CRC.
12. Roberts, G. O., Gelman, A. and Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Annals of Applied Probability, 7(1), 110–120.
13. Stan Development Team (2024). Stan Reference Manual. mc-stan.org.
14. Talts, S., Betancourt, M., Simpson, D., Vehtari, A. and Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. arXiv:1804.06788.
15. Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.