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

Question

我试图研究一个概率分布，它的矩是加泰罗尼亚数，并想出了

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, 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,4),
  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 版本更快？
• 在这种情况下，是否有关于如何最好地制作随机样本的一般建议？

Answer 1

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

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