R代码: Is there a way to make this Monte Carlo simulation quicker?

Question

想象一下我递给你一个上面印有“-1”的乒乓球。然后我告诉你从标有“第一袋”的袋子中取出另一个乒乓球。这个袋子里有 30,000 个球，有的标有“-1”，有的标有“0”，有的标有“+1”。无论您抽到哪个球，您都将其数字添加到您当前的“分数”-1 中。例如，如果您抽到 -1，您的新分数为 -2。

只要您的新分数低于零，您就可以再次从第一个袋子中抽取并再次添加到您的分数中。但是，如果您的分数达到零或更高，您就会从第二个袋子中抽奖，该第二个袋子的组成不同，包括 -1、0 和 +1。

我希望您从相应的袋子中总共抽取 1000 个乒乓球(即，取决于您当前的分数是否低于零)，然后在“的末尾写下您的总(累积)分数”放”。然后我希望你重复这个实验一百万次，并告诉我你最终得分高于零的组的百分比是多少。

有没有更快/更有效的方法来编码？由于绘制不是独立的，所以很难对循环进行矢量化，尽管也许我可以使用 ifelse 和 filter 的组合？我怀疑复制才是最昂贵的部分。

ptm <- proc.time()

###First bag
n=30000
s=155
f=255
z=n-s-f
first_bag=c(rep(0,z), rep(1,s), rep(-1,f))

###Second bag
n2=30000
s2=275
f2=285
z2=n2-s2-f2
second_bag=c(rep(0,z2), rep(1,s2), rep(-1,f2))

###Simulate draws
sim_draw=function(draws){

  score=-1

  for (i in 1:draws) {
    if (score < 0) {
      score=score + sample(first_bag, 1, replace=TRUE)} else {
      score=score + sample(second_bag, 1, replace=TRUE)}
  }
  score
}

###Repeat sims and find area above zero
samp_distribution=replicate(1000000, sim_draw(1000))
mean(samp_distribution>0)


print(proc.time() - ptm)

Answer 1

一些想法:

1. R 确实不擅长这种迭代过程。编译语言的性能会好得多。我在这里假设您想坚持使用基本的 R。

2. 学习使用探查器来查看您的实现在哪些方面浪费了时间。请参阅 ?summaryRprof 底部的示例使用方法只需替换example(glm)即可与您的代码:samp_distribution <- replicate(1000, sim_draw(1000)) 。您会发现，调用 sample 浪费了很多时间。一遍又一遍。因此，对代码的第一个改进可以是调用 sample只有几次:

   sim_draw_1 <- function(draws){
   
      s1 <- sample(bag1, draws, replace = TRUE)
      s2 <- sample(bag2, draws, replace = TRUE)
      score <- -1
   
      for (i in 1:draws)
         score <- score + if (score < 0) s1[i] else s2[i]
   
      score
   }
   

看到这快了近十倍(我发现微基准包是测量/比较计算时间的更可靠方法)

    library(microbenchmark)
    microbenchmark(sim_draw(1000), sim_draw_1(1000),
                   times = 1000)
    # Unit: microseconds
    #              expr      min       lq    median       uq       max neval
    #    sim_draw(1000) 5518.758 5845.465 6036.1375 6340.662 53023.483  1000
    #  sim_draw_1(1000)  690.796  730.292  743.8435  785.071  8248.163  1000

对于像您这样的迭代代码，总是值得尝试编译器:

library(compiler)
sim_draw_2 <- cmpfun(sim_draw_1)
library(microbenchmark)
microbenchmark(sim_draw_1(1000), sim_draw_2(1000), times = 1000)
# Unit: microseconds
#              expr     min       lq   median      uq      max neval
#  sim_draw_1(1000) 684.687 717.6640 748.3305 812.971 9412.936  1000
#  sim_draw_2(1000) 242.895 259.8125 268.3925 294.343 1710.290  1000

又提高了 3 倍，还不错。

最后，由于函数内最大的瓶颈仍然是 for 循环，因此您可以尝试重写它，这样您就可以使用矢量化(快速)函数来处理，而不是一次处理一个结果尽可能多的结果(迫使你转换帽子的结果的确切数量。)

sim_draw_3 <- function(draws, bag1 = first_bag,
                       bag2 = second_bag){

   s1 <- sample(bag1, draws, replace = TRUE)
   s2 <- sample(bag2, draws, replace = TRUE)

   score <- -1L
   idx   <- 1L
   while (idx <= draws) {
      bag         <- if (score < 0) s1 else s2
      switch.at   <- if (score < 0) 0L else -1L
      next.draws  <- bag[idx:draws]
      next.scores <- score + cumsum(next.draws)
      stop.idx    <- which(next.scores == switch.at)[1]
      if (is.na(stop.idx)) stop.idx <- length(next.draws)
      score <- next.scores[stop.idx]
      idx   <- idx + stop.idx
   } 
   score
}
sim_draw_4 <- cmpfun(sim_draw_3)

microbenchmark(sim_draw_2(1000), sim_draw_3(1000), sim_draw_4(1000), times = 1000)
# Unit: microseconds
#              expr     min      lq   median       uq      max neval
#  sim_draw_2(1000) 236.916 252.540 269.1355 293.7775 7819.841  1000
#  sim_draw_3(1000)  80.527  95.185 128.9840 162.7790  625.986  1000
#  sim_draw_4(1000)  79.486  92.378 123.5535 162.5085  518.594  1000

又一个 2 倍的改进。在这里您可以看到，编译器只给我们带来了微小的改进，因为迭代次数急剧下降，并且 R 代码中的其他所有内容都使用非常高效的(向量化)函数。

因此，每个函数调用的时间从 5845 微秒减少到 124 微秒，这是一个相当不错的改进。如果这仍然太慢，那么您可能必须切换到 c++(例如通过 Rcpp)。至少我希望这有助于向您展示一些有用的技巧。

最后但并非最不重要的一点是，由于您的函数调用都是独立的，因此您可以考虑并行运行它们。我会告诉你http://cran.r-project.org/web/views/HighPerformanceComputing.html并鼓励您四处搜索。