8.2 Markov chain Monte Carlo

Markov Chain Monte Carlo (MCMC) is a family of algorithms that generate samples from a target distribution by constructing a Markov chain whose stationary distribution is the posterior. The samples are correlated but, under regularity conditions, converge to the posterior.

8.2.1 Why MCMC?

Until now, we used quadratic approximation (quap) to estimate the posterior. This assumes the posterior is approximately multivariate normal. When this assumption fails — non-normal posteriors, discrete parameters, complex dependencies — we need MCMC.

MCMC works for any posterior, no matter its shape, because it does not approximate the posterior — it samples from it.

$p(\theta \mid D) = \frac{p(D \mid \theta) \cdot p(\theta)}{p(D)}$

MCMC avoids computing $p(D)$ (the normalizing constant) by using the ratio of unnormalized probabilities — just like King Markov.

# Compare quap (normal approximation) with MCMC samples
# for a posterior that is NOT normal

set.seed(42)

# Simulate from a skewed posterior (exponential-like)

# True posterior: theta ~ Gamma(2, 1)

n_samples <- 10000

# MCMC samples (here we just draw from the true distribution)

mcmc_samples <- rgamma(n_samples, shape = 2, rate = 1)

# Quadratic approximation: mean and SD of the samples

mean_q <- mean(mcmc_samples)
sd_q <- sd(mcmc_samples)
quap_samples <- rnorm(n_samples, mean = mean_q, sd = sd_q)

cat("Comparison: MCMC vs Quadratic approximation\n")
cat(" MCMC mean:", round(mean(mcmc_samples), 3), " SD:", round(sd(mcmc_samples), 3), "\n")
cat(" quap mean:", round(mean(quap_samples), 3), " SD:", round(sd(quap_samples), 3), "\n")
cat(" MCMC 95% CI: [", round(quantile(mcmc_samples, 0.025), 3), ",",
round(quantile(mcmc_samples, 0.975), 3), "]\n")
cat(" quap 95% CI: [", round(quantile(quap_samples, 0.025), 3), ",",
round(quantile(quap_samples, 0.975), 3), "]\n")

# Plot

par(mfrow = c(1, 2))
hist(mcmc_samples, breaks = 50, col = "steelblue", border = "white",
main = "MCMC (true posterior)", xlab = "theta", freq = FALSE)
hist(quap_samples, breaks = 50, col = "firebrick", border = "white",
main = "Quadratic approximation", xlab = "theta", freq = FALSE)
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

8.2.2 The Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm generalizes the Metropolis algorithm to allow arbitrary proposal distributions:

Propose: Draw $\theta' \sim q(\theta' \mid \theta)$ from a proposal distribution.
Accept ratio: $R = \frac{p(\theta' \mid D) \cdot q(\theta \mid \theta')}{p(\theta \mid D) \cdot q(\theta' \mid \theta)}$
Accept/Reject: Accept with probability $\min(1, R)$ .

# Metropolis-Hastings for a simple Bayesian model
# Model: y ~ Normal(mu, 1), prior: mu ~ Normal(0, 10)
# Posterior: mu ~ Normal(ybar, 1/n) (conjugate)

set.seed(42)
n <- 20
y <- rnorm(n, mean = 3, sd = 1) # true mu = 3
ybar <- mean(y)

# Metropolis-Hastings

n_iter <- 5000
mu <- 0 # starting value
samples <- numeric(n_iter)

for (t in 1:n_iter) {

# Propose from a normal centered at current value

proposal <- rnorm(1, mean = mu, sd = 0.5)

# Log posterior (up to constant): -n/2 _ (ybar - mu)^2 - mu^2 / (2 _ 100)

log*post <- function(m) {
-n/2 * (ybar - m)^2 - m^2 / (2 _ 100)
}

# Accept ratio

log_R <- log_post(proposal) - log_post(mu)

if (log(runif(1)) < log_R) {
mu <- proposal
}
samples[t] <- mu
}

cat("Metropolis-Hastings results:\n")
cat(" Posterior mean:", round(mean(samples), 3), "\n")
cat(" Posterior SD:", round(sd(samples), 3), "\n")
cat(" True posterior mean:", round(ybar, 3), "\n")

# Trace plot

plot(samples, type = "l", col = "steelblue",
main = "MCMC trace plot", xlab = "Iteration", ylab = "mu")

# Histogram

hist(samples, breaks = 50, col = "steelblue", border = "white",
main = "MCMC posterior samples", xlab = "mu", freq = FALSE)
curve(dnorm(x, mean = ybar, sd = 1/sqrt(n)), add = TRUE, lwd = 2, col = "firebrick")
legend("topright", legend = c("MCMC", "Exact"), col = c("steelblue", "firebrick"), lwd = 2)

Click Run to execute this code in your browser with WebR

8.2.3 Gibbs sampling

Gibbs sampling is a special case of Metropolis-Hastings where each parameter is updated conditional on all others. It is efficient when conditional distributions are easy to sample from.

Rethinking

MCMC is a tool, not a black box. The beauty of MCMC is that it works for any posterior. But this power comes with responsibility: you must check that the chain has converged. MCMC can get stuck, explore slowly, or miss regions of the posterior. Diagnostics (Section 8.4) help detect these problems. The key is to understand the mechanics — not to treat MCMC as a magic box that produces the answer.

Summary of Section 8.2

MCMC samples from the posterior without knowing the normalizing constant.
It works for any posterior shape, unlike quadratic approximation.
Metropolis-Hastings generalizes the Metropolis algorithm with arbitrary proposals.
Gibbs sampling updates each parameter conditional on the others.
MCMC is a tool, not a black box — always check convergence.

Good King Markov and His island kingdom

Easy HMC: map2stan