随机采样之接受拒绝采样

之前提到的逆变换采样（Inverse Transform Sampling）是一种生成随机样本的方法，它利用累积分布函数（CDF）的逆函数来生成具有特定分布的随机变量。以下是逆变换采样的缺点：

1. 计算复杂性：对于某些分布，找到累积分布函数（CDF）的逆函数可能是困难的，甚至是不可能的。
2. 效率问题：对于具有重尾分布的随机变量，逆变换采样可能非常低效，因为CDF的逆可能需要大量的计算。
3. 数值稳定性：在数值计算中，由于浮点数的精度限制，逆变换采样可能会引入误差，尤其是在CDF的值接近1时。

一、接受拒绝采样

接受-拒绝采样（Accept-Reject Sampling）方法是一种更为通用的采样方法，它可以用来生成具有任意分布的随机样本。这种方法不要求我们知道CDF的逆，而是利用一个简单的概率分布（称为提议分布）来生成样本，然后以一定的概率接受或拒绝这些样本。

接收-拒绝采样的基本步骤：

1. 选择提议分布$g(x)$：选择一个容易从中抽样的分布$g(x)$，并且确保对于所有的$x$，有$f(x)\le M\cdot g(x)$，其中$f(x)$是目标分布，$M$是一个正常数。

2. 抽样：从提议分布$g(x)$中抽取样本$x$和从均匀分布$U(0,1)$中抽取样本$u$。

3. 接受-拒绝条件：如果$u\le \frac{f(x)}{M\cdot g(x)}$，则接受$x$作为目标分布$f(x)$的一个样本；否则拒绝$x$。

接受拒绝采样可以使用下图进行表示（图片来源：【数之道】马尔可夫链蒙特卡洛方法是什么？十五分钟理解这个数据科学难点）。

二、接受拒绝采样证明

要证明接收-拒绝采样确实产生服从目标分布$f(x)$的样本，我们需要证明对于所有的$x$，有：
$$P(X=x)=f(x)\text{\hspace{1em}(1)}$$

其中$P(X=x)$是样本$x$被接受的概率。

证明：

1. 接受概率：样本$x$被接受的概率是$\frac{f(x)}{M\cdot g(x)}$，因为$u$是从$U(0,1)$中抽取的。

2. 联合概率：样本$x$从提议分布$g(x)$中抽取的概率是$g(x)$，并且$u$在$[0,\frac{f(x)}{M\cdot g(x)})$区间的概率是$\frac{f(x)}{M\cdot g(x)}$。因此，联合概率是：

   $$P(X=x,U\le \frac{f(x)}{M\cdot g(x)})=g(x)\cdot \frac{f(x)}{M\cdot g(x)}=\frac{f(x)}{M}\text{\hspace{1em}(2)}$$

3. 边缘概率：现在我们需要计算$X$的边缘概率$P(X=x)$，即样本$x$被接受的总概率。由于$u$是均匀分布的，我们可以将联合概率在$u$的所有可能值上积分：

   $$P(X=x)={\int }_0^{1}P(X=x,U=u)du={\int }_0^1\frac{f(x)}{M}du=\frac{f(x)}{M}\cdot {\int }_0^{1}du=\frac{f(x)}{M}\text{\hspace{1em}(3)}$$

4. 归一化常数：由于$M$是使得$f(x)\le M\cdot g(x)$对所有$x$成立的最小常数，我们可以将上式中的$M$移到$f(x)$的定义中，从而得到：

   $$P(X=x)=f(x)\text{\hspace{1em}(4)}$$
   这就证明了接收-拒绝采样确实产生了服从目标分布$f(x)$的样本。

三、接受拒绝采样模拟

借用作者anshuai_aw1的例子，设我们需要采样的pdf为：
$$f(x)=0.3\exp \left(-(x-0.3{)}^2\right)+0.7\exp \left(-(x-2{)}^2/0.3\right)\text{\hspace{1em}(5)}$$
其归一化常数为$Z=1.2113$, 参考分布为 $g(x)=N(\mu =1.4,{\sigma }^2=(1.2^2))$, $M=2.5$, 以确保 $M\cdot g(x)\ge f(x)$。采样的代码如下：

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return (0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3))/1.2113
x = np.arange(-4.,6.,0.01)
plt.plot(x,f(x),color = "red")

size = int(1e+07)
mu = 1.4
sigma = 1.2
M = 2.5

x = np.random.normal(loc = mu,scale = sigma, size = size)
g_x = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*(x-mu)**2/sigma**2)
u = np.random.uniform(low = 0, high = M*g_x, size = size)  #在[0,M*g_x]中均匀采样
fx =  0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3)
sample = x[u <= fx] # u < fx(x)
plt.hist(sample,bins=150, density=True, edgecolor='black')
plt.show()

结果如下，其中红色曲线的是公式(5)所示pdf的图像，蓝色区域是采样结果，可见采样结果跟真实分布几乎一致。

参考资料：

[1]【数之道】马尔可夫链蒙特卡洛方法是什么？十五分钟理解这个数据科学难点
[2] 逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解