8.3 Easy HMC: map2stan

Hamiltonian Monte Carlo (HMC) is a powerful MCMC algorithm that uses the gradient of the log-posterior to propose distant states with high acceptance probability. The rethinking package exposes HMC through the map2stan (or ulam) interface.

8.3.1 How HMC works

HMC simulates a physics analogy: a particle slides on the surface of the negative log-posterior. The particle's position represents parameter values, and its momentum helps it explore the posterior efficiently.

Assign random momentum to the particle.
Simulate dynamics using the gradient of the log-posterior (the "force").
Propose the new position after a fixed number of steps.
Accept/Reject using a Metropolis criterion.

The key advantage: HMC uses gradient information to move far in a single step, unlike random-walk Metropolis which moves slowly.

# Demonstrate HMC conceptually: gradient-based exploration
# Compare random walk vs gradient-informed proposals

set.seed(42)

# Target: standard normal

log_post <- function(x) -0.5 * x^2
grad_log_post <- function(x) -x

# Random walk Metropolis

rw_samples <- numeric(5000)
x <- 0
for (t in 1:5000) {
proposal <- x + rnorm(1, sd = 0.5)
if (log(runif(1)) < log_post(proposal) - log_post(x)) x <- proposal
rw_samples[t] <- x
}

# HMC-like (leapfrog integration, simplified)

hmc_samples <- numeric(5000)
x <- 0
step_size <- 0.3
n_leapfrog <- 10

for (t in 1:5000) {

# Random momentum

p <- rnorm(1)

# Leapfrog steps

x*new <- x
p_new <- p
for (i in 1:n_leapfrog) {
p_new <- p_new + 0.5 * step*size * grad*log_post(x_new)
x_new <- x_new + step_size * p*new
p_new <- p_new + 0.5 * step_size * grad_log_post(x_new)
}

# Metropolis accept

current*H <- -log_post(x) + 0.5 * p^2
proposed*H <- -log_post(x_new) + 0.5 * p_new^2
if (log(runif(1)) < current_H - proposed_H) x <- x_new
hmc_samples[t] <- x
}

cat("Random walk Metropolis:\n")
cat(" Mean:", round(mean(rw_samples), 3), " SD:", round(sd(rw_samples), 3), "\n")
cat(" Effective sample size (approx):", round(5000 / (1 + 2 * cor(rw_samples[-1], rw_samples[-5000]))), 0), "\n")

cat("\nHMC:\n")
cat(" Mean:", round(mean(hmc_samples), 3), " SD:", round(sd(hmc_samples), 3), "\n")
cat(" Effective sample size (approx):", round(5000 / (1 + 2 * cor(hmc_samples[-1], hmc_samples[-5000]))), 0), "\n")

cat("\nHMC has lower autocorrelation = higher effective sample size.\n")

Click Run to execute this code in your browser with WebR

8.3.2 Using Stan (or ulam)

In practice, we use Stan (via ulam in rethinking or cmdstanr directly) which implements HMC with automatic differentiation:

# Simulate a simple model and fit with HMC (using lm as proxy)
# In rethinking: ulam() or map2stan() would call Stan
# Here we demonstrate the workflow

set.seed(42)
n <- 100
x <- runif(n, 0, 10)
y <- 2 + 1.5 * x + rnorm(n, sd = 2)

# In rethinking, this would be:

# m8.1 <- ulam(

# alist(

# y ~ dnorm(mu, sigma),

# mu <- a + b * x,

# a ~ dnorm(0, 10),

# b ~ dnorm(0, 10),

# sigma ~ dcauchy(0, 2)

# ),

# data = list(y = y, x = x),

# warmup = 1000,

# iter = 4000,

# chains = 4

# )

# Using lm as a proxy (gives same point estimates with flat priors)

fit <- lm(y ~ x)

# Simulate posterior samples (as HMC would produce)

library(MASS)
post <- mvrnorm(4000, mu = coef(fit), Sigma = vcov(fit))

cat("HMC-style posterior (simulated):\n")
cat(" a: mean =", round(mean(post[,1]), 3),
" 95% CI = [", round(quantile(post[,1], 0.025), 3), ",",
round(quantile(post[,1], 0.975), 3), "]\n")
cat(" b: mean =", round(mean(post[,2]), 3),
" 95% CI = [", round(quantile(post[,2], 0.025), 3), ",",
round(quantile(post[,2], 0.975), 3), "]\n")

# Trace plot (simulated)

par(mfrow = c(2, 1))
plot(post[,1], type = "l", col = "steelblue",
main = "Trace plot: intercept", xlab = "Iteration", ylab = "a")
plot(post[,2], type = "l", col = "firebrick",
main = "Trace plot: slope", xlab = "Iteration", ylab = "b")
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

Rethinking

HMC is not magic — it is physics. The beauty of HMC is that it turns sampling into a physics problem: a particle sliding on a surface. The gradient tells the particle which way is "down" (higher probability), and the momentum carries it through flat regions. This is why HMC is so much more efficient than random-walk methods: it uses information about the shape of the posterior, not just its values. But this also means HMC can fail when the surface is too rough or has sharp curvature — these are divergent transitions (Section 8.4).

Summary of Section 8.3

HMC uses the gradient of the log-posterior to propose distant, high-acceptance moves.
It simulates a physics analogy: position = parameters, momentum = exploration.
HMC has lower autocorrelation than random-walk Metropolis, giving higher effective sample size.
In practice, use Stan (via ulam or cmdstanr) which implements HMC with automatic differentiation.
HMC can fail on rough posteriors — these are divergent transitions.

Markov chain Monte Carlo

Care and feeding of your Markov chain