从 R 中的分布随机抽样的速度

问题描述

我试图调查其矩为 Catalan numbers 的概率分布，并得出了

qcatmo <- function(p,k=4){ (qbeta(p/2+1/2,3/2,3/2)*2 - 1)^2 * k } 
colMeans(outer(qcatmo(ppoints(10^6)),0:10,"^"))
#      1     1     2     5    14    42   132   429  1430  4862 16796

效果很好。但是后来我尝试从这个分布中生成随机值，并找到了三种可能的方法（A 使用我已经知道的分位数函数应用于 runif，B 使用内置的 rbeta 函数稍微更直接,和 C 使用一种带有 runif) 的拒绝采样形式，在大样本上使用时具有明显不同的速度：

rcatmoA <- function(n,k=4){ qcatmo(runif(n),k) }
rcatmoB <- function(n,k=4){ (rbeta(n,3/2)*2 - 1)^2 * k }
rcatmoC <- function(n,k=4){
             n0 <- ceiling(n*4/pi + 7*sqrt(n) + 35)
             x0 <- runif(n0)^2 
             y0 <- runif(n0)^2 
             x0[x0 + y0 < 1][1:n] * k
             }

进行基准测试时给出了

library(microbenchmark)
n <- 10^4
microbenchmark(
  rcatmoA(n,4),rcatmoB(n,rcatmoC(n,4)
  ) 
#Unit: microseconds
#          expr     min       lq      mean   median       uq     max neval cld
# rcatmoA(n,4) 22817.2 23014.95 23259.688 23186.95 23322.80 25128.9   100   c
# rcatmoB(n,4)  1526.5  1534.40  1615.255  1541.30  1607.15  4952.1   100  b 
# rcatmoC(n,4)   781.5   788.70   884.339   795.00   813.80  7266.2   100 a  

我的问题是：

• 为什么 B 版比 A 版快这么多？
• 如果 B 版本更快是因为它避免将函数应用于 runif 数据，那么为什么 C 版本更快？
• 对于在这种情况下如何最好地制作随机样本有什么一般性建议吗？

解决方法

正如您所发现的，有不同的方法可以执行相同分布（在本例中为 beta 分布）的随机变量生成。

1. 我认为 A 版本是最慢的，因为 beta 分布的分位数没有封闭形式，因此 qbeta 必须求助于数值反演（变得更加困难，因为 beta 分布的 CDF，正则化 beta 函数，是通常仅称为积分）。事实上，正如 qbeta 的源代码所示，分位数函数远非微不足道。 Look at the R source code。
2. B 版本所依赖的 rbeta 的源代码使用了 Cheng (1978) 的 beta 分发算法。在这种特殊情况下，使用算法 BB，它使用拒绝采样器（并且每次迭代至少需要两个对数）。相比之下，C 版本采用了更简单的拒绝条件。 Cheng 的算法是多年来为 beta 分布提出的众多算法之一，L. Devroye 在 Non-Uniform Random Variate Generation 中提到了 1986 年及更早的算法。