Factor Metropolis X2X4 patch.py

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
Revision as of 23:01, 10 June 2024
Werner (Talk | contribs)

← Previous diff
Current revision
Werner (Talk | contribs)

Line 1: Line 1:
 +==Context==
 +This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).
 +
 +==Python program==
 +
 +
import math import math
import random import random
Line 36: Line 42:
plt.legend(loc='upper right') plt.legend(loc='upper right')
plt.show() plt.show()
 +
 +==Further information==
 +==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

Current revision

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 u2(pos): return 0.5 * pos ** 2

def u4(pos): return 0.25 * pos ** 4

x = 0.0
delta = 0.1

data = []
n_samples = 10 ** 6
for i in range(n_samples):
    new_x = x + random.uniform(-delta, delta)
    delta_u2 = u2(new_x) - u2(x)
    delta_u4 = u4(new_x) - u4(x)
    Upsilon2 = random.uniform(0.0, 1.0)
    Upsilon4 = random.uniform(0.0, 1.0)
    if Upsilon2 < math.exp(-delta_u2) and Upsilon4 < math.exp(-delta_u4):
        x = new_x
    data.append(x)
plt.title('Factorized Metropolis algorithm (patch), 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

Personal tools