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

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

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

一、逆采样(Inverse Sampling)

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

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

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

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

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

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

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

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

接下来，我们定义一下 $[0,1]$ 均匀分布的CDF,这个很好理解：
$P(a\le x) = H(x) = \begin{cases} 1 &,x \ge 1\\ x &,0\le x\le1\\ 0&, x\lt0 \end{cases}\qquad \text{(5)}$

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

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

据 $F(x)$ 的定义，它是exp分布的CDF，所以公式（7）的意思是 $F^{-1}(a)$ 符合exp分布，我们通过 $F$ 的反函数将一个0到1均匀分布的随机数转换成了符合exp分布的随机数，注意，以上推导对于CDF可逆的分布都是一样的。对于exp来说，它的反函数的形式是：
$F^{-1}_{exp}(a) = -\frac{1}{\lambda}*log(1-a)\qquad \text{(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/4圆内，则接受它，最后通过接受的占比来计算1/4圆形的面积，从而根据公式反算出预估 $\hat{\pi}$ 的值，随着采样点的增多，最后的结果 $\hat{\pi}$ 会越精准。

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

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

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

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

我们知道，每次采样的接受概率计算如下：
$p(accept)=\displaystyle\int{\dfrac{\tilde{p}(z)}{kq(z)}}q(z)dz = \dfrac{1}{k} \int{\tilde{p}(z)}dz$

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

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

2.2 实验

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

其归一化常数 $Z_p \approx 1.2113$ ,参考分布为高斯分布 $q(z)=Gassian(1.4,1.2)$ ，其中均值和方差是经过计算和尝试得到的，以满足 $kq(z)\geq\tilde{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()