pyhmc: Hamiltonain Monte Carlo in Python¶
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).
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')