逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解

在蒙特卡罗方法中，有一个关键的问题需要解决，即如何基于概率密度函数去采的 n n n个 x x x的样本集。

逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)就是用于解决这个问题的。

其中，对于常见的分布，如均匀分布，高斯分布，指数分布，t分布，F分布，Beta分布，Gamma分布等，可以采用逆采样的方法进行采样；不过很多时候，我们的 x x x的概率分布不是常见的分布，这些分布的概率分布函数CDF 不可逆，因此没有办法用逆采样来采样，这意味着我们没法方便的得到这些非常见的概率分布的样本集。拒绝采样就是用来解决这个问题。

一、逆采样(Inverse Sampling)

我们知道，对于常见的均匀分布 u n i f o r m ( 0 , 1 ) uniform(0,1) uniform(0,1)是非常容易采样样本的，一般通过线性同余发生器可以很方便的生成(0,1)之间的伪随机数样本。而其他常见的概率分布，无论是离散的分布还是连续的分布，它们的样本都可以通过 u n i f o r m ( 0 , 1 ) uniform(0,1) uniform(0,1)的样本转换而得。那么应该如何得到呢？这就是逆采样。

下面我们以指数分布为例，说明如何通过均匀分布来采样服从指数分布的样本集。指数分布的概率密度函数PDF为：
p e x p ( x ) = { λ e x p ( − λ x ) , x ≥ 0 0 , x &lt; 0 (1) p_{exp}(x) = \begin{cases} \lambda exp(-\lambda x) &amp;,x\ge0\\ 0&amp;, x\lt0 \end{cases}\qquad \text{(1)} pexp​(x)={λexp(−λx)0​,x≥0,x<0​(1)

那么它的概率分布函数CDF为：
F ( x ) = ∫ − ∞ x p ( x ) d x (2) F(x) = \int_{-\infty}^xp(x)dx\qquad \text{(2)} F(x)=∫−∞x​p(x)dx(2)

下图为指数分布和均匀分布的CDF图。从左图上看，在 x ≥ 0 x\ge0 x≥0的部分是一个单调递增的函数(在定义域上单调非减)，定义域和值域是 [ 0 , + ∞ ) → [ 0 , 1 ) [0,+\infty) \to [0,1) [0,+∞)→[0,1)，在 p ( x ) p(x) p(x)大的地方它增长快，反之亦然。

因为它是唯一映射的（在>0的部分，接下来我们只考虑这一部分），所以它的反函数可以表示为 F − 1 ( a ) ， a ∈ [ 0 , 1 ) , 值 域 为 [ 0 , + ∞ ) F^{-1}(a)，a\in [0,1), 值域为[0,+\infty) F−1(a)，a∈[0,1),值域为[0,+∞)。(上图的左图就是 F ( x ) F(x) F(x)的图)

因为 F F F单调递增，所以 F − 1 F^{-1} F−1也是单调递增的：
x ≤ y ⇒ F ( x ) ≤ F ( y ) a ≤ b ⇒ F − 1 ( a ) ≤ F − 1 ( b ) (3) \begin{aligned} x\le y &amp;\Rightarrow F(x) \le F(y)\\ a\le b &amp;\Rightarrow F^{-1}(a) \le F^{-1}(b)\qquad \text{(3)} \end{aligned} x≤ya≤b​⇒F(x)≤F(y)⇒F−1(a)≤F−1(b)(3)​

利用反函数的定义，我们有：
F − 1 ( a ) &lt; x , i f f    a &lt; F ( x ) (4) F^{-1}(a)&lt;x,iff\ \ a&lt;F(x)\qquad \text{(4)} F−1(a)<x,iff  a<F(x)(4)

接下来，我们定义一下 [ 0 , 1 ] [0,1] [0,1]均匀分布的CDF,这个很好理解：
P ( a ≤ x ) = H ( x ) = { 1 , x ≥ 1 x , 0 ≤ x ≤ 1 0 , x &lt; 0 (5) P(a\le x) = H(x) = \begin{cases} 1 &amp;,x \ge 1\\ x &amp;,0\le x\le1\\ 0&amp;, x\lt0 \end{cases}\qquad \text{(5)} P(a≤x)=H(x)=⎩⎪⎨⎪⎧​1x0​,x≥1,0≤x≤1,x<0​(5)

根据公式（4）和公式（5）,有：
P ( F − 1 ( a ) ≤ x ) = P ( a ≤ F ( x ) ) = H ( F ( x ) ) (6) P(F^{-1}(a) \le x) = P(a \le F(x)) = H(F(x))\qquad \text{(6)} P(F−1(a)≤x)=P(a≤F(x))=H(F(x))(6)

因为 F ( x ) F(x) F(x)的值域 [ 0 , 1 ) [0,1) [0,1)，根据公式（5），公式（6）可以改写为：
P ( F − 1 ( a ) ≤ x ) = H ( F ( x ) ) = F ( x ) (7) P(F^{-1}(a) \le x) = H(F(x)) = F(x)\qquad \text{(7)} P(F−1(a)≤x)=H(F(x))=F(x)(7)

据 F ( x ) F(x) F(x)的定义，它是exp分布的CDF，所以公式（7）的意思是 F − 1 ( a ) F^{-1}(a) F−1(a)符合exp分布，我们通过 F F F的反函数将一个0到1均匀分布的随机数转换成了符合exp分布的随机数，注意，以上推导对于CDF可逆的分布都是一样的。对于exp来说，它的反函数的形式是：
F e x p − 1 ( a ) = − 1 λ ∗ l o g ( 1 − a ) (8) F^{-1}_{exp}(a) = -\frac{1}{\lambda}*log(1-a)\qquad \text{(8)} Fexp−1​(a)=−λ1​∗log(1−a)(8)

具体的映射关系可以看下图(a)，我们从 y 轴 0-1 的均匀分布样本（绿色）映射得到了服从指数分布的样本（红色）。

我们写一点代码来看看效果:

def sampleExp(Lambda = 2,maxCnt = 50000):
    ys = []
    standardXaxis = []
    standardExp = []
    for i in range(maxCnt):
        u = np.random.random()
        y = -1/Lambda*np.log(1-u) #F-1(X)
        ys.append(y)
    for i in range(1000):
        t = Lambda * np.exp(-Lambda*i/100)
        standardXaxis.append(i/100)
        standardExp.append(t)
    plt.plot(standardXaxis,standardExp,'r')
    plt.hist(ys,1000,normed=True)
    plt.show()

sampleExp()

最后绘制出来的直方图可以看出来就是 exp 分布图，见上图(b)。可以看到随着采样数量的变多，概率直方图和真实的 CDF 就越接近。

以上就是逆采样的过程。我们的结论是：因为CDF是单调函数（累积的概率只能越来越大，直到为1），因此，只要某分布的CDF可逆，那么就可以通过均匀分布来采样服从该分布的样本集。

二、拒绝采样(Reject Sampling)

2.1 原理

对于那些复杂的不可求逆的分布，拒绝采样就是针对此类复杂问题的一种随机采样方法。

我们以求圆周率 π \pi π的例子入手，讲解拒绝采样的思想。通过采样的方法来计算 π \pi π值，也就是在一个 1 × 1 1×1 1×1的范围内随机采样一个点，如果它到原点的距离小于1,则说明它在1/4圆内，则接受它，最后通过接受的占比来计算1/4圆形的面积，从而根据公式反算出预估 π ^ \hat{\pi} π^的值，随着采样点的增多，最后的结果 π ^ \hat{\pi} π^会越精准。

上面这个例子里说明一个问题，我们想求一个空间里均匀分布的集合面积，可以尝试在更大范围内按照均匀分布随机采样，如果采样点在集合中，则接受，否则拒绝。最后的接受概率就是集合在”更大范围“的面积占比。

接下来，我们来形式化地说明拒绝采样。

给定一个概率分布 p ( z ) = 1 Z p p ~ ( z ) p(z)=\frac{1}{Z_p} \tilde{p}(z) p(z)=Zp​1​p~​(z)，其中， p ~ ( z ) \tilde{p}(z) p~​(z)已知， Z p Z_p Zp​为归一化常数，未知。要对该分布 p ( z ) p(z) p(z)进行拒绝采样，首先需要借用一个简单的参考分布(proposal distribution)，记为 q ( x ) q(x) q(x)，该分布的采样易于实现，如均匀分布、高斯分布。然后引入常数 k k k，使得对所有的的 z z z，满足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)，如下图所示，红色的曲线为 p ~ ( z ) \tilde{p}(z) p~​(z)，蓝色的曲线为 k q ( z ) kq(z) kq(z)。

在每次采样中，首先从 q ( z ) q(z) q(z)采样一个数值 z 0 z_0 z0​，然后在区间 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]进行均匀采样，得到 u 0 u_0 u0​。如果 u 0 &lt; p ~ ( z 0 ) u_0&lt;\tilde{p}(z_0) u0​<p~​(z0​)，则保留该采样值，否则舍弃该采样值。最后得到的数据就是对 p ( z ) {p}(z) p(z)分布的一个近似采样。 结合图，直观来理解上述的过程：在 x = z 0 x=z_0 x=z0​这条线上，从 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]均匀采样中一个值，如果这个值小于 p ~ ( z 0 ) \tilde{p}(z_0) p~​(z0​)，即这个均匀采样的这个值落在了 p ~ ( z 0 ) \tilde{p}(z_0) p~​(z0​)的下方，我们就接受 z 0 z_0 z0​这个采样值。

我们知道，每次采样的接受概率计算如下：
p ( a c c e p t ) = ∫ p ~ ( z ) k q ( z ) q ( z ) d z = 1 k ∫ p ~ ( z ) d z p(accept)=\displaystyle\int{\dfrac{\tilde{p}(z)}{kq(z)}}q(z)dz = \dfrac{1}{k} \int{\tilde{p}(z)}dz p(accept)=∫kq(z)p~​(z)​q(z)dz=k1​∫p~​(z)dz

所以，为了提高接受概率，防止舍弃过多的采样值而导致采样效率低下， k k k的选取应该在满足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)的基础上尽可能小。

拒绝采样问题可以这样理解， p ~ ( z ) \tilde{p}(z) p~​(z)与 x x x轴之间的区域为要估计的问题，类似于上面提到的圆形区域， k q ( z ) kq(z) kq(z)与 x x x轴之间的区域为参考区域，类似于上面提到的正方形。由于 k q ( z ) kq(z) kq(z)与 x x x轴之间的区域面积为 k k k(原来的概率密度函数的面积为1，现在扩大了 k k k倍，故 k × 1 k\times1 k×1)，所以， p ~ ( z ) \tilde{p}(z) p~​(z)与 x x x轴之间的区域面积除以 k k k即为对 p ( z ) p(z) p(z)的估计。在每一个采样点，以 [ 0 , k q ( z 0 ) ] [0,kq(z_0)] [0,kq(z0​)]为界限，落在 p ~ ( z ) \tilde{p}(z) p~​(z)曲线以下的点就是服从 p ( z ) p(z) p(z)分布的点。

2.2 实验

我们写一段代码来看看拒绝采样是如何对复杂分布进行采样的。假设，我们要采样的分布是：
p ~ ( z ) = 0.3 e x p ( − ( z − 0.3 ) 2 ) + 0.7 e x p ( − ( z − 2 ) 2 / 0.3 ) \tilde{p}(z)=0.3exp(-(z-0.3)^2)+0.7exp(-(z-2)^2/0.3) p~​(z)=0.3exp(−(z−0.3)2)+0.7exp(−(z−2)2/0.3)

其归一化常数 Z p ≈ 1.2113 Z_p \approx 1.2113 Zp​≈1.2113,参考分布为高斯分布 q ( z ) = G a s s i a n ( 1.4 , 1.2 ) q(z)=Gassian(1.4,1.2) q(z)=Gassian(1.4,1.2)，其中均值和方差是经过计算和尝试得到的，以满足 k q ( z ) ≥ p ~ ( z ) kq(z)\geq\tilde{p}(z) kq(z)≥p~​(z)。python代码如下：

import numpy as np
import matplotlib.pyplot as plt

def f1(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,f1(x),color = "red")

size = int(1e+07)
sigma = 1.2
z = np.random.normal(loc = 1.4,scale = sigma, size = size)
qz = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*(z-1.4)**2/sigma**2)
k = 2.5
#z = np.random.uniform(low = -4, high = 6, size = size)
#qz = 0.1
#k = 10
u = np.random.uniform(low = 0, high = k*qz, size = size)

pz =  0.3*np.exp(-(z-0.3)**2) + 0.7* np.exp(-(z-2.)**2/0.3)
sample = z[pz >= u]
plt.hist(sample,bins=150, normed=True, edgecolor='black')
plt.show()

得到的采样分布如下图：

可以看到，采样结果完全符合原分布。另外如果我们把参考分布换为均匀分布(代码中z,q,k换为注释部分)，仍然得到较好的采样结果，如下图:

可以看到，采样结果也是完全符合原分布。

以上的实验说明了拒绝采样的有效性。

2.3 注意事项

通过在2.1节的分析，我们知道：必须知道复杂分布的概率密度函数PDF，才可以进行拒绝采样。然而，在现实情况中：

1）对于一些二维分布 p ( x , y ) p(x,y) p(x,y)，有时候我们只能得到条件分布 p ( x ∣ y ) p(x|y) p(x∣y)和 p ( y ∣ x ) p(y|x) p(y∣x)，却很难得到二维分布的概率密度函数 p ( x , y ) p(x,y) p(x,y)的一般形式，这时我们无法用拒绝采样得到其样本集。

2）对于一些高维的复杂非常见分布 p ( x 1 , x 2 , . . . , x n ) p(x_1,x_2,...,x_n) p(x1​,x2​,...,xn​)，我们要找到一个合适的 q ( x ) q(x) q(x)和 k k k非常困难。

因此，实际上，我们仍然要找到一种方法可以解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。马尔科夫链有能力帮助找到这些复杂概率分布的对应的采样样本集。而这也是MCMC采样的基础，我们回头会讲解。

参考文献

【1】拒绝采样(reject sampling)原理详解

【2】采样方法（一）

【3】蒙特卡洛采样之拒绝采样（Reject Sampling）

【4】一文了解采样方法