Factor ZigZag X2X4.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 06:51, 11 June 2024 Werner (Talk | contribs) ← Previous diff |
Revision as of 07:03, 11 June 2024 Werner (Talk | contribs) Next diff → |
||
| 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 24: | Line 29: | ||
| sigma *= -1 | sigma *= -1 | ||
| + | ==Further information== | ||
| + | ==References== | ||
| plt.title('Factor-Zig-Zag algorithm, anharmonic oscillator') | plt.title('Factor-Zig-Zag algorithm, anharmonic oscillator') | ||
| plt.xlabel('$x$') | plt.xlabel('$x$') | ||
Revision as of 07:03, 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
x = 0.0
time_ev = 0.0
sigma = 1 if x <= 0.0 else -1
data = []
n_samples = 10 ** 6
while len(data) < n_samples:
delta_u2 = -math.log(random.uniform(0.0, 1.0))
delta_u4 = -math.log(random.uniform(0.0, 1.0))
new_x2 = math.sqrt(2.0 * delta_u2)
new_x4 = (4.0 * delta_u4) ** 0.25
new_x = sigma * min(abs(new_x2), abs(new_x4))
new_time_ev = time_ev + abs(new_x - x)
for t in range(math.ceil(time_ev), math.floor(new_time_ev) + 1):
data.append(x + sigma * (t - time_ev))
x = new_x
time_ev = new_time_ev
sigma *= -1
Further information
References
plt.title('Factor-Zig-Zag 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()
