生成无循环的随机数

Question

我试图尽可能多地减少一个函数的执行时间,该函数对伯努利序列序列的输出求和.

这是我工作但速度慢的方法:

set.seed(28100)
sim <- data.frame(result = rep(NA, 10))
for (i in 1:nrow(sim)) {
  sim$result[i] <- sum(rbinom(1200, size = 1, prob = 0.2))
}
sim
# result
# 1     268
# 2     230
# 3     223
# 4     242
# 5     224
# 6     218
# 7     237
# 8     254
# 9     227
# 10    247

如果没有for循环,我怎么能得到相同的结果？

我试过这个......

set.seed(28100)
sim <- data.frame(result = rep(sum(rbinom(1200, size = 1, prob = 0.2)), 10))
sim
# result
# 1     269
# 2     269
# 3     269
# 4     269
# 5     269
# 6     269
# 7     269
# 8     269
# 9     269
# 10    269

但很明显,这个论点rep()只执行了一次.

Answer 1

二项分布定义为伯努利试验的总和.

# this line from your question
sum(rbinom(1200, size = 1, prob = 0.2))
# is equivalent to this
rbinom(1, size = 1200, prob = 0.2)

# and replicating it
replicate(expr = sum(rbinom(1200, size = 1, prob = 0.2)), n = 10)
# is equivalent to setting n higher:

        ### This is the only line of code you need! ####
rbinom(10, size = 1200, prob = 0.2)

在我的(相当慢的)笔记本电脑上进行100,000次模拟需要大约0.01秒,对于1M模拟需要0.12秒.

修改@ eipi的漂亮基准测试,这比其他方法快700-900倍(现在有bug修复!)

          expr     min      lq       mean  median      uq     max neval cld
         binom   1.324   1.377   1.607959   1.413   1.931   2.306    10 a  
     replicate 716.300 737.200 756.288641 749.900 765.300 812.400    10  b 
        sapply 706.300 743.300 778.863587 763.800 853.500 860.300    10  b 
 matrixColSums 838.800 870.000 893.813083 894.800 907.500 978.200    10   c

基准代码:

nn = 10000
n_bern = 1200
library(microbenchmark)
print(
    microbenchmark::microbenchmark(
        replicate =
            replicate(nn, sum(rbinom(
                n_bern, size = 1, prob = 0.2
            )))
        ,
        matrixColSums =
            colSums(matrix(
                rbinom(n_bern * nn, size = 1, prob = 0.2), ncol = nn
            )),
        sapply = sapply(
            1:nn,
            FUN = function(x) {
                sum(rbinom(n_bern, size = 1, prob = 0.2))
            }
        ),
        binom = rbinom(nn, size = n_bern, prob = 0.2),
        times = 10
    ),
    order = "median",
    signif = 4
)