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

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

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

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

下文分为五个部分：

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

拒绝抽样简介

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

假设我们要对一个分布进行随机抽样，已知其p.d.f.\text{p.d.f.}p.d.f.是p(x)p(x)p(x)。

p(x)p(x)p(x)可能比较复杂以至于直接抽样有困难，我们可以选择一个容易抽样的分布，假设该分布p.d.f.\text{p.d.f.}p.d.f.是q(x)q(x)q(x)，再选择一个足够大的数MMM，使得Mq(x)Mq(x)Mq(x)始终比p(x)p(x)p(x)大。

根据q(x)q(x)q(x)抽取一个值xix_ixi。

从[0,1][0,1][0,1]均匀分布中抽取一个数uiu_iui，如果ui<p(xi)Mq(xi),u_i < \frac{p(x_i)}{Mq(x_i)},ui<Mq(xi)p(xi), 那么接受xix_ixi这个采样值，否则拒绝它。

重复第3步和第4步，知道采样的数量达到预定要求。

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

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

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

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

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

我们知道，目标分布Gamma(x,α,1)\text{Gamma}(x, \alpha, 1)Gamma(x,α,1)的p.d.f.\text{p.d.f.}p.d.f.是：
p(x)=e−xxα−1Γ(α)p(x) = \frac{e^{-x}x^{\alpha -1}}{\Gamma(\alpha)}p(x)=Γ(α)e−xxα−1

显然p(x)p(x)p(x)比较复杂，不易直接采样。我们选择一个较容易采样的分布，其p.d.f.\text{p.d.f.}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λ)2λαλxλ−1,
我们选择MMM：
M=4e−αααλΓ(α),M = \frac{4e^{-\alpha} \alpha^\alpha}{\lambda \Gamma(\alpha)},M=λΓ(α)4e−ααα,
所以：
Mq(x)=4e−ααλ+αxλ−1Γ(α)(αλ+xλ)2.Mq(x) = \frac{4e^{-\alpha} \alpha^{\lambda + \alpha} x^{\lambda-1}}{\Gamma(\alpha) (\alpha^\lambda + x^\lambda)^2}.Mq(x)=Γ(α)(αλ+xλ)24e−ααλ+αxλ−1.
可以证明，当α>1\alpha > 1α>1 并且 1≤λ≤2α−11 \le \lambda \le \sqrt{2\alpha - 1}1≤λ≤2α−1 时，对所有xxx都有p(x)<Mq(x)p(x) < Mq(x)p(x)<Mq(x)。并且可以看出，当λ\lambdaλ取值不同时，q(x)q(x)q(x)、MMM 以及 Mq(x)Mq(x)Mq(x) 都是在变化的，也会影响到接受率。

根据q(x)q(x)q(x)抽取一个值xix_ixi。看起来q(x)q(x)q(x)也是蛮复杂的，怎么对其抽样呢？我们可以通过变换均匀分布的方法来实现：
我们从[0,1][0,1][0,1]的均匀分布中抽取一个值rir_iri，然后令
xi=α(ri1−ri)1/λx_i = \alpha \left(\frac{r_i}{1-r_i} \right)^{1 / \lambda} xi=α(1−riri)1/λ
这种通过变换均匀分布的方式得到的xix_ixi等价于从q(x)q(x)q(x)中抽取xix_ixi。

从[0,1][0,1][0,1]均匀分布中抽取一个数uiu_iui，如果ui<p(xi)Mq(xi),u_i < \frac{p(x_i)}{Mq(x_i)},ui<Mq(xi)p(xi), 那么接受xix_ixi这个采样值，否则拒绝它。

重复第3步和第4步，知道采样的数量达到预定要求。

上文已经提及，λ\lambdaλ取值不同会影响MMM值，从而影响到接受率。为了验证这一点，笔者写了一段代码来做一个实验：
假设我们让α=5\alpha=5α=5，再分别让λ=3\lambda=3λ=3以及λ=1\lambda=1λ=1来看看λ\lambdaλ取值不同对接受率的影响。
代码如下：

# Step1. p(x)函数
px <- function(x, a) dgamma(x, shape=a, scale=1)  # Step2. q(x)函数
qx <- function(x, a, l) a ^ l * l * x ^ (l - 1) / ((a ^ l + x ^ l ) ^ 2)
getM <- function(a, l) 4 * exp(-a) * a ^ a / (gamma(a) * l)  # 计算M值
mqx <- function(x, a, l, M) M * qx(x, a, l)     # Mq(x)函数# Step3. 变换均匀分布以便对q(x)抽样的函数
transUnif <- function(x, a, l) a * (x / (1 - x)) ^ (1 / l)# 单纯为了画p(x)的理论曲线用的函数
px.1 <- function(x) px(x, alpha)
# 单纯为了画lambda=lambda.1时Mq(x)的理论曲线用的函数
mqx.1 <- function(x) mqx(x, alpha, lambda.1, M.1)
# 单纯为了画lambda=lambda.2时Mq(x)的理论曲线用的函数
mqx.2 <- function(x) mqx(x, alpha, lambda.2, M.2)ogamma <- function(a, l, M) {while (T) {x <- runif(1, 0, 1)y <- transUnif(x, a, l)   # 变换均匀分布实现对q(x)的抽样u <- runif(1, 0, 1) if (u < px(y, a) / mqx(y, a, l, M))   # Step4. 是否接受采样点breaknfail <<- nfail + 1   # 记录被拒绝的次数}y
}sgamma <- function(n, a, l, M) {set.seed(123)replicate(n, ogamma(a, l, M))
}N <- 100000
alpha <- 5# 第一种lambda值的效果
lambda.1 <- sqrt(2 * alpha - 1)
M.1 <- getM(alpha, lambda.1)
nfail <- 0
res.1 <- sgamma(N, alpha, lambda.1, M.1)
nfail / (nfail + N)   # 计算模拟的拒绝率，14.38%
1 - 1 / M.1         # 计算理论的拒绝率，14.51%
plot(density(res.1), xlim=c(0, 20), col="red", xlab="x",main=paste0("Reject Sampling for Gamma(x, ", alpha, ", 1)\nwith lambda=", lambda.1))
curve(px.1, 0, 20, col="blue", add=T, lty=2)
curve(mqx.1, 0, 20, col="black", add=T, lty=3)
legend("topright", legend=c("simulative", "theoretical p(x)", "theoretical Mq(x)"),col=c("red", "blue", "black"), lty=c(1,2,3), bty="n")# 第二种lambda值的效果
lambda.2 <- 1
M.2 <- getM(alpha, lambda.2)
nfail <- 0
res.2 <- sgamma(N, alpha, lambda.2, M.2)
nfail / (nfail + N)   # 计算模拟的拒绝率，71.54%
1 - 1 / M.2          # 计算理论的拒绝率，71.50%
plot(density(res.2), xlim=c(0, 20), ylim=c(0, 0.7), col="red", xlab="x",main=paste0("Reject Sampling for Gamma(x, ", alpha, ", 1)\nwith lambda=", lambda.2))
curve(px.1, 0, 20, col="blue", add=T, lty=2)
curve(mqx.2, 0, 20, col="black", add=T, lty=3)
legend("topright", legend=c("simulative", "theoretical p(x)", "theoretical Mq(x)"),col=c("red", "blue", "black"), lty=c(1,2,3), bty="n")

当α=5,λ=3\alpha=5, \lambda=3α=5,λ=3时，代码模拟出的接受率是85.62%85.62\%85.62%，理论上的接受率是1M=85.49%\frac{1}{M}=85.49\%M1=85.49%，很接近。采样点的具体分布如下：

当α=5,λ=1\alpha=5, \lambda=1α=5,λ=1时，代码模拟出的接受率是28.46%28.46\%28.46%，理论上的接受率是1M=28.5%\frac{1}{M}=28.5\%M1=28.5%，也很接近。采样点的具体分布如下：

我们可以看出，当λ\lambdaλ分别取333和111时，接受率从85.49%85.49\%85.49%下降到了28.5%28.5\%28.5%。所以，λ\lambdaλ对接受率的影响是非常大的。

如何利用rand7对rand10抽样

这来源于一个常见的算法题：如果已知一个rand7()函数，它可以从1到7这7个数字中随机地、等概率地抽取一个数，那么如何利用rand7()函数实现rand10()函数呢？所谓rand10()函数，就是这个函数可以从1到10这10个数字中随机地、等概率地抽取一个数。

这个问题的一个常见解答步骤是：

首先利用[rand7()−1]×7+rand7()[rand7() - 1] \times 7 + rand7()[rand7()−1]×7+rand7()这个表达式实现从1-49这49个数中随机而等概率地抽取一个数xxx；

如果x≤40x \le 40x≤40，那么接受xxx，并将x%10+1x \% 10 + 1x%10+1的值作为一个采样点；否则拒绝xxx；那么这样得到的一系列采样点就是符合预期要求的。

这个解答步骤为什么是正确的，证明也很简单，和文末的证明过程类似，所以在此略过。那为什么要介绍这个题目呢？因为我们看到抽样过程中也涉及到了“接受-拒绝”的过程，所以笔者认为这可以看作文初所述的拒绝抽样过程的一种变型。

这个问题的代码也很简单：

rand7 <- function() {sample(7, 1)    # 从1-7中随机抽取一个整数
}rand10 <- function() {while (T) {x <- (rand7() - 1) * 7 + rand7()   # 1-49 uniformly.if (x <= 40)break}x %% 10 + 1    # 一个服从均匀分布的rand10的采样值
}N <- 100000  # 样本大小
set.seed(123)
res <- replicate(N, rand10())
hist(res, prob=T, breaks = 0:10, ylim=c(0, 0.15), xlab="x",main="Simulation of rand10")
abline(h=0.1, col="red", lty=2)

结果如下：

从上图中可以看出，rand10()的抽样结果符合均匀分布，达到预期。

蓄水池抽样算法

如同上面rand10()的题目一样，蓄水池算法的实现过程中也涉及到了“接受-拒绝”过程，所以在此一并介绍：

所谓蓄水池抽样算法，主要应用于如下问题：即我们要从1,2,…,N1,2,\ldots,N1,2,…,N这NNN个样本中随机抽取mmm个样本，NNN非常大或者未知，且m<<Nm << Nm<<N。
蓄水池抽样算法的步骤是：

首先将1,2,…,m1,2,\ldots,m1,2,…,m这mmm个样本选进来；

从i=m+1i=m+1i=m+1开始，从[1,i][1, i][1,i]中随机等概率抽取一个数rrr，如果1≤r≤m1 \le r \le m1≤r≤m，那么用第iii个样本替代第rrr个样本。否则不做处理

i=i+1i=i+1i=i+1。

重复第2步和第3步，直至处理完全部的样本。

其证明过程也很简单，用归纳法即可。笔者曾经写过一篇文章介绍了该算法在测序数据reads抽取方面的应用：《算法（二）蓄水池抽样算法快速随机抽取reads》。

对拒绝抽样的补充说明

为什么按照拒绝抽样的过程进行抽样所得到的采样点就是符合p(x)p(x)p(x)的分布的？

证明：
首先我们给出一些记号：假设我们将目标分布p(x)p(x)p(x)中的点记为X^\hat{X}X^，从q(x)q(x)q(x)中抽取的任意一个值记为XXX；并且让事件AAA代表XXX这个值被接受了。如果当X=xiX=x_iX=xi时被接受了，那么就是X^=X=xi\hat{X}=X=x_iX^=X=xi。

有了这些记号，我们可以很容易地知道：
Pr⁡(A∣X=xi)=Pr⁡[ui<p(xi)Mq(xi)]=p(xi)Mq(xi)\Pr(A|X=x_i)=\Pr \left[u_i < \frac{p(x_i)}{Mq(x_i)} \right] = \frac{p(x_i)}{Mq(x_i)}Pr(A∣X=xi)=Pr[ui<Mq(xi)p(xi)]=Mq(xi)p(xi)

利用全概率公式，我们知道：
Pr⁡(A)=∫Pr⁡(A∣X=xi)q(xi)dxi=∫p(xi)Mq(xi)q(xi)dxi=1M∫p(xi)dxi=1M\begin{aligned} \Pr(A) & = \int \Pr(A|X=x_i) \, q(x_i) {\rm d}x_i \\ & = \int \frac{p(x_i)}{Mq(x_i)}q(x_i) {\rm d}x_i \\ & = \frac{1}{M} \int p(x_i) {\rm d}x_i \\ & = \frac{1}{M} \end{aligned}Pr(A)=∫Pr(A∣X=xi)q(xi)dxi=∫Mq(xi)p(xi)q(xi)dxi=M1∫p(xi)dxi=M1

最终，我们想要知道的X^\hat{X}X^的分布：
Pr⁡(X^≤xi)=Pr⁡(X≤xi∣A)\Pr(\hat{X} \le x_i) = \Pr(X \le x_i|A)Pr(X^≤xi)=Pr(X≤xi∣A)

利用贝叶斯公式，我们有：
Pr⁡(X^≤xi)=Pr⁡(X≤xi∣A)=Pr⁡(X≤xi,A)Pr⁡(A)=∫xiq(t)p(t)Mq(t)dtPr⁡(A)=1M∫xip(t)dt1M=∫xip(t)dt\begin{aligned} \Pr(\hat{X} \le x_i) & = \Pr(X \le x_i|A) \\ & = \frac{\Pr(X \le x_i, A)}{\Pr(A)} \\ & = \frac{\displaystyle \int^{x_i} q(t)\frac{p(t)}{Mq(t)} {\rm d}t}{\Pr(A)} \\ & = \cfrac{\cfrac{1}{M} \displaystyle \int^{x_i} p(t) {\rm d}t}{\cfrac{1}{M}} \\ & = \int^{x_i} p(t) {\rm d}t \end{aligned}Pr(X^≤xi)=Pr(X≤xi∣A)=Pr(A)Pr(X≤xi,A)=Pr(A)∫xiq(t)Mq(t)p(t)dt=M1M1∫xip(t)dt=∫xip(t)dt

所以可以证明，依照拒绝抽样过程得到的采样点符合p(x)p(x)p(x)分布。

拒绝采样的接受率如何计算？

从上面的证明过程中可以看出，接受率就是Pr⁡(A)\Pr(A)Pr(A)，等于1M\frac{1}{M}M1。上述证明过程是笔者按照自己的理解写的，如果有错误，请朋友们指正。

小结

本文介绍了拒绝采样，并以Gamma\text{Gamma}Gamma分布、rand10以及蓄水池抽样算法等三个例子加以说明。拒绝采样不同于变换均匀分布的方法，它是直接根据p.d.f.\text{p.d.f.}p.d.f.进行抽样。其缺点是要找到合适的q(x)q(x)q(x)和MMM值并不容易，尤其是在高维的情况下，MMM值往往偏大，从而导致拒绝率很高。

参考

《生物序列分析》第11章

（公众号：生信了）