metropolis.py

import seaborn as sns

def action_fn(x):
    return 0.5*x*x  

num_of_samples = 100000
step_size = 0.5
x=0.
num_of_accepted = 0

sampled_x = [0]
for i in range(num_of_samples):
    dx = np.random.rand()
    dx = (dx - 0.5)*step_size*2.
    pre_ = action_fn(x)
    next_ = action_fn(x+dx)
    r = np.random.rand()
    if r < np.exp(pre_ - next_):
        x = x + dx
        num_of_accepted += 1
        
    sampled_x.append(x)
    
print(f"Prop. Accepted {num_of_accepted/num_of_samples}")
print(f"E[x] {sum(sampled_x)/num_of_samples}")
print(f"E[x^2] {sum([i*i for i in sampled_x])/num_of_samples}")
sns.displot(sampled_x)