Levy harmonic.py

From Werner KRAUTH

Contents 1 Context 2 Python program 3 Further information 4 References

Context

This page is part of my public lectures on Algorithms and Computations in theoretical physics at the University of Oxford (see the announcement for a syllabus). For more information, see the main page of the public lectures. It presents the hamiltonian Monte Carlo algorithm invented by Duane et al. (1988). For an in-depth discussion, see my below reference. Other programs in the same series are the Metropolis algorithm and the Lévy construction.

Python program

import math, random

def U(x):
    U = 0.0
    for k in range(N):
        k_minus = (k - 1) % N
        x_minus = x[k_minus]
        if k == 0: x_minus -= L
        U += (x[k] - x_minus) ** 2 / 2.0
    return U

N = 8
L = 16
delta = 1.0

Iter = 1000000

for iter in range(Iter):
    xstart = random.uniform(0.0, L)
    xend = xstart + L
    x = [xstart]
    for k in range(1, N):        # loop over internal slices
        x_mean = ((N - k) * x[k - 1] + xend) / (1.0 + N - k)
        sigma = math.sqrt(1.0 / (1.0  + 1.0 / (N - k) ))
        x.append(random.gauss(x_mean, sigma))
    x.append(xend)

Further information

References

Krauth, W. Hamiltonian Monte Carlo vs. event-chain Monte Carlo: an appraisal of sampling strategies beyond the diffusive regime. ArXiv: 2411.11690