R代码:有没有办法让这个蒙特卡罗模拟更快？

Question

想象一下,我递给你一个乒乓球,上面印有"-1".然后我告诉你从一个标有"First Bag"的包中画另一个乒乓球.这个包有30,000个球,一些标有"-1",一些标有"0",有些标有"+1".无论你绘制哪个球,你都要将它的数字加到你当前的"得分"-1.例如,如果你画-1,你的新分数是-2.

只要您的新分数低于零,您就可以从第一个包中再次绘制并再次添加到您的分数中.但是,如果你的分数达到零或更高,你可以从第二个包中抽取,它具有-1s 0s和+ 1s的不同组成.

我希望你从相应的袋子里抽出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

3. 对于像你这样的非常迭代的代码,总是值得尝试编译器:

   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倍的改善,还不错.

4. 最后,因为函数内部仍然是最大的瓶颈是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并鼓励您搜索.