Bounded Lifted Metropolis X2X4.py

From Werner KRAUTH

Revision as of 06:39, 11 June 2024 Werner (Talk | contribs) ← Previous diff		Revision as of 06:40, 11 June 2024 Werner (Talk | contribs) Next diff →
Line 4:		Line 4:
	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 17:
	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 37:		Line 37:
	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 51:
	plt.legend(loc='upper right')		plt.legend(loc='upper right')
	plt.show()		plt.show()
-	~

---

Revision as of 06:40, 11 June 2024

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)
    fil = math.exp(-delta_u)
    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()