8.6 Practice

The practice problems for Chapter 8 ask you to run MCMC on simple models, diagnose convergence, and compare HMC with grid approximation.

Practice

Easy problems

8E1. Draw a DAG for the Metropolis algorithm. What are the nodes and what are the arrows?

8E2. Why does the Metropolis algorithm only need the ratio of probabilities, not the absolute values?

8E3. What is the difference between Metropolis, Metropolis-Hastings, and Gibbs sampling?

8E4. What does R-hat measure? What value indicates convergence?

Medium problems

8M1. Implement the Metropolis algorithm for a simple posterior and compare the samples to the exact posterior.

# 8M1: Metropolis algorithm for a simple posterior
# Model: y ~ Normal(mu, 1), prior: mu ~ Normal(0, 10)
# Exact posterior: mu ~ Normal(ybar, 1/n)

set.seed(42)
n <- 50
y <- rnorm(n, mean = 5, sd = 1)
ybar <- mean(y)

# Metropolis

n_iter <- 10000
mu <- 0
samples <- numeric(n_iter)

for (t in 1:n_iter) {
proposal <- mu + rnorm(1, sd = 0.3)
log_post <- function(m) -n/2 * (ybar - m)^2 - m^2 / 200
if (log(runif(1)) < log_post(proposal) - log_post(mu)) mu <- proposal
samples[t] <- mu
}

cat("Metropolis:\n")
cat(" Mean:", round(mean(samples), 3), "\n")
cat(" SD:", round(sd(samples), 3), "\n")

cat("\nExact posterior:\n")
cat(" Mean:", round(ybar, 3), "\n")
cat(" SD:", round(1/sqrt(n), 3), "\n")

# Compare histograms

hist(samples, breaks = 50, col = adjustcolor("steelblue", 0.5),
freq = FALSE, main = "8M1: Metropolis vs exact posterior",
xlab = "mu")
curve(dnorm(x, mean = ybar, sd = 1/sqrt(n)), add = TRUE, lwd = 3, col = "firebrick")
legend("topright", legend = c("Metropolis", "Exact"), col = c("steelblue", "firebrick"), lwd = 2)

Click Run to execute this code in your browser with WebR

8M2. Run the Metropolis algorithm with different proposal standard deviations (0.1, 0.5, 1.0). How does the proposal affect the chain?

# 8M2: Effect of proposal SD on Metropolis
set.seed(42)
n <- 50
y <- rnorm(n, mean = 5, sd = 1)
ybar <- mean(y)

log_post <- function(m) -n/2 * (ybar - m)^2 - m^2 / 200

run_metropolis <- function(prop_sd, n_iter = 3000) {
mu <- 0
samples <- numeric(n_iter)
accepts <- 0
for (t in 1:n_iter) {
proposal <- mu + rnorm(1, sd = prop_sd)
if (log(runif(1)) < log_post(proposal) - log_post(mu)) {
mu <- proposal
accepts <- accepts + 1
}
samples[t] <- mu
}
return(list(samples = samples, accept_rate = accepts / n_iter))
}

par(mfrow = c(3, 1))
for (sd in c(0.1, 0.5, 1.0)) {
result <- run_metropolis(sd)
plot(result$samples, type = "l", col = "steelblue",
     main = sprintf("Proposal SD = %.1f (accept rate = %.2f)", sd, result$accept_rate),
xlab = "Iteration", ylab = "mu")
}
par(mfrow = c(1, 1))

cat("Small proposal: high accept rate, slow exploration\n")
cat("Large proposal: low accept rate, fast but jumpy\n")
cat("Medium proposal: balanced accept rate and exploration\n")

Click Run to execute this code in your browser with WebR

Hard problems

8H1. Implement a simple HMC sampler for a 2D Gaussian posterior. Check convergence with R-hat.

# 8H1: Simple HMC for a 2D Gaussian
set.seed(42)

# Target: 2D standard normal

log_post <- function(theta) -0.5 * sum(theta^2)
grad_log_post <- function(theta) -theta

# HMC sampler

hmc_2d <- function(n_iter, step_size = 0.2, n_leapfrog = 15) {
theta <- c(0, 0) # starting value
samples <- matrix(0, nrow = n_iter, ncol = 2)

for (t in 1:n_iter) { # Random momentum
p <- rnorm(2)

  # Leapfrog
  theta_new <- theta
  p_new <- p
  for (i in 1:n_leapfrog) {
    p_new <- p_new + 0.5 * step_size * grad_log_post(theta_new)
    theta_new <- theta_new + step_size * p_new
    p_new <- p_new + 0.5 * step_size * grad_log_post(theta_new)
  }

  # Metropolis accept
  current_H <- -log_post(theta) + 0.5 * sum(p^2)
  proposed_H <- -log_post(theta_new) + 0.5 * sum(p_new^2)
  if (log(runif(1)) < current_H - proposed_H) theta <- theta_new
  samples[t, ] <- theta

}
return(samples)
}

# Run 4 chains

chains <- list()
for (c in 1:4) {
chains[[c]] <- hmc_2d(2000)
}

# Compute R-hat for each parameter

compute*rhat <- function(chain_list) {
m <- length(chain_list)
n <- nrow(chain_list[[1]])
mat <- sapply(chain_list, function(x) x[, 1])
W <- mean(apply(mat, 2, var))
B <- n * var(colMeans(mat))
var*hat <- (n - 1) / n * W + B / n
return(sqrt(var_hat / W))
}

cat("HMC 2D Gaussian:\n")
cat(" Parameter 1: mean =", round(mean(sapply(chains, function(x) mean(x[,1]))), 3),
" R-hat =", round(compute_rhat(chains), 4), "\n")

# Plot

plot(chains[[1]][,1], chains[[1]][,2], type = "l", col = "steelblue",
main = "8H1: HMC 2D Gaussian", xlab = "theta1", ylab = "theta2")
for (c in 2:4) {
lines(chains[[c]][,1], chains[[c]][,2], col = "steelblue")
}
points(0, 0, pch = 16, col = "firebrick", cex = 2)
legend("topright", legend = c("HMC samples", "True mean"),
col = c("steelblue", "firebrick"), pch = c(NA, 16), lwd = c(1, NA))

Click Run to execute this code in your browser with WebR

8H2. Compare the effective sample size of random-walk Metropolis and HMC for the same target distribution.

8H3. Simulate data from a model with a non-normal posterior (e.g., a mixture). Fit it with both quap and HMC. How do the results differ?

Summary

Maximum entropy