Python&数学基礎
メトロポリス法（MCMC）
プログラム（Python 3.10.4）
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import pandas as pd
import random
def metro_f(x):
return 1/math.sqrt(2*math.pi)*math.exp(-x**2/2)
def metro():
fig_sample = plt.figure(figsize=(4, 3))
times = 10000
skip_times = 10
t = times * skip_times
cnt = 0
theta = 0
s_data = []
for i in range(t):
theta_new = theta + random.uniform(-1,1)
alpha = min(1, metro_f(theta_new)/metro_f(theta))
r = random.uniform(0,1)
if r > alpha:
theta_new = theta
theta = theta_new
cnt += 1
if cnt % skip_times == 0:
s_data.append(theta)
x_data = [i/100 for i in range(-500,500)]
y_data = [metro_f(x_data[i]) for i in range(len(x_data))]
plt.plot(x_data,y_data,color="green")
plt.hist(s_data, bins=100, density=True, color="red")
plt.xlabel("x")
plt.ylabel("metro_f(x)")
plt.legend
plt.grid(True)
plt.show()