连续性分布的MH抽样法python

Metropolis-Hastings(MH)抽样法是一种用于从任意连续分布中进行抽样的方法。在这个算法中，我们构建一个马尔可夫链，使得其平稳分布是我们想要抽样的分布。具体步骤如下：

1. 选择一个合适的提议分布$q(x|y)$，其中$x$是我们想要抽样的变量，$y$是当前状态。
2. 从提议分布$q(x|y)$中抽样得到一个候选状态$x'$。
3. 计算接受概率$\alpha=\min\left(1,\frac{p(x')q(y|x')}{p(y)q(x'|y)}\right)$，其中$p(x)$是我们想要抽样的分布，$p(y)$是当前状态下的概率密度函数。
4. 以概率$\alpha$接受候选状态$x'$，否则保持当前状态$y$不变。
5. 重复步骤2-4直到得到足够多的样本。

下面是一个使用MH抽样法抽样标准正态分布的Python代码：

import numpy as np
import scipy.stats as stats

def mh_sampler(p, q, x0, n):
    """
    p: 目标分布的概率密度函数
    q: 提议分布的概率密度函数
    x0: 初始状态
    n: 抽样次数
    """
    samples = [x0]
    x = x0
    for i in range(n):
        x_prime = q(x)
        alpha = min(1, p(x_prime)/p(x) * q(x, x_prime)/q(x_prime, x))
        u = np.random.uniform()
        if u < alpha:
            x = x_prime
        samples.append(x)
    return np.array(samples[1:])

# 目标分布：标准正态分布
p = stats.norm.pdf
# 提议分布：以当前状态为均值，标准差为1的正态分布
q = lambda x, x_prime: stats.norm.pdf(x_prime, loc=x, scale=1)
# 初始状态
x0 = 0
# 抽样次数
n = 10000

samples = mh_sampler(p, q, x0, n)

运行以上代码后，我们得到了10000个标准正态分布的样本。可以使用这些样本来估计期望、方差等统计量。注意，使用MH抽样法抽样连续分布时，提议分布$q(x|y)$的选择对抽样效率有很大影响。通常需要通过试验来确定一个合适的提议分布。