生成没有循环的随机数答案

【问题标题】：Generate random numbers without loops生成没有循环的随机数
【发布时间】：2016-08-07 15:51:39
【问题描述】：

我正在尝试尽可能多地减少函数的执行时间，该函数对一系列伯努利试验的输出求和。

这是我的工作但缓慢的方法：

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() 的参数只执行一次。

【问题讨论】：

我会把你的答案指向我的解释，但单线解决方案是rbinom(10, size = 1200, prob = 0.2)。

标签： r performance loops random

【解决方案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 倍（现在修复了错误！）

          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
)

【讨论】：

如果您在答案中查看我的replicate 代码，您会发现我没有正确的参数值（1 和 12 而不是 1200 和 1）。我正在朝着类似于您的答案的方向前进，但我想我一定是在中间进行计时，而不是事先进行计时。无论如何，replicate 并不比其他两种方法快，而你的方法显然是要走的路。我只是想让您知道，以便您可以更正 replicate 方法的代码和时间（我已经更正了我的答案）。
这很有趣，当我参数化 nn 进行基准测试时，我也开始提取 n_bernoulli = 1200，但是当我拿到你的代码时，你只有 12 个 - 我以为你在做一些花哨的东西在其他地方解释它 - 我没有花时间去考虑它。

【解决方案2】：

set.seed(28100)
nsim=10
sim = data.frame(result=replicate(nsim, sum(rbinom(1200, size=1, prob=0.2))))

sim

以下是具有 10,000 次模拟的各种方法的一些时间安排：

microbenchmark::microbenchmark(
  replicate = {nsim=10000
  data.frame(result=replicate(nsim, sum(rbinom(1200, size=1, prob=0.2))))},
  matrixColSums = {
    sims <- 10000
    n <- 1200
    r <- rbinom(n*sims, size = 1, prob = 0.2)
    r <- matrix(r, ncol=sims)
    data.frame(result=colSums(r)) },
  sapply = data.frame(result=sapply(1:10000, FUN = function(x) {sum(rbinom(1200, size = 1, prob = 0.2))})),
  times=10
)

Unit: milliseconds
         expr      min       lq     mean   median       uq      max neval cld
    replicate 584.2389 597.5571 615.7545 614.0977 630.7354 648.8328    10  a 
matrixColSums 655.0608 664.2053 684.0069 682.1868 702.1426 713.0240    10   b
       sapply 589.9830 610.5784 626.8738 629.2161 642.2589 660.6092    10  a

【讨论】：

【解决方案3】：

矢量化是关键。

主要的节省时间（至少对于大n）是使用sample。

例如对于

n <- 1e7
sample(0:1, n, replace=TRUE)

大约需要 0.2 秒，而

for(i in 1:n) sample(0:1, 1)

大约需要 24 秒。向量化操作通常可以替代循环，但知道何时何地取决于您是否熟悉满足您需求的可用函数。

【讨论】：

【解决方案4】：

这样怎么样：

set.seed(28100)
sims <- 10
n <- 1200
r <- rbinom(n*sims, size = 1, prob = 0.2)
r <- matrix(r, ncol=sims)
colSums(r)

对我而言，100,000 次模拟的速度大约是其两倍（6 秒对 13 秒），但 R. Schifini 和 eipi10 的解决方案要快一些（约 5.5 秒）

【讨论】：

【解决方案5】：

执行以下操作：

sim = rep(NA, 10)
sapply(sim,FUN = function(x) {sum(rbinom(1200, size = 1, prob = 0.2))})

结果：

[1] 216 231 234 249 249 236 255 251 231 244

然后转换为数据框

【讨论】：