为什么 Beta(2,2) 的接受概率为 2/3 的拒绝抽样不比 `rbeta(N,2,2)` 慢？

机器算法验证 模拟 贝塔分布 拒绝抽样
2022-03-23 15:59:33

我试图说明当有一种不丢弃样本的替代方法可用时，拒绝抽样是低效的。（这篇文章显然是迁移到 SO 的候选者，但我觉得它可能还有一个统计方面。）

我的例子是模拟Beta(2,2)随拒绝算法而变化，并且 rbeta(N,2,2). 我的代码如下：

# densityatmode <- dbeta(.5,2,2)=1.5
rejectionsampling_beta22 <- function(N){
     px <- runif(N,min=0,max=1)
     py <- runif(N,min=0,max=1.5)
     ikeep <- (py < 6*px*(1-px))
     return(px[ikeep])
}

samples <- 1500000

接受概率area under density/area in box当然area under density等于 1，因此是 1/1.5，或者等效地，当使用拒绝时，我们需要 1.5 个候选样本。

library(microbenchmark)
result <- microbenchmark(rejectionsampling_beta22(1.5*samples),rbeta(samples,2,2))

结果令我惊讶，因为我觉得拒绝 1/3 的平局应该是浪费：

Warning message:
In microbenchmark(rejectionsampling_beta22(1.5 * samples), rbeta(samples,  :
  Could not measure a positive execution time for 33 evaluations.
> result
Unit: nanoseconds
                                    expr       min        lq         mean    median        uq       max neval cld
 rejectionsampling_beta22(1.5 * samples) 180171443 254891504 248482870.58 257005186 258866373 329533988   100  b 
                    rbeta(samples, 2, 2) 253091893 254717640 261602707.69 257391248 259902160 333055032   100   c
                                   neval         0         0        63.97         1         1       604   100 a
2个回答

这是Cheng & Feast Gamma 生成器代码，R rbeta 和 rgamma 函数基于该代码：

function x=gamrnd_cheng(alpha)
% Gamma(alpha,1) generator using Cheng--Feast method
% Algorithm 4.35
c1=alpha-1; c2=(alpha-1/(6*alpha))/c1; c3=2/c1; c4=1+c3;
c5=1/sqrt(alpha);
flag=0;
while flag==0;
    U1=rand; U2=rand;
    if alpha>2.5
        U1=U2+c5*(1-1.86*U1);
    end
       W=c2*U2/U1;
     flag=(U1<1)&&(U1>0)&&(((c3*U1+W+1/W)<c4)||((c3*log(U1)-log(W)+W)<1));
end
x=c1*W;

可以以大约相同的成本将其回收到 Beta 生成器中。它使用两个制服，加上一个拒绝条件，所以对于(α,β)您选择的，即拒绝概率为1/3，接受-拒绝算法可能同样有效。但是，您还应该对较大的非整数值运行比较(α,β)检查 Cheng & Feast Gamma 发生器是否仍然有效。

例如，Joe Whittaker 的 BetaB(α,β)生成器具有以下形式的拒绝条件

U11/α+U21/β>1
这随着频率的增加而发生α和β增加。我记得 Luc Devroye 提到过，因为G(α,1)分布，不可能找到独立于计算时间的界限α...

按照西安的建议，问题的结果确实似乎是所选特定（小）参数的产物。特别是选择更高的α和β, 使得 β 密度变得更集中，模式处的密度相应更高，因此在基本拒绝算法中需要围绕 β 密度构建一个更大的框，使得接受概率降低，结果消失，rbeta并且更高效。

这是一个示例α=β=10.2：

# densityatmode <- dbeta(.5,10.2,10.2)=3.559877
# precompute the beta density:
# gamma(20.4)/gamma(10.2)^2=1231365
rejectionsampling_betaab <- function(N){
  px <- runif(N,min=0,max=1)
  py <- runif(N,min=0,max=3.559877)
  #ikeep <- (py < dbeta(px,10.2,10.2))
  ikeep <- (py < 1231365*px^9.2*(1-px)^9.2)
  return(px[ikeep])
}

AcceptanceProbability <- 1/(3.559877)
samples <- 1e5

library(microbenchmark)
result <- microbenchmark(rejectionsampling_betaab(1/AcceptanceProbability*samples),rbeta(samples,10.2,10.2))

现在的结果显然有利于rbeta：

> result
Unit: milliseconds
                                                        expr       min        lq      mean    median        uq      max neval
 rejectionsampling_betaab(1/AcceptanceProbability * samples) 417.46861 426.16041 469.10049 430.37921 485.60090 635.9086   100
                                  rbeta(samples, 10.2, 10.2)  90.39748  90.94824  94.24759  91.52765  92.51329 264.5119   100
其它你可能感兴趣的问题
上一篇简单的 kNN 示例 下一篇为什么 SGD 和反向传播适用于 ReLU？