import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
def final_dist_prob(value):
return norm.pdf(value, loc=3, scale=2)
num_samples = 1000
t = 0
pi = 0 # 初始状态值
samples = []
while t < num_samples:
pi_star = norm.rvs(loc=pi, scale=1, size=1)[0] # 从条件概率分布中采样
alpha = min(1, final_dist_prob(pi_star) / (final_dist_prob(pi) + 10e-8)) # 计算接受率
u = np.random.rand()
if u < alpha:
samples.append(pi_star)
pi = pi_star
t += 1
x_range = np.linspace(-3, 9, 200)
plt.plot(x_range, final_dist_prob(x_range))
plt.hist(samples, bins=50, normed=1, alpha=0.5)
plt.show()