生成MCMC 的python代码

以下是一个简单的 Metropolis-Hastings MCMC 的 Python 代码示例，用于从给定分布中生成样本：

import numpy as np
import matplotlib.pyplot as plt

# 定义目标分布的概率密度函数，这里以标准正态分布为例
def target_distribution(x):
    return np.exp(-0.5 * x ** 2) / np.sqrt(2 * np.pi)

# 定义转移概率分布的概率密度函数，这里以标准正态分布为例
def proposal_distribution(x, sigma=1.0):
    return np.exp(-0.5 * ((x - sigma) / sigma) ** 2) / np.sqrt(2 * np.pi * sigma ** 2)

# 定义 MCMC 的 Metropolis-Hastings 步骤
def metropolis_hastings_step(x, target_distribution, proposal_distribution, sigma=1.0):
    x_new = np.random.normal(x, sigma)
    acceptance_ratio = (target_distribution(x_new) * proposal_distribution(x, sigma)) / (target_distribution(x) * proposal_distribution(x_new, sigma))
    u = np.random.uniform()
    if u < acceptance_ratio:
        return x_new
    else:
        return x

# 运行 MCMC 的 Metropolis-Hastings 算法
def metropolis_hastings(n_samples, initial_value, target_distribution, proposal_distribution, sigma=1.0):
    samples = [initial_value]
    for i in range(n_samples):
        sample = metropolis_hastings_step(samples[-1], target_distribution, proposal_distribution, sigma)
        samples.append(sample)
    return np.array(samples)

# 生成样本
samples = metropolis_hastings(n_samples=10000, initial_value=0.0, target_distribution=target_distribution, proposal_distribution=proposal_distribution)

# 绘制样本分布
plt.hist(samples, bins=50, density=True, alpha=0.5)
x = np.linspace(-5, 5, 100)
plt.plot(x, target_distribution(x), 'r-', lw=2, label='Target distribution')
plt.plot(x, proposal_distribution(x), 'g-', lw=2, label='Proposal distribution')
plt.legend()
plt.show()

这个代码示例中实现了一个简单的 Metropolis-Hastings MCMC 算法，用于从标准正态分布中生成样本。在实际应用中，需要根据实际情况选择适当的目标分布和转移概率分布，以及调整参数，才能获得好的采样效果。