メトロポリス法でガウス積分を求める

import random
import math
import matplotlib.pyplot as plt
from tqdm import tqdm

niters = [int(1e3),int(1e5),int(1e7)]
# niter = 1000
for niter in tqdm(niters):
    step_size = 0.5

    x = 0
    naccept = 0
    xs = []
    for i in range(niter):
        backup_x = x
        action_init = 0.5*x**2
        dx = random.random()
        dx = (dx-0.5)*step_size*2
        x = x+dx
        action_fin = 0.5*x**2
        metropolis = random.random()
        if (math.exp(action_init-action_fin) > metropolis):
            # accept
            naccept = naccept+1
        else:
            x = backup_x
        # print(f"{x},{naccept}/{i}")
        xs.append(x)
    plt.figure()
    plt.hist(xs,100)

plt.show()

f:id:seinzumtode:20220117171741p:plain