Bounded Lifted Metropolis X2X4.py

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
Revision as of 07:06, 11 June 2024
Werner (Talk | contribs)

← Previous diff
Revision as of 13:50, 11 June 2024
Werner (Talk | contribs)
(References)
Next diff →
Line 56: Line 56:
==Further information== ==Further information==
==References== ==References==
 +* Tartero, G., Krauth, W. Concepts in Monte Carlo sampling, Am. J. Phys. 92, 65–77 (2024) [https://arxiv.org/pdf/2309.03136 arXiv:2309.03136]

Revision as of 13:50, 11 June 2024

Contents

Context

This page is part of my 2024 Beg Rohu Lectures on "The second Markov chain revolution" at the Summer School "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

Python program

import math
import random
import matplotlib.pyplot as plt
def u(x):
    return x ** 2 / 2.0 + x ** 4 / 4.0

def u_bound(pos):
    floor = math.floor(abs(pos))
    ceiling = math.ceil(abs(pos))
    u_floor = 0
    for n in range(floor + 1):
        u_floor += n + n ** 3
    u_ceiling = u_floor + ceiling + ceiling ** 3
    u_pos = u_floor
    if floor != abs(pos):
        m = (u_ceiling - u_floor) / (ceiling - floor)
        u_pos = m * (abs(pos) - floor) + u_floor
    return u_pos

x = 0.0
delta = 0.1
sigma = random.choice([-1, 1])
data = []
n_samples = 10 ** 6
for i in range(n_samples):
    new_x = x + sigma * random.uniform(0.0, delta)
    delta_u = u(new_x) - u(x)
    delta_u_tilde = u_bound(new_x) - u_bound(x)
    if random.uniform(0.0, 1.0) < min(1.0, math.exp(-delta_u_tilde)):
        x = new_x
    else:
        if random.uniform(0.0, 1.0) > (1.0 - math.exp(-delta_u)) / (1.0 - math.exp(-delta_u_tilde)):
            x = new_x
        else:
            sigma *= -1
    data.append(x)

plt.title('Bounded-Lifted Metropolis algorithm, anharmonic oscillator')
plt.xlabel('$x$')
plt.ylabel('$\pi(x)$')
plt.hist(data, bins=100, density=True, label='data')
XValues = []
YValues = []
for i in range(-1000,1000):
    x = i / 400.0
    XValues.append(x)
    YValues.append(math.exp(- u(x)) / 1.93525)
plt.plot(XValues, YValues, label='theory')
plt.legend(loc='upper right')
plt.show()

Further information

References

  • Tartero, G., Krauth, W. Concepts in Monte Carlo sampling, Am. J. Phys. 92, 65–77 (2024) arXiv:2309.03136
Personal tools