Код R: Есть ли способ ускорить симуляцию Монте-Карло?

Несколько идей:

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
   }
   

Обратите внимание, что это почти в десять раз быстрее (я считаю пакет microbenchmark более надежным методом измерения/сравнения времени вычислений).

    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, вы можете попытаться переписать его так, чтобы вместо обработки одного результата за раз вы использовали векторизованные (быстрые) функции для обработки как можно большего количества результатов (цикл 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
   

Еще одно двукратное улучшение. Здесь вы видите, что компилятор дает нам лишь незначительное улучшение, потому что количество итераций резко сократилось, а все остальное в нашем коде R использует очень эффективные (векторизованные) функции.

Таким образом, мы получили от 5845 микросекунд до 124 на вызов функции, что является довольно хорошим улучшением. Если это все еще слишком медленно, вам, вероятно, придется переключиться на c++ (например, через Rcpp). По крайней мере, я надеюсь, что это помогло показать вам некоторые полезные приемы.

И последнее, но не менее важное: я хотел бы упомянуть, что, поскольку все ваши вызовы функций независимы, вы можете рассмотреть их параллельное выполнение. Я укажу вам на http://cran.r-project.org/web/views/HighPerformanceComputing.html и предложить вам поискать.