8.3 HMC fácil: map2stan
Hamiltonian Monte Carlo (HMC) es un potente algoritmo MCMC que usa el gradiente de la log-distribución a posteriori para proponer estados lejanos con alta probabilidad de aceptación. El paquete rethinking expone HMC a través de la interfaz map2stan (o ulam).
8.3.1 Cómo funciona HMC
HMC simula una analogía física: una partícula se desliza sobre la superficie del negativo de la log-distribución a posteriori. La posición de la partícula representa los valores de los parámetros, y su momento le ayuda a explorar la distribución a posteriori eficientemente.
- Asignar momento aleatorio a la partícula.
- Simular la dinámica usando el gradiente de la log-distribución a posteriori (la "fuerza").
- Proponer la nueva posición tras un número fijo de pasos.
- Aceptar/Rechazar usando un criterio de Metropolis.
La ventaja clave: HMC usa información del gradiente para moverse lejos en un solo paso, a diferencia del Metropolis de caminata aleatoria que se mueve lentamente.
# 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")8.3.2 Usando Stan (o ulam)
En la práctica, usamos Stan (vía ulam en rethinking o cmdstanr directamente) que implementa HMC con diferenciación automática:
# 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))HMC no es magia — es física. La belleza de HMC es que convierte el muestreo en un problema de física: una partícula deslizándose sobre una superficie. El gradiente le dice a la partícula hacia dónde está "cuesta abajo" (mayor probabilidad), y el momento la lleva a través de regiones planas. Por esto HMC es mucho más eficiente que los métodos de caminata aleatoria: usa información sobre la forma de la distribución a posteriori, no solo sus valores. Pero esto también significa que HMC puede fallar cuando la superficie es demasiado rugosa o tiene curvatura pronunciada — estas son transiciones divergentes (Sección 8.4).
Resumen de la Sección 8.3
- HMC usa el gradiente de la log-distribución a posteriori para proponer movimientos lejanos y de alta aceptación.
- Simula una analogía física: posición = parámetros, momento = exploración.
- HMC tiene menor autocorrelación que el Metropolis de caminata aleatoria, dando mayor tamaño efectivo de muestra.
- En la práctica, use Stan (vía
ulamocmdstanr) que implementa HMC con diferenciación automática. - HMC puede fallar en distribuciones a posteriori rugosas — estas son transiciones divergentes.