Lifted Metropolis X2X4.py

From Werner KRAUTH

Revision as of 13:47, 11 June 2024 Werner (Talk | contribs) ← Previous diff		Current revision Werner (Talk | contribs) (→References)
Line 39:		Line 39:
	==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	+	* 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]

---

Current revision

Contents 1 Context 2 Python program 3 Further information 4 References

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

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)
    if random.random() < math.exp(-delta_u): x = new_x
    else: sigma *= -1
    data.append(x)

plt.title('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