MCMC and Hamiltonian Monte Carlo

The posterior \(p(\theta \mid x) = p(x \mid \theta)p(\theta) / p(x)\) has an intractable denominator: an integral over the whole parameter space. Markov chain Monte Carlo makes that irrelevant. Build a Markov chain whose stationary distribution is the posterior, run it, and treat the states as samples. Every acceptance decision involves a ratio of posterior densities, and the normalising constant cancels.

Metropolis-Hastings is the base case. Propose \(\theta'\) from \(q(\theta' \mid \theta)\), accept with probability

and otherwise stay put. Detailed balance guarantees the posterior is stationary, and the chain converges to it under mild conditions. Gibbs sampling is the special case where you cycle through the coordinates sampling each from its exact conditional, with acceptance always 1.

Why random walks fail in high dimensions

A random-walk proposal has to be tuned to a step size \(\epsilon\). Too large and almost everything is rejected; too small and the chain crawls. In \(d\) dimensions the optimal acceptance rate for a random-walk Metropolis on a Gaussian target is about 0.234, and the chain's efficiency scales roughly as \(1/d\): exploring the posterior takes \(O(d^2)\) steps because it is diffusive, covering distance like \(\sqrt{\text{steps}}\).

The deeper problem is the geometry of the typical set. In high dimensions, almost all posterior mass sits in a thin shell away from the mode, not at the mode. A proposal that steps in a random direction almost always steps out of that shell into a region of much lower density, and gets rejected. The number of directions that stay inside the shell shrinks as \(d\) grows.

What Hamiltonian dynamics changes

Hamiltonian Monte Carlo introduces an auxiliary momentum \(r\) and treats the negative log posterior as a potential energy:

Sample a fresh momentum from a Gaussian, then simulate Hamiltonian dynamics for some time using the leapfrog integrator, which is symplectic and time-reversible. Because Hamiltonian flow conserves \(H\) exactly and leapfrog conserves it nearly, the proposal at the end of a long trajectory has an acceptance probability close to 1 despite being far from the start. Gradient information means the trajectory follows the contours of the typical set instead of stepping blindly out of it (Neal, 2011, MCMC using Hamiltonian dynamics, Handbook of Markov Chain Monte Carlo).

The remaining tuning problem is trajectory length: too short is a random walk again, too long wastes gradient evaluations retracing. The No-U-Turn Sampler builds the trajectory adaptively, doubling it until it starts to double back, and adapts the step size during warmup (Hoffman and Gelman, 2014, The No-U-Turn Sampler, JMLR 15). This is what Stan, PyMC and NumPyro run by default, and it is why HMC became usable without an expert tuning it.

Diagnostics that actually detect failure

There is no test for convergence, only tests for detectable non-convergence. Three are standard.

\(\hat{R}\) compares within-chain and between-chain variance across several chains started from dispersed points; values above about 1.01 indicate the chains have not mixed. Effective sample size accounts for autocorrelation, and 1,000 draws with an ESS of 40 carry roughly 40 draws' worth of information, which is what should be reported.

Divergent transitions in HMC are the most informative signal and have no analogue in random-walk methods. A divergence means the leapfrog integrator's energy error exploded, which happens where the posterior has high curvature the step size cannot resolve. Divergences clustered in one region mean the sampler is systematically failing to explore it, and the resulting bias is not visible in the trace plot.

When it breaks

Multimodality defeats it. A chain in one mode separated by a low-density region will not cross it in any reasonable time, and every diagnostic looks healthy because the chain is mixing beautifully within its mode. Running several chains from dispersed inits is the only cheap detection; parallel tempering is the expensive fix.

Funnel geometry needs reparameterisation. Hierarchical models produce a funnel: the scale parameter controls the width of the group effects, so the posterior has a narrow neck no single step size can handle. Non-centred parameterisation, sampling standardised effects and scaling them afterwards, changes the geometry rather than the model and is the standard resolution.

Cost scales with the likelihood, not with the sampler. Each leapfrog step needs a gradient of the full log posterior, so a 10,000-step trajectory budget over a million-row dataset is a million-row pass per step. This is why variational inference and subsampling methods exist, and why Bayesian treatment of large neural networks remains largely approximate.

Discrete parameters have no gradient. HMC cannot sample them, so they must be marginalised analytically where possible, which is why Stan requires it, or handled by a separate mechanism inside a mixed sampler.