pyhmc: Hamiltonain Monte Carlo in Python¶

Introduction¶

Hamiltonian Monte Carlo or Hybrid Monte Carlo (HMC) is a Markov chain Monte Carlo (MCMC) algorithm. Hamiltonian dynamics can be used to produce distant proposals for the Metropolis algorithm, thereby avoiding the slow exploration of the state space that results from the diffusive behaviour of simple random-walk proposals. It does this by taking a series of steps informed by first-order gradient information.

This feature allows it to converge much more quickly to high-dimensional target distributions compared to simpler methods such as Metropolis, Gibbs sampling (and derivatives).

References¶

Neal, Radford. “MCMC using Hamiltonian dynamics.” Handbook of Markov Chain Monte Carlo 2 (2011).

Example¶

If you wanted to draw samples from a 5 dimensional Gaussian, you would do something like:

# define your probability distribution
import numpy as np
def logprob(x, ivar):
    logp = -0.5 * np.sum(ivar * x**2)
    grad = -ivar * x
    return logp, grad

# run the sampler
from pyhmc import hmc
ivar = 1. / np.random.rand(5)
samples = hmc(logprob, x0=np.random.randn(5), args=(ivar,), n_samples=1e4)

# Optionally, plot the results (requires an external package)
import triangle  # pip install triangle_plot
figure = triangle.corner(samples)
figure.savefig('triangle.png')