R-概率统计与模拟(四)拒绝抽样

本文详细介绍了拒绝抽样方法,包括其基本原理、如何利用拒绝抽样对Gamma分布抽样,以及结合rand7实现rand10抽样和蓄水池抽样算法。文中通过实例和代码展示了拒绝抽样的过程,并讨论了如何计算接受率及其影响因素。
摘要由CSDN通过智能技术生成

本文介绍了拒绝抽样(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.,直接抽样是不可取的,需要用其它方法。拒绝抽样就是其中一种。

下文分为五个部分:

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

拒绝抽样简介

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

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

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

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

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

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

  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 − x x α − 1 Γ ( α ) p(x) = \frac{e^{-x}x^{\alpha -1}}{\Gamma(\alpha)} p(x)=Γ(α)exxα1
  2. 显然 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,
    我们选择 M M M
    M = 4 e − α α α λ Γ ( α ) , M = \frac{4e^{-\alpha} \alpha^\alpha}{\lambda \Gamma(\alpha)}, M=λΓ(α)4eααα,
    所以:
    M q ( x ) = 4 e − α α λ + α 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 α − 1 1 \le \lambda \le \sqrt{2\alpha - 1} 1λ2α1 时,对所有 x x x都有 p ( x ) < M q ( x ) p(x) < Mq(x) p(x)<Mq(x)。并且可以看出,当 λ \lambda λ取值不同时, q ( x ) q(x) q(x) M M M 以及 M q ( x ) Mq(x) Mq(x) 都是在变化的,也会影响到接受率。
  3. 根据 q ( x ) q(x) q(x)抽取一个值 x i x_i xi。看起来 q ( x ) q(x) q(x)也是蛮复杂的,怎么对其抽样呢?我们可以通过变换均匀分布的方法来实现:
    我们从 [ 0 , 1 ] [0,1] [0,1]的均匀分布中抽取一个值 r i r_i ri,然后令
    x i = α ( r i 1 − r i ) 1 / λ x_i = \alpha \left(\frac{r_i}{1-r_i} \right)^{1 / \lambda} xi=α(1riri)1/λ
    这种通过变换均匀分布的方式得到的 x i x_i xi等价于从 q ( x ) q(x) q(x)中抽取 x i x_i xi
  4. [ 0 , 1 ] [0,1] [0,1]均匀分布中抽取一个数 u i u_i ui,如果 u i < p ( x i ) M q ( x i ) , u_i < \frac{p(x_i)}{Mq(x_i)}, ui<Mq(xi)p(xi), 那么接受 x i x_i xi这个采样值,否则拒绝它。
  5. 重复第3步和第4步,知道采样的数量达到预定要求。

上文已经提及, λ \lambda λ取值不同会影响 M M M值,从而影响到接受率。为了验证这一点,笔者写了一段代码来做一个实验:
假设我们让 α = 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. 是否接受采样点
      break
    nfail <<- 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%,理论上的接受率是 1 M = 85.49 % \frac{1}{M}=85.49\% M1=85.49%,很接近。采样点的具体分布如下:
在这里插入图片描述
α = 5 , λ = 1 \alpha=5, \lambda=1 α=5,λ=1时,代码模拟出的接受率是 28.46 % 28.46\% 28.46%,理论上的接受率是 1 M = 28.5 % \frac{1}{M}=28.5\% M1=28.5%,也很接近。采样点的具体分布如下:
在这里插入图片描述
我们可以看出,当 λ \lambda λ分别取 3 3 3 1 1 1时,接受率从 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个数字中随机地、等概率地抽取一个数。

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

  1. 首先利用 [ r a n d 7 ( ) − 1 ] × 7 + r a n d 7 ( ) [rand7() - 1] \times 7 + rand7() [rand7()1]×7+rand7()这个表达式实现从1-49这49个数中随机而等概率地抽取一个数 x x x
  2. 如果 x ≤ 40 x \le 40 x40,那么接受 x x x,并将 x % 10 + 1 x \% 10 + 1 x%10+1的值作为一个采样点;否则拒绝 x x x;那么这样得到的一系列采样点就是符合预期要求的。

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

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

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 , … , N 1,2,\ldots,N 1,2,,N N N N个样本中随机抽取 m m m个样本, N N N非常大或者未知,且 m < < N m << N m<<N
蓄水池抽样算法的步骤是:

  1. 首先将 1 , 2 , … , m 1,2,\ldots,m 1,2,,m m m m个样本选进来;
  2. i = m + 1 i=m+1 i=m+1开始,从 [ 1 , i ] [1, i] [1,i]中随机等概率抽取一个数 r r r,如果 1 ≤ r ≤ m 1 \le r \le m 1rm,那么用第 i i i个样本替代第 r r r个样本。否则不做处理
  3. i = i + 1 i=i+1 i=i+1
  4. 重复第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)中抽取的任意一个值记为 X X X;并且让事件 A A A代表 X X X这个值被接受了。如果当 X = x i X=x_i X=xi时被接受了,那么就是 X ^ = X = x i \hat{X}=X=x_i X^=X=xi

有了这些记号,我们可以很容易地知道:
Pr ⁡ ( A ∣ X = x i ) = Pr ⁡ [ u i < p ( x i ) M q ( x i ) ] = p ( x i ) M q ( x i ) \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(AX=xi)=Pr[ui<Mq(xi)p(xi)]=Mq(xi)p(xi)

利用全概率公式,我们知道:
Pr ⁡ ( A ) = ∫ Pr ⁡ ( A ∣ X = x i )   q ( x i ) d x i = ∫ p ( x i ) M q ( x i ) q ( x i ) d x i = 1 M ∫ p ( x i ) d x i = 1 M \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(AX=xi)q(xi)dxi=Mq(xi)p(xi)q(xi)dxi=M1p(xi)dxi=M1

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

利用贝叶斯公式,我们有:
Pr ⁡ ( X ^ ≤ x i ) = Pr ⁡ ( X ≤ x i ∣ A ) = Pr ⁡ ( X ≤ x i , A ) Pr ⁡ ( A ) = ∫ x i q ( t ) p ( t ) M q ( t ) d t Pr ⁡ ( A ) = 1 M ∫ x i p ( t ) d t 1 M = ∫ x i p ( t ) d t \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(XxiA)=Pr(A)Pr(Xxi,A)=Pr(A)xiq(t)Mq(t)p(t)dt=M1M1xip(t)dt=xip(t)dt

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

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

从上面的证明过程中可以看出,接受率就是 Pr ⁡ ( A ) \Pr(A) Pr(A),等于 1 M \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) M M M值并不容易,尤其是在高维的情况下, M M M值往往偏大,从而导致拒绝率很高。

参考

《生物序列分析》第11章

(公众号:生信了)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值