Introduction to Hamiltonian Monte Carlo

What is HMC?

Hamiltonian Monte Carlo (HMC) is a powerful MCMC method that suppresses random walk behavior by introducing an auxiliary momentum variable and using gradient information to guide sampling.

Key Problems HMC Solves:

Slow mixing in standard MCMC
High autocorrelation in Gibbs samplers
Difficult tuning of proposal distributions in Metropolis-Hastings
Inefficient exploration of high-dimensional posteriors

Core Intuition: Physics Analogy

Think of sampling from the posterior as simulating a frictionless puck sliding over a landscape:

Parameters \(\theta\) = position of the puck
Posterior log-density \(L(\theta)\) = negative potential energy
Auxiliary momentum r = momentum of the puck
Total energy = Potential + Kinetic = \(-L(\theta) + r^Tr/2\)

The Hamiltonian dynamics cause the puck to slide along equal-energy contours, efficiently exploring the posterior surface.

HMC Conceptual Diagram

Mathematical Foundation

The Joint Distribution

We augment the parameter space with an auxiliary momentum variable:

\[p(\theta, r | y) \propto p(\theta | y) \times \mathcal{N}(r | 0, I) \propto \exp\left(L(\theta) - \frac{r^T r}{2}\right)\]

Components:

\(p(\theta | y)\) = Target posterior distribution
\(L(\theta) = \log p(\theta | y)\) = Log-posterior
\(\mathcal{N}(r | 0, I)\) = Momentum distribution (independent)
\(-\frac{r^T r}{2}\) = Kinetic energy term

The Hamiltonian:

\[H(\theta, r) = -L(\theta) + \frac{r^T r}{2} = U(\theta) + K(r)\]

Potential Energy: \(U(\theta) = -L(\theta) = -\log p(\theta | y)\)
Kinetic Energy: \(K(r) = \frac{r^T r}{2}\)

Derivation of the HMC Joint Distribution

1. Understanding the Components

Let’s break down the equation step by step:

\[p(\theta, r | y) \propto p(\theta | y) \times \mathcal{N}(r | 0, I) \propto \exp\left(L(\theta) - \frac{r^T r}{2}\right)\]

Component 1: The Target Posterior \(p(\theta | y)\)

By Bayes’ theorem:

\[p(\theta | y) = \frac{p(y | \theta) \cdot p(\theta)}{p(y)} \propto p(y | \theta) \cdot p(\theta)\]

Where: - \(p(y | \theta)\) = Likelihood of data given parameters - \(p(\theta)\) = Prior distribution for parameters - \(p(y)\) = Marginal likelihood (normalizing constant)

Component 2: The Log-Posterior \(L(\theta)\)

We define:

\[L(\theta) = \log p(\theta | y)\]

This is the log-posterior (up to an additive constant). Since:

\[p(\theta | y) \propto p(y | \theta) \cdot p(\theta)\]

We can write:

\[L(\theta) = \log p(y | \theta) + \log p(\theta) + \text{constant}\]

Component 3: The Momentum Distribution \(\mathcal{N}(r | 0, I)\)

The momentum \(r\) is drawn from a multivariate standard normal:

\[\mathcal{N}(r | 0, I) = \frac{1}{(2\pi)^{d/2}} \exp\left(-\frac{r^T r}{2}\right)\]

Where \(d\) is the dimension of \(\theta\).

2. The Full Derivation

Step 1: Start with the augmented joint distribution

\[p(\theta, r | y) = p(\theta | y) \cdot p(r)\]

This assumes independence between \(\theta\) and \(r\) given \(y\). The momentum \(r\) is completely independent of the data \(y\).

Step 2: Substitute the momentum distribution

\[p(\theta, r | y) = p(\theta | y) \cdot \frac{1}{(2\pi)^{d/2}} \exp\left(-\frac{r^T r}{2}\right)\]

Step 3: Express in terms of log-posterior

Since \(p(\theta | y) \propto \exp(L(\theta))\), we can write:

\[p(\theta, r | y) \propto \exp(L(\theta)) \cdot \exp\left(-\frac{r^T r}{2}\right)\]

Step 4: Combine exponentials

\[p(\theta, r | y) \propto \exp\left(L(\theta) - \frac{r^T r}{2}\right)\]

Why Energy = Negative Log-Posterior

The Direct Mathematical Connection

This is a crucial concept that often causes confusion. Let me explain it step by step.

Step 1: Start with the Joint Distribution

From Bayes’ theorem and our augmentation:

\[p(\theta, r | y) \propto p(\theta | y) \times \mathcal{N}(r | 0, I)\]

Step 2: Express in Exponential Form

Any probability distribution can be written as an exponential:

\[p(\theta | y) \propto \exp(\log p(\theta | y))\]

So:

\[p(\theta, r | y) \propto \exp(\log p(\theta | y)) \times \exp\left(-\frac{r^T r}{2}\right)\]

Step 3: Combine Exponentials

\[p(\theta, r | y) \propto \exp\left(\log p(\theta | y) - \frac{r^T r}{2}\right)\]

Step 4: Define the Hamiltonian

We define:

\[H(\theta, r) = -\log p(\theta | y) + \frac{r^T r}{2}\]

This makes:

\[p(\theta, r | y) \propto \exp(-H(\theta, r))\]

The Physical Intuition

Think of the posterior as a landscape:

High Posterior = Low Potential Energy (Valleys)

Posterior: Large value (e.g., 0.4)
Log-posterior: Less negative (e.g., -0.9)
Energy: Small (e.g., 0.9)

Low Posterior = High Potential Energy (Mountains)

Posterior: Small value (e.g., 0.01)
Log-posterior: More negative (e.g., -4.6)
Energy: Large (e.g., 4.6)

# Simple 1D example with explicit numbers
theta <- seq(-3, 3, length.out = 100)

# Define posterior: Normal(0, 1)
log_posterior <- dnorm(theta, 0, 1, log = TRUE)
posterior <- exp(log_posterior)

# Energy (using our definition)
# H(θ) = -log p(θ) + r²/2
# At r = 0: H(θ) = -log p(θ)
energy_at_r0 <- -log_posterior

# At r = 1: H(θ) = -log p(θ) + 0.5
energy_at_r1 <- -log_posterior + 0.5

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# 1. Posterior
plot(theta, posterior, type = "l", lwd = 2,
     xlab = expression(theta), ylab = "p(θ | y)",
     main = "Posterior Distribution",
     col = "blue")
abline(v = 0, col = "red", lty = 2)
text(0, 0.35, "Highest density\n= Lowest energy", col = "red", cex = 0.8)

# 2. Log-Posterior
plot(theta, log_posterior, type = "l", lwd = 2,
     xlab = expression(theta), ylab = "log p(θ | y)",
     main = "Log-Posterior",
     col = "blue")
abline(v = 0, col = "red", lty = 2)
text(0, -0.5, "Highest log-posterior", col = "red", cex = 0.8)

# 3. Energy (at r = 0)
plot(theta, energy_at_r0, type = "l", lwd = 2,
     xlab = expression(theta), ylab = "H(θ, r=0)",
     main = "Energy H(θ, r=0) = -log p(θ|y)",
     col = "darkgreen")
abline(v = 0, col = "red", lty = 2)
text(0, 0.5, "Lowest energy = Highest density", col = "red", cex = 0.8)

# 4. Energy at different momentum values
colors <- c("blue", "green", "red", "purple")
r_values <- c(0, 0.5, 1, 1.5)

plot(theta, energy_at_r0, type = "l", lwd = 2,
     xlab = expression(theta), ylab = "H(θ, r)",
     main = "Energy for Different r Values",
     col = colors[1], ylim = c(0, 5))

for (i in 2:length(r_values)) {
    lines(theta, -log_posterior + r_values[i]^2/2, 
          col = colors[i], lwd = 2)
}
legend("topright", 
       paste("r =", r_values), 
       col = colors, lty = 1, lwd = 2, cex = 0.8)

Energy vs Log-Posterior Example

cat("\n=== Numeric Example ===\n")

## 
## === Numeric Example ===

cat("At θ = 0 (the mode):\n")

## At θ = 0 (the mode):

cat("  log p(θ | y) =", round(dnorm(0, 0, 1, log = TRUE), 3), "\n")

##   log p(θ | y) = -0.919

cat("  H(θ, r=0) =", round(-dnorm(0, 0, 1, log = TRUE), 3), "\n")

##   H(θ, r=0) = 0.919

cat("  H(θ, r=1) =", round(-dnorm(0, 0, 1, log = TRUE) + 0.5, 3), "\n\n")

##   H(θ, r=1) = 1.419

cat("At θ = 2 (low density):\n")

## At θ = 2 (low density):

cat("  log p(θ | y) =", round(dnorm(2, 0, 1, log = TRUE), 3), "\n")

##   log p(θ | y) = -2.919

cat("  H(θ, r=0) =", round(-dnorm(2, 0, 1, log = TRUE), 3), "\n\n")

##   H(θ, r=0) = 2.919

cat("Notice: When posterior density is high, energy is low!\n")

## Notice: When posterior density is high, energy is low!

cat("        When posterior density is low, energy is high!\n")

##         When posterior density is low, energy is high!

The Complete Picture

Mathematical Relationship

For any point (θ, r):

Joint distribution: \(p(\theta, r | y)\)
Hamiltonian (Energy): \(H(\theta, r) = -\log p(\theta, r | y)\)
Joint distribution in terms of energy: \(p(\theta, r | y) = \exp(-H(\theta, r))\)

The Components

Component	Formula	Physical Meaning	Statistical Meaning
Potential Energy	\(U(\theta) = -\log p(\theta \| y)\)	Height in landscape	Negative log-posterior
Kinetic Energy	\(K(r) = \frac{r^T r}{2}\)	Motion energy	Squared momentum/2
Total Energy	\(H(\theta, r) = U(\theta) + K(r)\)	Total energy	Negative log-joint

Why Is This Important?

1. Conservation of Energy

In Hamiltonian dynamics, total energy is conserved:

\[\frac{dH}{dt} = 0\]

This means: - The trajectory moves along contours of equal energy - Equal energy = Equal joint probability - So we explore the posterior efficiently!

2. The Acceptance Probability

The Metropolis acceptance uses:

\[\alpha = \min\left(1, \frac{\exp(-H(\theta^*, r^*))}{\exp(-H(\theta, r))}\right)\]

Which simplifies to:

\[\alpha = \min(1, \exp(-(H(\theta^*, r^*) - H(\theta, r))))\]

If energy is conserved perfectly, \(\alpha = 1\)!

3. The Gradient Connection

The leapfrog uses:

\[\frac{\partial H}{\partial \theta} = -\frac{\partial}{\partial \theta} \log p(\theta | y) = -\nabla L(\theta)\]

So the gradient of the log-posterior drives the dynamics!

# Create a nice visualization
theta <- seq(-3, 3, length.out = 200)
posterior <- dnorm(theta, 0, 1)
log_posterior <- dnorm(theta, 0, 1, log = TRUE)
energy <- -log_posterior

par(mfrow = c(1, 3), mar = c(4, 4, 3, 2))

# 1. Posterior (as a hill)
plot(theta, posterior, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "p(θ | y)",
     main = "Posterior (Hill)",
     col = "blue")
polygon(c(theta, rev(theta)), c(rep(0, length(theta)), rev(posterior)),
        col = rgb(0, 0, 1, alpha = 0.3))
text(0, 0.3, "High\nDensity", col = "darkblue", cex = 1.2)
text(2, 0.05, "Low\nDensity", col = "darkblue", cex = 1.2)

# 2. Log-posterior
plot(theta, log_posterior, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "log p(θ | y)",
     main = "Log-Posterior",
     col = "blue")
abline(h = -2, col = "red", lty = 2)
text(0, -0.5, "High", col = "red", cex = 1.2)
text(2, -3, "Low", col = "red", cex = 1.2)

# 3. Energy (as a valley)
plot(theta, energy, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "H(θ, r=0)",
     main = "Energy (Valley)",
     col = "darkgreen")
polygon(c(theta, rev(theta)), c(rep(0, length(theta)), rev(energy)),
        col = rgb(0, 0.5, 0, alpha = 0.3))
text(0, 0.5, "Low\nEnergy", col = "darkgreen", cex = 1.2)
text(2, 2.5, "High\nEnergy", col = "darkgreen", cex = 1.2)

Posterior as Energy Landscape

cat("\n=== Analogy ===\n")

## 
## === Analogy ===

cat("Posterior is like a mountain (high = good).\n")

## Posterior is like a mountain (high = good).

cat("Energy is like a valley (low = good).\n")

## Energy is like a valley (low = good).

cat("They are inverses of each other!\n")

## They are inverses of each other!

The Boltzmann Distribution

What is the Boltzmann Distribution?

The Boltzmann distribution (also called the Gibbs distribution) is a fundamental concept in statistical physics that describes how energy is distributed among particles in a system at thermal equilibrium.

The Mathematical Form

The Boltzmann distribution states that the probability of a system being in a particular state is:

\[P(\text{state}) \propto \exp\left(-\frac{E(\text{state})}{kT}\right)\]

Where: - E(state) = Energy of the state - k = Boltzmann constant (≈ 1.38 × 10⁻²³ J/K) - T = Absolute temperature (in Kelvin)

Or more simply, when kT = 1:

\[P(\text{state}) \propto \exp(-E(\text{state}))\]

Visualizing the Boltzmann Distribution

# Create a clear visualization of the Boltzmann distribution
energy_levels <- seq(0, 5, length.out = 100)

# Different temperatures
temperature_levels <- c(0.5, 1, 2, 5)  # In units where k = 1
colors <- c("red", "blue", "green", "purple")

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# 1. Basic Boltzmann distribution at T = 1
probability <- exp(-energy_levels)
plot(energy_levels, probability, type = "l", lwd = 3,
     xlab = "Energy E", ylab = "P(state)",
     main = "Boltzmann Distribution (T = 1)",
     col = "blue")
abline(v = 0, col = "red", lty = 2)
text(0, 0.5, "Highest probability", col = "red", cex = 0.9)
text(3, 0.1, "Low probability", col = "red", cex = 0.9)

# 2. Effect of temperature
plot(energy_levels, exp(-energy_levels/0.5), type = "l", 
     lwd = 3, col = colors[1],
     xlab = "Energy E", ylab = "P(state)",
     main = "Effect of Temperature",
     ylim = c(0, 1))

for (i in 1:length(temperature_levels)) {
    lines(energy_levels, exp(-energy_levels/temperature_levels[i]), 
          col = colors[i], lwd = 2)
}
legend("topright", 
       paste("T =", temperature_levels), 
       col = colors, lty = 1, lwd = 2, cex = 0.8)

# 3. Energy levels (discrete states)
states <- 1:10
energies <- c(0, 0.5, 0.5, 1, 1.5, 2, 2.5, 3, 4, 5)
probabilities <- exp(-energies) / sum(exp(-energies))

barplot(probabilities, names.arg = states,
        xlab = "State", ylab = "Probability",
        main = "Discrete Boltzmann Distribution",
        col = rainbow(10))
abline(h = 0, col = "gray")

# 4. Compare with our HMC energy
theta <- seq(-3, 3, length.out = 100)
posterior <- dnorm(theta, 0, 1)
log_posterior <- dnorm(theta, 0, 1, log = TRUE)
hmc_energy <- -log_posterior

plot(theta, posterior, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "p(θ | y)",
     main = "HMC Connection",
     col = "blue")
# Add Boltzmann distribution with energy = -log(posterior)
points(theta, exp(-hmc_energy), col = "red", pch = 16, cex = 0.3)

legend("topright", 
       c("Posterior", "exp(-H(θ)) = exp(-log p(θ|y))"),
       col = c("blue", "red"), 
       lty = c(1, NA), pch = c(NA, 16),
       cex = 0.8)

Visualizing the Boltzmann Distribution

cat("=== Key Points ===\n")

## === Key Points ===

cat("1. Lower energy states have HIGHER probability\n")

## 1. Lower energy states have HIGHER probability

cat("2. Higher energy states have LOWER probability\n")

## 2. Higher energy states have LOWER probability

cat("3. Temperature controls how spread out the distribution is\n")

## 3. Temperature controls how spread out the distribution is

cat("4. In HMC, we set kT = 1 for convenience\n")

## 4. In HMC, we set kT = 1 for convenience

Simple Examples of the Boltzmann Distribution

Example 1: Two-State System

# Two-state example
energy_diff <- seq(0, 5, length.out = 100)
prob_state1 <- 1 / (1 + exp(-energy_diff))  # Probability of being in state 1
prob_state2 <- exp(-energy_diff) / (1 + exp(-energy_diff))  # Probability of state 2

par(mfrow = c(1, 2), mar = c(4, 4, 3, 2))

# 1. Probabilities as function of energy difference
plot(energy_diff, prob_state1, type = "l", lwd = 3,
     xlab = "Energy Difference (E₂ - E₁)", ylab = "Probability",
     main = "Two-State Boltzmann Distribution",
     col = "blue", ylim = c(0, 1))
lines(energy_diff, prob_state2, col = "red", lwd = 3)
abline(h = 0.5, col = "gray", lty = 2)
abline(v = 0, col = "gray", lty = 2)

legend("topright", 
       c("P(state 1)", "P(state 2)"), 
       col = c("blue", "red"), lty = 1, lwd = 2, cex = 0.8)

# 2. Visualize the two states
barplot(c(0.7, 0.3), names.arg = c("State 1 (Low E)", "State 2 (High E)"),
        ylab = "Probability", main = "Example: Low Energy State is More Likely",
        col = c("blue", "red"))

Two-State Boltzmann Distribution

cat("When E₂ - E₁ > 0 (state 2 has higher energy):\n")

## When E₂ - E₁ > 0 (state 2 has higher energy):

cat("  P(state 1) > P(state 2)\n")

##   P(state 1) > P(state 2)

cat("  The lower energy state is MORE probable!\n\n")

##   The lower energy state is MORE probable!

cat("When E₂ - E₁ < 0 (state 1 has higher energy):\n")

## When E₂ - E₁ < 0 (state 1 has higher energy):

cat("  P(state 2) > P(state 1)\n")

##   P(state 2) > P(state 1)

cat("  The lower energy state is MORE probable!\n")

##   The lower energy state is MORE probable!

Example 2: Three-State System

# Three-state example
states <- c("A", "B", "C")
energies <- c(0, 1, 3)  # Different energy levels

# Compute probabilities at different temperatures
temperatures <- c(0.5, 1, 2)

par(mfrow = c(1, 3), mar = c(4, 4, 3, 2))

for (T in temperatures) {
    probs <- exp(-energies / T) / sum(exp(-energies / T))
    barplot(probs, names.arg = states,
            ylab = "Probability", 
            main = paste("T =", T),
            col = c("blue", "green", "red"),
            ylim = c(0, 1))
    # Add energy labels
    text(0.7, 0.95, paste("E =", energies[1]), cex = 0.8)
    text(1.9, 0.95, paste("E =", energies[2]), cex = 0.8)
    text(3.1, 0.95, paste("E =", energies[3]), cex = 0.8)
}

Three-State Boltzmann Distribution

State A: E = 0 (lowest energy)

State B: E = 1 (medium energy)

State C: E = 3 (highest energy)

At low temperature (T = 0.5):

State A is overwhelmingly most likely!

The system freezes into the lowest energy state.

At high temperature (T = 2):

States are more evenly distributed.

Thermal energy allows exploring higher energy states.

The Connection to HMC

How HMC Uses the Boltzmann Distribution

In HMC, we create an analogy between:

Statistical Physics → Bayesian Statistics
Energy of a state → Negative log-posterior
Temperature → Set to 1 (kT = 1)

This gives us:

# Show the complete connection
theta <- seq(-4, 4, length.out = 200)

# Create a complex posterior (mixture)
posterior <- 0.3 * dnorm(theta, -2, 0.5) + 0.4 * dnorm(theta, 0, 0.8) + 0.3 * dnorm(theta, 2, 0.6)
posterior <- posterior / sum(posterior) * diff(theta)[1]

# Energy = -log(posterior)
energy <- -log(posterior + 1e-10)

# Boltzmann probability with energy
boltzmann_prob <- exp(-energy)
boltzmann_prob <- boltzmann_prob / sum(boltzmann_prob) * diff(theta)[1]

par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))

# 1. Posterior
plot(theta, posterior, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "p(θ | y)",
     main = "Posterior Distribution",
     col = "blue")

# 2. Energy (negative log-posterior)
plot(theta, energy, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "H(θ) = -log p(θ|y)",
     main = "Energy Landscape",
     col = "darkgreen")

# 3. Boltzmann distribution with energy = -log(posterior)
plot(theta, boltzmann_prob, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "exp(-H(θ))",
     main = "Boltzmann Distribution",
     col = "red")

# 4. Compare posterior vs Boltzmann
plot(theta, posterior, type = "l", lwd = 3,
     xlab = expression(theta), ylab = "Density",
     main = "Posterior = Boltzmann Distribution",
     col = "blue")
lines(theta, boltzmann_prob, col = "red", lwd = 2, lty = 2)
legend("topright", 
       c("Posterior", "exp(-H(θ))"), 
       col = c("blue", "red"), 
       lty = c(1, 2), lwd = 2, cex = 0.8)

HMC-Boltzmann Connection

cat("=== The HMC Connection ===\n")

## === The HMC Connection ===

cat("In HMC, we define: H(θ, r) = -log p(θ | y) + r²/2\n\n")

## In HMC, we define: H(θ, r) = -log p(θ | y) + r²/2

cat("Then: p(θ, r | y) = exp(-H(θ, r)) / Z\n\n")

## Then: p(θ, r | y) = exp(-H(θ, r)) / Z

cat("This is EXACTLY the Boltzmann distribution with:\n")

## This is EXACTLY the Boltzmann distribution with:

cat("  - Energy = H(θ, r)\n")

##   - Energy = H(θ, r)

cat("  - Temperature set so that kT = 1\n")

##   - Temperature set so that kT = 1

cat("  - States = (θ, r) pairs\n\n")

##   - States = (θ, r) pairs

cat("So HMC is using the Boltzmann distribution to sample from the posterior!\n")

## So HMC is using the Boltzmann distribution to sample from the posterior!

Physical Interpretation

What Does Temperature Mean in HMC?

In statistical physics: - Low temperature → System stays in low-energy states - High temperature → System explores more states

In HMC: - We effectively set T = 1/k (so kT = 1) - This is a “natural” temperature that makes the posterior the target distribution - If we wanted different behavior, we could set different temperatures (this is called “tempered” or “annealed” MCMC)

The Role of the Momentum

The momentum adds kinetic energy:

\[H(\theta, r) = U(\theta) + K(r)\]

Where: - U(θ) = Potential energy (from the posterior) - K(r) = Kinetic energy (from the momentum)

This is exactly like a physical system where: - Particles move in a potential energy landscape - Temperature determines how much they move - Kinetic energy allows them to explore different regions

Why Is This So Powerful?

1. It Provides a Theoretical Foundation

The Boltzmann distribution is well-understood in physics. By mapping our problem to it, we can use all the tools of statistical mechanics!

2. It Explains Why HMC Works

The Boltzmann distribution naturally samples states according to their probability. Lower energy = more probable. This is exactly what we want for Bayesian inference!

3. It Gives Us Conservation Laws

In Hamiltonian dynamics, energy is conserved. This means: - We move along constant energy contours - We explore regions of equal probability efficiently - No random walk behavior!

4. It Connects to Temperature-Based Methods

Methods like: - Simulated annealing (decreasing temperature to find modes) - Parallel tempering (running chains at different temperatures) - Thermodynamic integration (computing normalizing constants)

All come naturally from the Boltzmann distribution framework.

Complete Example: Sampling from the Boltzmann Distribution

# Demonstrate sampling from the Boltzmann distribution
# This is essentially what HMC does!

# Define a potential energy (negative log-posterior)
potential <- function(theta) {
    # Two wells with a barrier between them
    return((theta^2 - 1)^2 + 0.5 * theta)
}

# Define gradient
grad_potential <- function(theta) {
    h <- 1e-6
    (potential(theta + h) - potential(theta - h)) / (2 * h)
}

# Simple Metropolis-Hastings sampler for Boltzmann distribution
sample_boltzmann <- function(potential, grad_potential, 
                            n_samples, epsilon, n_steps) {
    
    samples <- numeric(n_samples)
    theta <- 0
    accept <- 0
    
    for (i in 1:n_samples) {
        # Propose using gradient (simplified HMC-like move)
        theta_star <- theta + rnorm(1, 0, epsilon)
        
        # Accept/reject based on energy difference
        dE <- potential(theta_star) - potential(theta)
        
        if (log(runif(1)) < -dE) {
            theta <- theta_star
            accept <- accept + 1
        }
        
        samples[i] <- theta
    }
    
    return(list(samples = samples, acceptance = accept / n_samples))
}

# Run the sampler
set.seed(123)
theta_range <- seq(-2.5, 2.5, length.out = 200)
n_samples <- 5000

# Sample at different temperatures (kT values)
temperatures <- c(0.2, 1.0, 2.0)
results <- list()

par(mfrow = c(2, 3), mar = c(4, 4, 3, 2))

for (i in 1:length(temperatures)) {
    # Sample from Boltzmann distribution
    # We incorporate temperature by scaling the energy
    scaled_potential <- function(theta) potential(theta) / temperatures[i]
    scaled_grad <- function(theta) grad_potential(theta) / temperatures[i]
    
    result <- sample_boltzmann(scaled_potential, scaled_grad, 
                              n_samples, epsilon = 0.1, n_steps = 1)
    results[[i]] <- result
    
    # Plot histogram
    hist(result$samples, breaks = 30, prob = TRUE,
         xlab = expression(theta), 
         main = paste("T =", temperatures[i]),
         col = rgb(0, 0, 1, alpha = 0.5),
         xlim = c(-2.5, 2.5))
    
    # Add theoretical Boltzmann distribution
    energy_vals <- sapply(theta_range, potential) / temperatures[i]
    boltzmann_vals <- exp(-energy_vals)
    boltzmann_vals <- boltzmann_vals / sum(boltzmann_vals) / diff(theta_range)[1]
    lines(theta_range, boltzmann_vals, col = "red", lwd = 2)
    
    # Add potential
    pot_vals <- sapply(theta_range, potential)
    lines(theta_range, pot_vals / max(pot_vals) * 0.5, 
          col = "darkgreen", lwd = 2, lty = 2)
    
    cat("Temperature T =", temperatures[i], "\n")
    cat("  Acceptance rate:", round(result$acceptance, 3), "\n")
    cat("  Mean theta:", round(mean(result$samples), 3), "\n\n")
}

## Temperature T = 0.2 
##   Acceptance rate: 0.801 
##   Mean theta: -0.943

## Temperature T = 1 
##   Acceptance rate: 0.935 
##   Mean theta: -0.176

## Temperature T = 2 
##   Acceptance rate: 0.945 
##   Mean theta: -0.909

cat("=== Effect of Temperature ===\n")

## === Effect of Temperature ===

cat("T = 0.2: Low temperature - samples concentrated in low energy regions\n")

## T = 0.2: Low temperature - samples concentrated in low energy regions

cat("T = 1.0: Moderate temperature - explores more of the landscape\n")

## T = 1.0: Moderate temperature - explores more of the landscape

cat("T = 2.0: High temperature - samples spread out more\n")

## T = 2.0: High temperature - samples spread out more

cat("\nIn HMC, we typically use T = 1 (kT = 1)!\n")

## 
## In HMC, we typically use T = 1 (kT = 1)!

Sampling from the Boltzmann Distribution

Summary of the Boltzmann Distribution in HMC

The Key Equation to Remember:

\[p(\theta, r | y) \propto \exp(-H(\theta, r))\]

Where \(H(\theta, r) = -\log p(\theta | y) + \frac{r^T r}{2}\)

This is the Boltzmann distribution with: - Energy = H(θ, r) - Temperature = 1/k (so kT = 1) - States = (θ, r) pairs - Probability = The joint posterior we want to sample from

Key Insights:

Lower energy = Higher probability: The Boltzmann distribution naturally samples from high-probability regions
Temperature controls exploration: Higher temperature = more exploration
Energy is conserved: Hamiltonian dynamics preserve energy, leading to efficient exploration
Physical intuition: We can think of the posterior as an energy landscape

The beauty of HMC is that it simulates a physical system (Boltzmann distribution) to solve a statistical problem (Bayesian inference)!

The Leapfrog Integrator

Why Leapfrog?

The leapfrog integrator is a numerical method for approximating Hamiltonian dynamics with three key properties:

Volume preserving (symplectic)
Time-reversible
Energy conserving (approximately)

Leapfrog Algorithm

The leapfrog updates proceed in three steps:

Half-step momentum:
\[r(t + \varepsilon/2) = r(t) + \frac{\varepsilon}{2}\nabla L(\theta(t))\]
Full-step position:
\[\theta(t + \varepsilon) = \theta(t) + \varepsilon \cdot r(t + \varepsilon/2)\]
Half-step momentum:
\[r(t + \varepsilon) = r(t + \varepsilon/2) + \frac{\varepsilon}{2}\nabla L(\theta(t + \varepsilon))\]

#' Leapfrog Integrator Implementation
leapfrog <- function(grad_log_posterior, theta, r, epsilon) {
    # Half-step update for momentum
    r_half <- r + (epsilon / 2) * grad_log_posterior(theta)
    
    # Full-step update for position
    theta_new <- theta + epsilon * r_half
    
    # Half-step update for momentum using new position
    r_new <- r_half + (epsilon / 2) * grad_log_posterior(theta_new)
    
    return(list(theta = theta_new, r = r_new))
}

#' Leapfrog with Multiple Steps
leapfrog_multi <- function(grad_log_posterior, theta, r, epsilon, L) {
    trajectory <- matrix(NA, nrow = L + 1, ncol = length(theta))
    trajectory[1, ] <- theta
    
    for (i in 1:L) {
        result <- leapfrog(grad_log_posterior, theta, r, epsilon)
        theta <- result$theta
        r <- result$r
        trajectory[i + 1, ] <- theta
    }
    
    return(list(
        theta = theta,
        r = r,
        trajectory = trajectory
    ))
}

#' Demonstrate Leapfrog Properties
demonstrate_leapfrog <- function() {
    # Simple 1D example: Normal(0,1) target
    grad_log <- function(theta) -theta
    
    # Initial conditions
    theta0 <- 2
    r0 <- 0.5
    epsilon <- 0.3
    L <- 20
    
    # Run leapfrog
    result <- leapfrog_multi(grad_log, theta0, r0, epsilon, L)
    
    # Plot
    par(mfrow = c(2, 2), mar = c(4, 4, 2, 2))
    
    # Phase space trajectory
    plot(1, type = "n", xlim = c(-3, 3), ylim = c(-2, 2),
         xlab = expression(theta), ylab = "r",
         main = "Phase Space Trajectory")
    points(result$trajectory[, 1], rep(r0, L+1), 
           type = "o", col = "blue", pch = 16, cex = 0.8)
    points(result$trajectory[, 1], rep(r0, L+1), 
           type = "o", col = "red", pch = 16, cex = 0.8)
    
    # Evolution of theta
    plot(0:L, result$trajectory[, 1], type = "o",
         xlab = "Leapfrog Step", ylab = expression(theta),
         main = "Position Evolution", col = "blue")
    abline(h = 0, col = "red", lty = 2)
    
    # Energy conservation
    H <- function(theta, r) theta^2/2 + r^2/2
    energy <- sapply(1:(L+1), function(i) H(result$trajectory[i, 1], r0))
    plot(0:L, energy, type = "o",
         xlab = "Leapfrog Step", ylab = "Hamiltonian H",
         main = "Energy Conservation", col = "darkgreen")
    
    # Error in energy
    rel_error <- abs(energy - energy[1]) / abs(energy[1])
    plot(0:L, rel_error, type = "o",
         xlab = "Leapfrog Step", ylab = "Relative Error",
         main = "Energy Error", col = "red")
    abline(h = 0, col = "gray", lty = 2)
}

demonstrate_leapfrog()

Leapfrog Integration Demonstration

The HMC Algorithm

Complete HMC Algorithm

For iteration \(t = 1\) to \(T\):

Sample \(r^* \sim \mathcal{N}(0, I)\)
Set \((\theta_0, r_0) = (\theta^{(t-1)}, r^*)\)
For \(l = 1\) to \(L\): \[(\theta_l, r_l) = \text{Leapfrog}(\theta_{l-1}, r_{l-1})\]
Accept \((\theta_L, r_L)\) with probability: \[\min\left(1, \exp\left(L(\theta_L) - \frac{r_L^T r_L}{2} - L(\theta^{(t-1)}) + \frac{r^{*T} r^*}{2}\right)\right)\]
If accepted: \(\theta^{(t)} = \theta_L\)
Else: \(\theta^{(t)} = \theta^{(t-1)}\)

R Implementation

#' Hamiltonian Monte Carlo Sampler
#' 
#' @param log_posterior Function computing log of target posterior
#' @param grad_log_posterior Function computing gradient of log posterior
#' @param theta_init Initial parameter values
#' @param n_samples Number of samples to generate
#' @param L_steps Number of leapfrog steps per iteration
#' @param epsilon Step size
#' @param warmup Number of warmup iterations
#' @param verbose Print progress
#' @return List containing samples, acceptance rates, and diagnostics
hmc <- function(log_posterior, grad_log_posterior,
                theta_init, n_samples, L_steps, epsilon,
                warmup = 500, verbose = TRUE) {
    
    # Input validation
    if (!is.numeric(theta_init)) stop("theta_init must be numeric")
    d <- length(theta_init)
    if (d == 0) stop("theta_init must have at least one element")
    if (n_samples <= 0) stop("n_samples must be positive")
    if (L_steps < 1) stop("L_steps must be at least 1")
    if (epsilon <= 0) stop("epsilon must be positive")
    
    # Initialize storage
    total_iterations <- n_samples + warmup
    samples <- matrix(NA, nrow = n_samples, ncol = d)
    acceptance <- logical(total_iterations)
    theta <- theta_init
    n_accepted <- 0
    
    if (verbose) {
        cat("Starting HMC with", d, "parameters\n")
        cat("Step size:", epsilon, "| Leapfrog steps:", L_steps, "\n")
        cat("Total iterations:", total_iterations, "\n")
    }
    
    for (t in 1:total_iterations) {
        # Step 1: Sample momentum from standard normal
        r_star <- rnorm(d)
        r_current <- r_star
        theta_current <- theta
        
        # Step 2: Simulate Hamiltonian dynamics
        for (l in 1:L_steps) {
            result <- leapfrog(grad_log_posterior, theta, r_star, epsilon)
            theta <- result$theta
            r_star <- result$r
        }
        
        # Step 3: Metropolis acceptance
        # Compute log probabilities (joint density)
        current_log_prob <- log_posterior(theta_current) - sum(r_current^2) / 2
        proposed_log_prob <- log_posterior(theta) - sum(r_star^2) / 2
        
        # Log acceptance ratio
        log_accept_ratio <- proposed_log_prob - current_log_prob
        
        # Accept or reject
        if (is.finite(log_accept_ratio) && 
            log(runif(1)) < min(0, log_accept_ratio)) {
            # Accept proposal
            if (t > warmup) {
                samples[t - warmup, ] <- theta
                n_accepted <- n_accepted + 1
            }
            acceptance[t] <- TRUE
        } else {
            # Reject proposal - revert to previous state
            if (t > warmup) {
                samples[t - warmup, ] <- theta_current
            }
            theta <- theta_current
            acceptance[t] <- FALSE
        }
        
        # Progress reporting
        if (verbose && t %% 1000 == 0) {
            current_rate <- if (t > 0) sum(acceptance[1:t]) / t else 0
            cat(sprintf("Iteration %d/%d | Acceptance: %.3f\n", 
                       t, total_iterations, current_rate))
        }
    }
    
    # Compute final acceptance rate
    final_rate <- if (n_samples > 0) n_accepted / n_samples else 0
    
    if (verbose) {
        cat("\nHMC Complete\n")
        cat("Acceptance rate:", round(final_rate, 3), "\n")
        cat("Effective samples:", n_samples, "\n")
    }
    
    return(list(
        samples = samples,
        acceptance_rate = final_rate,
        acceptance = acceptance,
        theta_final = theta,
        n_accepted = n_accepted
    ))
}

Examples

Example 1: Sampling from Multivariate Normal

# Define target: Bivariate Normal with correlation
set.seed(42)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0.7, 0.7, 1), nrow = 2)
Sigma_inv <- solve(Sigma)

log_posterior_mvn <- function(theta) {
    if (length(theta) != 2) stop("theta must be length 2")
    diff <- theta - mu
    return(-0.5 * t(diff) %*% Sigma_inv %*% diff)
}

grad_log_posterior_mvn <- function(theta) {
    if (length(theta) != 2) stop("theta must be length 2")
    diff <- theta - mu
    return(-Sigma_inv %*% diff)
}

# Run HMC
theta_init <- c(2, 2)
n_samples <- 5000
L_steps <- 25
epsilon <- 0.15

cat("=== Running HMC for Bivariate Normal ===\n")

## === Running HMC for Bivariate Normal ===

hmc_result <- hmc(log_posterior_mvn, grad_log_posterior_mvn,
                  theta_init, n_samples, L_steps, epsilon,
                  warmup = 1000, verbose = TRUE)

## Starting HMC with 2 parameters
## Step size: 0.15 | Leapfrog steps: 25 
## Total iterations: 6000 
## Iteration 1000/6000 | Acceptance: 0.997
## Iteration 2000/6000 | Acceptance: 0.998
## Iteration 3000/6000 | Acceptance: 0.997
## Iteration 4000/6000 | Acceptance: 0.997
## Iteration 5000/6000 | Acceptance: 0.997
## Iteration 6000/6000 | Acceptance: 0.997
## 
## HMC Complete
## Acceptance rate: 0.997 
## Effective samples: 5000

# Extract samples
samples <- hmc_result$samples
acceptance_rate <- hmc_result$acceptance_rate

# Visualize results
par(mfrow = c(2, 3), mar = c(4, 4, 3, 2))

# Trace plots
plot(samples[, 1], type = "l", col = "blue",
     xlab = "Iteration", ylab = expression(theta[1]),
     main = paste("Theta 1 Trace\nAccept:", round(acceptance_rate, 3)))
abline(h = 0, col = "red", lty = 2)

plot(samples[, 2], type = "l", col = "blue",
     xlab = "Iteration", ylab = expression(theta[2]),
     main = "Theta 2 Trace")
abline(h = 0, col = "red", lty = 2)

# Scatter plot
plot(samples[, 1], samples[, 2],
     xlab = expression(theta[1]), ylab = expression(theta[2]),
     main = "Posterior Samples",
     col = rgb(0, 0, 1, alpha = 0.3), pch = 16)
points(theta_init[1], theta_init[2], col = "red", pch = 17, cex = 1.5)
legend("topright", "Start", col = "red", pch = 17, cex = 1.2)

# Histograms
hist(samples[, 1], breaks = 30, prob = TRUE,
     xlab = expression(theta[1]), main = "Marginal Theta 1",
     col = "lightblue")
curve(dnorm(x, 0, 1), add = TRUE, col = "red", lwd = 2)

hist(samples[, 2], breaks = 30, prob = TRUE,
     xlab = expression(theta[2]), main = "Marginal Theta 2",
     col = "lightblue")
curve(dnorm(x, 0, 1), add = TRUE, col = "red", lwd = 2)

# Autocorrelation
acf(samples[, 1], main = "Autocorrelation Theta 1", 
    ylab = "ACF", col = "blue", lwd = 2)

HMC Sampling from Bivariate Normal

Example 2: Bayesian Logistic Regression

# Simulate data
set.seed(456)
n <- 200
beta_true <- c(1, -0.5, 1.5)
X <- cbind(1, matrix(rnorm(n * 2), ncol = 2))
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(n, 1, p)

# Log posterior
log_posterior_logistic <- function(beta, X, y, prior_var = 10) {
    eta <- X %*% beta
    log_lik <- sum(y * eta - log(1 + exp(eta)))
    log_prior <- -0.5 * sum(beta^2) / prior_var
    return(log_lik + log_prior)
}

# Gradient
grad_log_posterior_logistic <- function(beta, X, y, prior_var = 10) {
    eta <- X %*% beta
    p <- 1 / (1 + exp(-eta))
    grad_lik <- t(X) %*% (y - p)
    grad_prior <- -beta / prior_var
    return(grad_lik + grad_prior)
}

# Wrapper functions
log_post_wrapper <- function(theta) {
    log_posterior_logistic(theta, X, y)
}

grad_log_post_wrapper <- function(theta) {
    grad_log_posterior_logistic(theta, X, y)
}

# Run HMC
theta_init <- c(0, 0, 0)
hmc_logistic <- hmc(
    log_posterior = log_post_wrapper,
    grad_log_posterior = grad_log_post_wrapper,
    theta_init = theta_init,
    n_samples = 5000,
    L_steps = 30,
    epsilon = 0.08,
    warmup = 1000,
    verbose = TRUE
)

## Starting HMC with 3 parameters
## Step size: 0.08 | Leapfrog steps: 30 
## Total iterations: 6000 
## Iteration 1000/6000 | Acceptance: 0.971
## Iteration 2000/6000 | Acceptance: 0.975
## Iteration 3000/6000 | Acceptance: 0.974
## Iteration 4000/6000 | Acceptance: 0.973
## Iteration 5000/6000 | Acceptance: 0.972
## Iteration 6000/6000 | Acceptance: 0.973
## 
## HMC Complete
## Acceptance rate: 0.974 
## Effective samples: 5000

cat(sprintf("\nLogistic Regression Acceptance rate: %.3f\n", 
            hmc_logistic$acceptance_rate))

## 
## Logistic Regression Acceptance rate: 0.974

# Results
par(mfrow = c(2, 3), mar = c(4, 4, 3, 2))
for (j in 1:3) {
    plot(hmc_logistic$samples[, j], type = "l",
         xlab = "Iteration", ylab = bquote(beta[.(j-1)]),
         main = paste("Beta", j-1), col = "blue")
    abline(h = beta_true[j], col = "red", lty = 2)
    
    hist(hmc_logistic$samples[, j], breaks = 30, prob = TRUE,
         xlab = bquote(beta[.(j-1)]), 
         main = paste("Posterior Beta", j-1),
         col = "lightblue")
    abline(v = beta_true[j], col = "red", lwd = 2)
}

HMC for Logistic Regression

cat("\nTrue beta:", beta_true, "\n")

## 
## True beta: 1 -0.5 1.5

cat("Posterior means:", colMeans(hmc_logistic$samples), "\n")

## Posterior means: 0.9498134 -0.602402 1.56464

cat("Posterior SDs:", apply(hmc_logistic$samples, 2, sd), "\n")

## Posterior SDs: 0.1926422 0.2001623 0.2619369

Example 3: Tuned HMC with Adaptive Step Size

#' Tuned HMC with Adaptive Step Size
hmc_tuned <- function(log_posterior, grad_log_posterior,
                      theta_init, n_samples, L_steps,
                      epsilon_init = 0.1, target_accept = 0.65,
                      adapt_steps = 1000, verbose = TRUE) {
    
    d <- length(theta_init)
    theta <- theta_init
    epsilon <- epsilon_init
    
    total_steps <- adapt_steps + n_samples
    samples <- matrix(NA, nrow = n_samples, ncol = d)
    acceptance <- logical(total_steps)
    epsilon_history <- numeric(total_steps)
    
    accept_count <- 0
    accept_rate_window <- numeric(100)
    window_idx <- 1
    
    if (verbose) {
        cat("Starting tuned HMC\n")
        cat("Initial epsilon:", epsilon, "\n")
        cat("Target acceptance:", target_accept, "\n")
    }
    
    for (t in 1:total_steps) {
        # Sample momentum
        r_star <- rnorm(d)
        theta_current <- theta
        
        # Leapfrog steps
        for (l in 1:L_steps) {
            result <- leapfrog(grad_log_posterior, theta, r_star, epsilon)
            theta <- result$theta
            r_star <- result$r
        }
        
        # Metropolis
        current_log_prob <- log_posterior(theta_current) - sum(r_star^2) / 2
        proposed_log_prob <- log_posterior(theta) - sum(r_star^2) / 2
        
        log_accept <- proposed_log_prob - current_log_prob
        
        if (is.finite(log_accept) && log(runif(1)) < min(0, log_accept)) {
            if (t > adapt_steps) {
                samples[t - adapt_steps, ] <- theta
            }
            acceptance[t] <- TRUE
            accept_count <- accept_count + 1
        } else {
            if (t > adapt_steps) {
                samples[t - adapt_steps, ] <- theta_current
            }
            theta <- theta_current
            acceptance[t] <- FALSE
        }
        
        # Adaptation
        if (t <= adapt_steps) {
            accept_rate_window[window_idx] <- as.numeric(acceptance[t])
            window_idx <- (window_idx %% 100) + 1
            
            recent_rate <- mean(accept_rate_window[1:min(t, 100)])
            
            if (t > 10 && is.finite(recent_rate)) {
                epsilon <- epsilon * exp(0.02 * (recent_rate - target_accept))
                epsilon <- max(min(epsilon, 1.0), 0.001)
            }
            
            epsilon_history[t] <- epsilon
        }
        
        if (verbose && t %% 500 == 0) {
            cat(sprintf("Iteration %d/%d | Epsilon: %.4f | Accept: %.3f\n",
                       t, total_steps, epsilon,
                       mean(acceptance[1:t])))
        }
    }
    
    final_rate <- if (n_samples > 0) sum(acceptance[(adapt_steps+1):total_steps]) / n_samples else 0
    
    return(list(
        samples = samples,
        acceptance_rate = final_rate,
        epsilon_final = epsilon,
        epsilon_history = epsilon_history,
        acceptance = acceptance
    ))
}

# Test tuned HMC
set.seed(123)
tuned_result <- hmc_tuned(
    log_posterior = log_posterior_mvn,
    grad_log_posterior = grad_log_posterior_mvn,
    theta_init = c(3, 3),
    n_samples = 4000,
    L_steps = 20,
    epsilon_init = 0.2
)

## Starting tuned HMC
## Initial epsilon: 0.2 
## Target acceptance: 0.65 
## Iteration 500/5000 | Epsilon: 0.5481 | Accept: 0.742
## Iteration 1000/5000 | Epsilon: 0.6579 | Accept: 0.713
## Iteration 1500/5000 | Epsilon: 0.6579 | Accept: 0.723
## Iteration 2000/5000 | Epsilon: 0.6579 | Accept: 0.729
## Iteration 2500/5000 | Epsilon: 0.6579 | Accept: 0.735
## Iteration 3000/5000 | Epsilon: 0.6579 | Accept: 0.740
## Iteration 3500/5000 | Epsilon: 0.6579 | Accept: 0.739
## Iteration 4000/5000 | Epsilon: 0.6579 | Accept: 0.739
## Iteration 4500/5000 | Epsilon: 0.6579 | Accept: 0.738
## Iteration 5000/5000 | Epsilon: 0.6579 | Accept: 0.737

cat("\n=== Tuned HMC Results ===\n")

## 
## === Tuned HMC Results ===

cat("Final acceptance rate:", round(tuned_result$acceptance_rate, 3), "\n")

## Final acceptance rate: 0.743

cat("Final epsilon:", round(tuned_result$epsilon_final, 4), "\n")

## Final epsilon: 0.6579

# Plot epsilon adaptation
par(mfrow = c(1, 2), mar = c(4, 4, 3, 2))
plot(tuned_result$epsilon_history, type = "l", col = "blue",
     xlab = "Iteration", ylab = "Epsilon",
     main = "Step Size Adaptation")
abline(h = tuned_result$epsilon_final, col = "red", lty = 2)

# Trace plot
plot(tuned_result$samples[, 1], type = "l", col = "blue",
     xlab = "Iteration", ylab = expression(theta[1]),
     main = "Trace After Adaptation")
abline(h = 0, col = "red", lty = 2)

Adaptive HMC Step Size

Diagnostics and Evaluation

Effective Sample Size (ESS)

#' Compute Effective Sample Size
effective_sample_size <- function(samples, max_lag = NULL) {
    if (is.matrix(samples)) {
        return(apply(samples, 2, effective_sample_size, max_lag))
    }
    
    n <- length(samples)
    if (is.null(max_lag)) max_lag <- min(100, n - 1)
    
    samples <- samples[!is.na(samples)]
    n <- length(samples)
    
    if (n < 2) return(NA)
    
    mean_sample <- mean(samples)
    
    acf_vals <- numeric(max_lag + 1)
    acf_vals[1] <- 1
    
    for (lag in 1:max_lag) {
        numerator <- sum((samples[1:(n-lag)] - mean_sample) * 
                        (samples[(lag+1):n] - mean_sample))
        denominator <- sum((samples - mean_sample)^2)
        acf_vals[lag + 1] <- if (denominator > 0) numerator / denominator else 0
    }
    
    acf_vals <- pmax(acf_vals, 0)
    ess <- n / (1 + 2 * sum(acf_vals[-1]))
    return(ess)
}

#' Gelman-Rubin R-hat Diagnostic
rhat <- function(chains) {
    n_chains <- length(chains)
    if (n_chains < 2) stop("Need at least 2 chains")
    
    n_samples <- nrow(chains[[1]])
    d <- ncol(chains[[1]])
    
    warmup <- floor(n_samples / 2)
    chains <- lapply(chains, function(x) x[(warmup+1):n_samples, , drop = FALSE])
    n_samples <- nrow(chains[[1]])
    
    rhat_vals <- numeric(d)
    
    for (j in 1:d) {
        within_var <- mean(sapply(chains, function(chain) var(chain[, j])))
        chain_means <- sapply(chains, function(chain) mean(chain[, j]))
        between_var <- n_samples * var(chain_means)
        
        var_hat <- (n_samples - 1) / n_samples * within_var + between_var / n_samples
        rhat_vals[j] <- sqrt(var_hat / within_var)
    }
    
    return(rhat_vals)
}

#' HMC Summary
hmc_summary <- function(hmc_result, true_values = NULL) {
    samples <- hmc_result$samples
    d <- ncol(samples)
    
    means <- colMeans(samples, na.rm = TRUE)
    sds <- apply(samples, 2, sd, na.rm = TRUE)
    quantiles <- apply(samples, 2, quantile, 
                      probs = c(0.025, 0.25, 0.5, 0.75, 0.975), 
                      na.rm = TRUE)
    
    ess <- effective_sample_size(samples)
    
    summary_table <- data.frame(
        Mean = means,
        SD = sds,
        ESS = ess,
        Q2.5 = quantiles[1, ],
        Q25 = quantiles[2, ],
        Q50 = quantiles[3, ],
        Q75 = quantiles[4, ],
        Q97.5 = quantiles[5, ]
    )
    
    if (!is.null(true_values)) {
        summary_table$True = true_values
        summary_table$Bias = summary_table$Mean - true_values
        summary_table$RMSE = sqrt(summary_table$Bias^2 + summary_table$SD^2)
    }
    
    return(summary_table)
}

# Compute diagnostics
summary <- hmc_summary(hmc_result, true_values = c(0, 0))
cat("\n=== HMC Summary Statistics ===\n")

## 
## === HMC Summary Statistics ===

print(round(summary, 4))

##      Mean     SD      ESS    Q2.5     Q25     Q50    Q75  Q97.5 True    Bias
## 1  0.0308 0.9963 270.1772 -1.9049 -0.6256  0.0153 0.7065 1.9771    0  0.0308
## 2 -0.0321 0.9913 272.8997 -2.0011 -0.6881 -0.0351 0.6311 1.9060    0 -0.0321
##     RMSE
## 1 0.9967
## 2 0.9918

# Multiple chains for R-hat
set.seed(123)
chain1 <- hmc(log_posterior_mvn, grad_log_posterior_mvn,
              theta_init = c(2, 2), n_samples = 1000, 
              L_steps = 20, epsilon = 0.15, verbose = FALSE)

set.seed(456)
chain2 <- hmc(log_posterior_mvn, grad_log_posterior_mvn,
              theta_init = c(-2, -2), n_samples = 1000, 
              L_steps = 20, epsilon = 0.15, verbose = FALSE)

set.seed(789)
chain3 <- hmc(log_posterior_mvn, grad_log_posterior_mvn,
              theta_init = c(0, 3), n_samples = 1000, 
              L_steps = 20, epsilon = 0.15, verbose = FALSE)

chains <- list(chain1$samples, chain2$samples, chain3$samples)
rhat_values <- rhat(chains)
cat("\nR-hat values:", round(rhat_values, 4), "\n")

## 
## R-hat values: 1.0003 0.9992

cat("R-hat < 1.1 indicates convergence:", all(rhat_values < 1.1), "\n")

## R-hat < 1.1 indicates convergence: TRUE

Comparison with Random Walk Metropolis

compare_samplers <- function(target_log_posterior, grad_target,
                             theta_init, n_samples, d) {
    
    # HMC
    cat("Running HMC...\n")
    hmc_result <- hmc(target_log_posterior, grad_target,
                       theta_init, n_samples,
                       L_steps = 20, epsilon = 0.15,
                       warmup = 500, verbose = FALSE)
    
    # Random Walk Metropolis
    cat("Running Random Walk Metropolis...\n")
    rwm_samples <- matrix(NA, nrow = n_samples, ncol = d)
    theta <- theta_init
    proposal_sd <- 0.5
    
    accepted <- 0
    for (i in 1:n_samples) {
        theta_prop <- theta + rnorm(d) * proposal_sd
        log_accept <- target_log_posterior(theta_prop) - target_log_posterior(theta)
        
        if (log(runif(1)) < min(0, log_accept)) {
            theta <- theta_prop
            accepted <- accepted + 1
        }
        rwm_samples[i, ] <- theta
    }
    
    # Compare
    par(mfrow = c(2, 2), mar = c(4, 4, 3, 2))
    
    # Trace plots
    plot(hmc_result$samples[, 1], type = "l", col = "blue",
         xlab = "Iteration", ylab = "Theta",
         main = paste("HMC (Accept:", round(hmc_result$acceptance_rate, 3), ")"))
    
    plot(rwm_samples[, 1], type = "l", col = "red",
         xlab = "Iteration", ylab = "Theta",
         main = paste("RWM (Accept:", round(accepted/n_samples, 3), ")"))
    
    # Autocorrelation
    acf(hmc_result$samples[, 1], main = "HMC ACF", 
        col = "blue", lwd = 2)
    acf(rwm_samples[, 1], main = "RWM ACF", 
        col = "red", lwd = 2)
    
    # Efficiency
    hmc_ess <- effective_sample_size(hmc_result$samples[, 1])
    rwm_ess <- effective_sample_size(rwm_samples[, 1])
    
    cat("\n=== Comparison ===\n")
    cat("HMC ESS:", round(hmc_ess, 0), "\n")
    cat("RWM ESS:", round(rwm_ess, 0), "\n")
    cat("HMC Efficiency:", round(hmc_ess / n_samples * 100, 1), "%\n")
    cat("RWM Efficiency:", round(rwm_ess / n_samples * 100, 1), "%\n")
    
    return(list(hmc = hmc_result, rwm = rwm_samples))
}

# Run comparison
comparison <- compare_samplers(
    target_log_posterior = log_posterior_mvn,
    grad_target = grad_log_posterior_mvn,
    theta_init = c(2, 2),
    n_samples = 5000,
    d = 2
)

## Running HMC...
## Running Random Walk Metropolis...

HMC vs Random Walk Metropolis

## 
## === Comparison ===
## HMC ESS: 1164 
## RWM ESS: 145 
## HMC Efficiency: 23.3 %
## RWM Efficiency: 2.9 %

Summary and Best Practices

Key Takeaways

HMC improves upon standard MCMC by using gradient information to guide sampling
The leapfrog integrator preserves volume and time-reversibility
The Metropolis step corrects for numerical errors in the integrator
HMC scales well to high dimensions compared to random walk methods
The Boltzmann distribution connection provides theoretical foundation
Energy = -log(posterior) explains why HMC explores high-probability regions

Tuning Guidelines

Parameter	Effect	Recommended Range
Step size (ε)	Controls integration accuracy	0.01 - 0.3
Number of steps (L)	Controls trajectory length	10 - 50
Target acceptance	Balances exploration vs computation	0.65 - 0.80

When to Use HMC

Good for: - High-dimensional posteriors (d > 10) - Problems where gradients are available - Complex correlated posteriors - Automatic inference engines

Avoid for: - Discrete parameters - Very rough or discontinuous posteriors - Problems where gradients are too expensive

Advanced Topics

NUTS (No-U-Turn Sampler): Automatically tunes L
Adaptive HMC: Tunes ε during warmup
Riemannian HMC: Uses position-dependent mass matrix
Multivariate HMC: Uses full covariance for momentum

Hamiltonian Monte Carlo (HMC) - Complete Guide with R Implementation

2026-07-09

Introduction to Hamiltonian Monte Carlo

What is HMC?

Key Problems HMC Solves:

Core Intuition: Physics Analogy

Mathematical Foundation

The Joint Distribution

Components:

The Hamiltonian:

Derivation of the HMC Joint Distribution

1. Understanding the Components

Component 1: The Target Posterior \(p(\theta | y)\)

Component 2: The Log-Posterior \(L(\theta)\)

Component 3: The Momentum Distribution \(\mathcal{N}(r | 0, I)\)

2. The Full Derivation

Step 1: Start with the augmented joint distribution

Step 2: Substitute the momentum distribution

Step 3: Express in terms of log-posterior

Step 4: Combine exponentials

Why Energy = Negative Log-Posterior

The Direct Mathematical Connection

Step 1: Start with the Joint Distribution

Step 2: Express in Exponential Form

Step 3: Combine Exponentials

Step 4: Define the Hamiltonian

The Physical Intuition

High Posterior = Low Potential Energy (Valleys)

Low Posterior = High Potential Energy (Mountains)

The Complete Picture

Mathematical Relationship

The Components

Why Is This Important?

1. Conservation of Energy

2. The Acceptance Probability

3. The Gradient Connection

The Boltzmann Distribution

What is the Boltzmann Distribution?

The Mathematical Form

Visualizing the Boltzmann Distribution

Simple Examples of the Boltzmann Distribution

Example 1: Two-State System

Example 2: Three-State System

The Connection to HMC

How HMC Uses the Boltzmann Distribution

Physical Interpretation

What Does Temperature Mean in HMC?

The Role of the Momentum

Why Is This So Powerful?

1. It Provides a Theoretical Foundation

2. It Explains Why HMC Works

3. It Gives Us Conservation Laws

4. It Connects to Temperature-Based Methods

Complete Example: Sampling from the Boltzmann Distribution

Summary of the Boltzmann Distribution in HMC

The Key Equation to Remember:

Key Insights:

The Leapfrog Integrator

Why Leapfrog?

Leapfrog Algorithm

The HMC Algorithm

Complete HMC Algorithm

R Implementation

Examples

Example 1: Sampling from Multivariate Normal

Example 2: Bayesian Logistic Regression

Example 3: Tuned HMC with Adaptive Step Size

Diagnostics and Evaluation

Effective Sample Size (ESS)

Comparison with Random Walk Metropolis

Summary and Best Practices

Key Takeaways

Tuning Guidelines

When to Use HMC

Advanced Topics