Hamiltonian Monte Carlo

The Hamiltonian Monte Carlo algorithm (originally known as hybrid Monte Carlo) is a Markov chain Monte Carlo method for obtaining a sequence of random samples whose distribution converges to a target probability distribution that is difficult to sample directly. This sequence can be used to estimate integrals of the target distribution, such as expected values and moments.

Hamiltonian Monte Carlo corresponds to an instance of the Metropolis–Hastings algorithm, with a Hamiltonian dynamics evolution simulated using a time-reversible and volume-preserving numerical integrator (typically the leapfrog integrator) to propose a move to a new point in the state space. Compared to using a Gaussian random walk proposal distribution in the Metropolis–Hastings algorithm, Hamiltonian Monte Carlo reduces the correlation between successive sampled states by proposing moves to distant states which maintain a high probability of acceptance due to the approximate energy conserving properties of the simulated Hamiltonian dynamic when using a symplectic integrator. The reduced correlation means fewer Markov chain samples are needed to approximate integrals with respect to the target probability distribution for a given Monte Carlo error.

The algorithm was originally proposed by Simon Duane, Anthony Kennedy, Brian Pendleton and Duncan Roweth in 1987 for calculations in lattice quantum chromodynamics. In 1996, Radford M. Neal showed how the method could be used for a broader class of statistical problems, in particular artificial neural networks. But the burden of having to provide the algorithm with gradients of the model graph delayed its wider adoption in statistics and other quantitative disciplines, until in the mid-2010s the developers of Stan implemented HMC in combination with automatic differentiation.

Algorithm

Suppose the target distribution to sample is

f

(

x

)

{\displaystyle f(\mathbf {x} )}

for

x

∈

R

d

{\displaystyle \mathbf {x} \in \mathbb {R} ^{d}}

(

d

≥

1

{\displaystyle d\geq 1}

) and a chain of samples

X

0

,

X

1

,

X

2

,

…

{\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }

is required.

Hamilton's equations are

d

x

i

d

t

=

∂

H

∂

p

i

and

d

p

i

d

t

=

−

∂

H

∂

x

i

{\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}\quad {\text{and}}\quad {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}

where

x

i

{\displaystyle x_{i}}

and

p

i

{\displaystyle p_{i}}

are the

i

{\displaystyle i}

th component of the position and momentum vector respectively and

H

{\displaystyle H}

is the Hamiltonian.