R-概率统计与模拟（四）拒绝抽样

本文介绍了拒绝抽样（Reject Sampling）。

前文《R-概率统计与模拟（三）变换均匀分布对特定分布进行抽样》介绍了通过“变换均匀分布”来对特定分布进行抽样的方法，但是该方法需要知道累积分布的解析表达式及其反函数，所以有一定的限制。

其实，我们最常接触的还是 $\text{p.d.f.}$，根据$\text{p.d.f.}$抽样往往更直接。比如，均匀分布的$\text{p.d.f.}$就很简单，对其抽样也很方便。但是，对于一些复杂的$\text{p.d.f.}$，直接抽样是不可取的，需要用其它方法。拒绝抽样就是其中一种。

下文分为五个部分：

1. 拒绝抽样简介
2. 利用拒绝抽样方法对$\text{Gamma}$分布抽样
3. 如何利用rand7对rand10抽样
4. 蓄水池抽样算法
5. 对拒绝抽样的补充说明

拒绝抽样简介

所谓拒绝抽样，是针对一个已知$\text{p.d.f.}$的概率分布进行抽样的方法。具体步骤如下：

1. 假设我们要对一个分布进行随机抽样，已知其$\text{p.d.f.}$是$p(x)$。
2. $p(x)$可能比较复杂以至于直接抽样有困难，我们可以选择一个容易抽样的分布，假设该分布$\text{p.d.f.}$是$q(x)$，再选择一个足够大的数$M$，使得$Mq(x)$始终比$p(x)$大。
3. 根据$q(x)$抽取一个值$x_i$。
4. 从$[0,1]$均匀分布中抽取一个数$u_i$，如果$$u_i<\frac{p(x_i)}{Mq(x_i)},$$ 那么接受$x_i$这个采样值，否则拒绝它。
5. 重复第3步和第4步，知道采样的数量达到预定要求。

通过这些步骤获取到的采样值就是符合$p(x)$分布的。需要注意的是$M$的取值，如果$M$过大，会导致抽样过程中被拒绝的点增多，降低采样点的接受率。事实上，我们可以证明采样点的接受率是$\frac{1}{M}$，所以$M$的取值应该在保证$Mq(x)>p(x),\forall x$的前提下尽可能得小。（相关的证明在文末。）

像上图一样，$M$的取值要保证$Mq(x)>p(x),\forall x$。

利用拒绝抽样方法对 $\text{Gamma}$分布抽样

上面的步骤比较抽象，我们举一个实际的例子来说明：我们怎么通过拒绝抽样的方法对$\text{Gamma}(x,\alpha ,1)$进行抽样？
（选这个例子是因为它也用到了我们在前文提到的“变换均匀分布”的方法。）

我们对照上面的步骤一步一步来说：

1. 我们知道，目标分布$\text{Gamma}(x,\alpha ,1)$的$\text{p.d.f.}$是：
   $$p(x)=\frac{e^{-x}{x}^{\alpha -1}}{\Gamma (\alpha )}$$
2. 显然$p(x)$比较复杂，不易直接采样。我们选择一个较容易采样的分布，其$\text{p.d.f.}$是：
   q ( x ) = λ α λ x λ − 1 ( α λ + x λ ) 2 , q(x) = \frac{\lambda \alpha^{\lambda} x^{\lambda-1}}{(\alpha^\lambda + x^\lambda)^2}, q(x)=(αλ+xλ)