拒绝抽样

拒绝抽样（rejection sampling，又称接受拒绝抽样）从一个容易抽样的提议分布产生候选，再按目标密度与包络密度的比例决定是否接受。只要包络条件正确，接受的结果就服从目标分布；各次候选和接受随机数独立时，依次接受的样本也独立同分布。

例如，希望在0至1之间抽取随机数，让较大的值更常出现，具体密度为 $f(x)=2x$。均匀抽一个数并不符合要求；但先均匀抽候选x，再以概率x接受它，就能得到所需分布。下面从这个三角形密度出发，证明接受规则，并计算它要花多少次尝试。

单位正方形中的接受区域

独立抽取 $X\sim \mathrm{Uniform}(0,1)$、$U\sim \mathrm{Uniform}(0,1)$，把点画在单位正方形内。保留 $U\le X$ 的点，丢弃其余点。对角线右下侧面积为二分之一，所以候选平均有一半被接受。

横坐标X是候选值，纵坐标U是独立的接受随机数。候选对在正方形内均匀分布；接受后取点的横坐标。

在很窄的区间 $[x,x+\mathrm{d}x]$，候选出现的概率约为 $\mathrm{d}x$，接受的比例约为x，所以“落在该区间而且接受”的概率约为 $x\mathrm{d}x$。除以总接受概率二分之一，得到条件密度 $2x$。这里增加大数出现机会的过程，是以正确的比例选择候选，而不是把均匀样本简单排序或重复较大的数。

包络条件与完整证明

设目标密度f、提议密度g相对于同一参考测度定义，均归一化为一。要求有可验证的有限常数 $M\ge 1$，使 $$f(x)\le Mg(x)\text{几乎处处成立}.$$ 因此在f为正的区域，g也必须为正。独立抽取候选X∼g及U∼Uniform(0,1)，在 $U\le f(X)/(Mg(X))$ 时接受。g为零的点不会被提议，对这部分零概率点不必计算比值。

记接受事件为A。由条件期望， $$Pr(A)=\int g(x)\frac{f(x)}{Mg(x)}\mathrm{d}x=\frac{1}{M}.$$ 对任意可测集合B，同样有 $$Pr(X\in B,A)=\frac{1}{M}{\int }_{B}f(x)\mathrm{d}x.$$ 两式相除，得到 $Pr(X\in B\mid A)={\int }_{B}f(x)\mathrm{d}x$，这就证明了接受后分布正确。每次尝试使用全新的独立候选对；一次接受之后，后面的候选对与此前独立，所以重新开始的接受样本仍独立且服从同一个f。

三角形例中g=1，f=2x，可取M=2，于是接受阈值正好是x。接受后 $F(x)=x^2$（0≤x≤1），故 $$𝔼X={\int }_0^{1}x(2x)\mathrm{d}x=\frac{2}{3},Pr(X>1/2)=1-(1/2{)}^2=\frac{3}{4}.$$ 均匀候选在二分之一右边的概率原来只有二分之一；保留规则把接受样本中的该比例提高为四分之三。

尝试次数与未归一化密度

一次成功概率为 $a=1/M$。取得一个接受样本所需的尝试次数T满足 $Pr(T=t)=(1-a{)}^{t-1}a$。从第一次尝试分解，也能直接列出 $$𝔼T=1+(1-a)𝔼T,𝔼T=1/a=M.$$ 这给出了平均代价，不是每个样本的固定尝试次数。取得n个接受样本，期望共需nM次尝试；在候选生成和密度计算的单位成本固定时，总期望时间为 $O(nM)$。

有时只知道未归一化的 $h(x)=Zf(x)$，其中 $0<Z<\mathrm{\infty }$ 未知。若能找到c使 $h(x)\le cg(x)$，仍可按 $U\le h(X)/(cg(X))$ 接受。接受概率变为 $Z/c$，接受后密度为 $h/Z$；实现不需要知道Z。构造有效包络c仍不可省略。

把上界估小会改变分布

对三角形例，若错误取M=1，计算出的阈值2x在x>1/2时超过一。程序将这些候选全部接受，实际接受概率相当于 $min(2x,1)$。总接受率成为 $${\int }_0^{1/2}2x\mathrm{d}x+{\int }_{1/2}^{1}1\mathrm{d}x=\frac{3}{4},$$ 因此输出密度为 $$f_{\mathrm{w}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{g}}(x)=\{\begin{array}{cc}(8/3)x,& 0<x<1/2,\\ 4/3,& 1/2<x<1.\end{array}$$ 它的右半段成了平台，均值为 $11/18$，不再是 $2/3$。较高的接受率来自抽错分布，不能当作效率改进。

青线为目标2x，金色虚线为错误取M=1后的实际输出密度；右半区错误密度为4/3。

三角形密度的Python实现

下面实现已证明的U≤X规则，返回n个接受样本和总尝试次数。它不接收任意函数，因而不暗示已经替其他密度验证了包络条件。

Python 3

import random

def triangle_rejection(n, seed=26):
    if type(n) is not int or n < 0:
        raise ValueError("n must be a nonnegative integer")
    rng = random.Random(seed)
    samples, attempts = [], 0
    while len(samples) < n:
        x = rng.random()
        u = rng.random()
        attempts += 1
        if u <= x:
            samples.append(x)
    return samples, attempts

values, attempts = triangle_rejection(10000)
print("sample mean:", sum(values) / len(values))
print("accepted fraction:", len(values) / attempts)

循环不变量是：列表只包含已经通过接受检验的候选，attempts计入所有尝试。对于n>0，它在理想独立均匀随机数模型下几乎必然停止，但没有固定的最坏尝试次数；期望时间为O(n)，返回列表占O(n)空间。n=0返回空列表及零次尝试，示例中的均值表达式只适用于非空输出。伪随机数与浮点计算是这一概率模型的数值实现，模拟均值应在三分之二附近波动。

变式。 仍使用均匀提议，却把有效包络改为M=4，接受阈值成为x/2。输出密度仍为2x，但接受率从二分之一降至四分之一，平均尝试次数由二变四。若目标恰好也是均匀分布，取M=1便可接受全部候选。

高维问题的困难通常在于有效M可能很大。例如在d维边长二的立方体中均匀提议，目标为内接单位球上的均匀分布，接受率等于球体积除以立方体体积；大量候选可能被丢弃。若提议尾部比目标衰减得更快，密度比还可能没有有限上界。此时应改变提议，或考虑重要性抽样、马尔可夫链蒙特卡洛等方法，并重新核验其条件。

参考资料

• Gautam Iyer：Basic Sampling Algorithms，“Rejection sampling”，算法14、命题16与引理17。