Road three

Hamiltonian Monte Carlo, one trajectory at a time

HMC treats the negative log posterior as a landscape and the coefficients as a particle rolling across it. Each iteration gives the particle a random push, simulates its motion with the leapfrog integrator, then accepts or rejects where it lands. This is the original HMC.fn, on the corrected two-parameter posterior.

Click or tap the plot to drop the particle somewhere else. While this panel has focus: n next proposal, p auto-run.

This proposal

Press Next proposal. A random momentum is drawn, the particle slides along the log-posterior surface for L leapfrog steps, and a Metropolis test decides whether to keep the end point.

Chain so far

Iteration: 1
Current position: (0.393, -0.160)
Acceptance rate: –

Seed 90125, L = 20, ε = 1/20: every point here is the same draw R's original HMC.fn gives on the corrected data (scripts/r-reference.R), burn-in included.

Integrator

Leapfrog steps L20

Keep ε·L = 1, as the original code does

Step size ε (set by L)0.050

Push ε past about 0.14 and the leapfrog integrator goes unstable along the narrow direction of the posterior: energy is no longer conserved and proposals get rejected.

What one iteration does

Draw a momentum $\phi \sim N(0, M)$ with the identity mass matrix. Then repeat L times (L = 20, ε = 1/L in the original):

\begin{aligned} \phi &\leftarrow \phi + \tfrac{\varepsilon}{2}\, \nabla_\beta \log p(\beta \mid y) \\ \beta &\leftarrow \beta + \varepsilon\, M^{-1} \phi \\ \phi &\leftarrow \phi + \tfrac{\varepsilon}{2}\, \nabla_\beta \log p(\beta \mid y) \end{aligned}

For logistic regression with a flat prior the gradient is simply $X^\top (y - n \odot p)$ , the line crossprod(X, y - n * p.b) in the R code.

Why it is accepted so often

The leapfrog scheme almost conserves the total energy $H(\beta, \phi) = -\log p(\beta \mid y) + \tfrac12 \phi^\top M^{-1} \phi$ , so the Metropolis ratio $\exp(H_0 - H_L)$ is usually close to one. In the corrected run 95.9% of proposals were accepted, against 35.8% for random-walk Metropolis, and the effective sample size per draw is more than ten times higher.

The energy plot in the “This proposal” panel shows this directly: the dashed line is the starting energy, and the further the curve ends from it, the less likely the jump is to be kept. Compare all four methods in the explorer.