David MacKay. Information Theory, Inference, and Learning Algorithms. Algorithm 30.1.

HMC.py

import numpy as np
import matplotlib
matplotlib.use('agg')
from pylab import *

def findE(x):
    return 250.25 * (x[0] * x[0] + x[1] * x[1]) \
        - 449.5 * x[0] * x[1]

def gradE(x):
    return np.array([500.5 * x[0] - 449.5 * x[1],
                     500.5 * x[1] - 449.5 * x[0]])

L = 2000
D = 2
Tau = 19
epsilon = 0.055

x = np.array([-0.5, 0.5])

g = gradE(x)
E = findE(x)

for l in range(L):
    p = np.random.normal(0.0, 1.0, D)
    H = np.dot(p, p) / 2 + E

    xnew = x
    gnew = g
    for tau in range(Tau):

        p = p - epsilon * gnew / 2
        xnew = xnew + epsilon * p
        gnew = gradE(xnew)
        p = p - epsilon * gnew / 2

    Enew = findE(xnew)
    Hnew = np.dot(p, p) / 2 + Enew
    dH = Hnew - H

    if dH < 0:
        accept = True
    elif np.random.uniform() < np.exp(- dH):
        accept = True
    else:
        accept = False

    if accept:
        g = gnew
        x = xnew
        E = Enew
        plt.plot(x[0], x[1], 'kx')

xlim(-0.6, 0.6)
ylim(-0.6, 0.6)
savefig('test.png')