Bounded Lifted Metropolis X2X4.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 06:39, 11 June 2024 Werner (Talk | contribs) ← Previous diff |
Current revision Werner (Talk | contribs) (→Python program) |
||
| 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 4: | Line 9: | ||
| def u(x): | def u(x): | ||
| return x ** 2 / 2.0 + x ** 4 / 4.0 | return x ** 2 / 2.0 + x ** 4 / 4.0 | ||
| - | + | ||
| def u_bound(pos): | def u_bound(pos): | ||
| floor = math.floor(abs(pos)) | floor = math.floor(abs(pos)) | ||
| Line 17: | Line 22: | ||
| u_pos = m * (abs(pos) - floor) + u_floor | u_pos = m * (abs(pos) - floor) + u_floor | ||
| return u_pos | return u_pos | ||
| - | + | ||
| x = 0.0 | x = 0.0 | ||
| delta = 0.1 | delta = 0.1 | ||
| sigma = random.choice([-1, 1]) | sigma = random.choice([-1, 1]) | ||
| - | |||
| data = [] | data = [] | ||
| n_samples = 10 ** 6 | n_samples = 10 ** 6 | ||
| Line 28: | Line 32: | ||
| delta_u = u(new_x) - u(x) | delta_u = u(new_x) - u(x) | ||
| delta_u_tilde = u_bound(new_x) - u_bound(x) | delta_u_tilde = u_bound(new_x) - u_bound(x) | ||
| - | fil = math.exp(-delta_u) | ||
| if random.uniform(0.0, 1.0) < min(1.0, math.exp(-delta_u_tilde)): | if random.uniform(0.0, 1.0) < min(1.0, math.exp(-delta_u_tilde)): | ||
| x = new_x | x = new_x | ||
| Line 37: | Line 40: | ||
| sigma *= -1 | sigma *= -1 | ||
| data.append(x) | data.append(x) | ||
| - | + | ||
| plt.title('Bounded-Lifted Metropolis algorithm, anharmonic oscillator') | plt.title('Bounded-Lifted Metropolis algorithm, anharmonic oscillator') | ||
| plt.xlabel('$x$') | plt.xlabel('$x$') | ||
| Line 51: | Line 54: | ||
| plt.legend(loc='upper right') | plt.legend(loc='upper right') | ||
| plt.show() | plt.show() | ||
| - | ~ | + | |
| + | This program takes a [[Bernoulli two pebbles patch.py| two-pebble decision]] that we already encountered as an "advanced" algorithm to sample the Bernoulli distribution. | ||
| + | |||
| + | ==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 |
[edit]
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).
[edit]
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()
This program takes a two-pebble decision that we already encountered as an "advanced" algorithm to sample the Bernoulli distribution.
[edit]
Further information
[edit]
References
- Tartero, G., Krauth, W. Concepts in Monte Carlo sampling, Am. J. Phys. 92, 65–77 (2024) arXiv:2309.03136
