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

Question

我试图研究一个矩为 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, 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 分布）的随机变量生成。

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

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

Speed of random sampling from distribution in R

random

performance

r

distribution