Levy harmonic.py

From Werner KRAUTH

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 Lévy construction, in other words the direct sampling algorithm invented by P. Lévy in 1940. For an in-depth discussion, see my below reference. Other programs in the same series are the Metropolis algorithm and the Hamiltonian Monte Carlo algorithm, both for the harmonic chain.

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