第19章 马尔可夫链蒙特卡洛法 - M-H采样python实现 - 实现norm_dist_prob函数

import numpy as np
import random

def norm_dist_prob(theta):
    return np.exp(-theta**2 / 2) / np.sqrt(2 * np.pi)

def mh_sampling(n_samples=10000, sigma=1.0):
    samples = []
    theta = 0.0
    for _ in range(n_samples):
        theta_new = theta + np.random.normal(0, sigma)
        alpha = min(1, norm_dist_prob(theta_new) / norm_dist_prob(theta))
        if random.random() < alpha:
            theta = theta_new
        samples.append(theta)
    return samples

samples = mh_sampling(1000)
print("采样完成，均值:", np.mean(samples))