如何加速蒙特卡罗模拟？

Question

每当我在S-Plus中运行大规模蒙特卡罗模拟时，我总是要等待很长时间，期间我会留起胡须。

在R中运行蒙特卡罗模拟的最佳技巧是什么？有没有好的分布式处理进程的示例？

- griffin

Answer 1

• 如果你只是使用并行独立复制，那么使用多个核心/机器应该很简单，但要注意常见的随机数生成器缺陷（例如，如果使用当前时间作为种子，在每个进程中生成一个随机数生成器可能会产生相关的随机数，这会导致无效结果 - 参见例如这篇论文）。

• 您可能需要使用方差缩减来减少所需的复制数量，即缩小所需样本的大小。更高级的方差缩减技术可以在许多教科书中找到，例如在这本书中。

- Roland Ewald

Answer 2

预分配您的向量！

> nsims <- 10000
> n <- 100
> 
> system.time({
     res <- NULL
     for (i in 1:nsims) {
         res <- c(res,mean(rnorm(n)))
     }
 })
   user  system elapsed 
  0.761   0.015   0.783 
> 
> system.time({
     res <- rep(NA, nsims)
     for (i in 1:nsims) {
         res[i] <- mean(rnorm(n))
     }
 })
   user  system elapsed 
  0.485   0.001   0.488 
> 

- Eduardo Leoni

---

我昨天在一个模型中使用了这个，它将我的运行时间缩短了超过15％。绝对值得写一行代码。 - JD Long

Answer 3

我知道这个帖子已经很老了，但如果有人偶然发现并正在寻找一种更快的方法，我认为以下方法是可行的：

library(data.table)
library(microbenchmark)

nsims <- 10000
n <- 100

# Answer from @Eduardo_Leoni:
preallocate<-function(nsims, n) {
  res <- rep(NA, nsims)
  for (i in 1:nsims) {
    res[i] <- mean(rnorm(n))
  }
  return(res)
}

# Answer using data.table:
datatable<-function(nsims,n) {
  dt <- data.table(i=1:nsims)[,list(res=mean(rnorm(1:n))),by=i]
  return(dt)
}

# Timing benchmark:
microbenchmark(preallocate(nsims,n), datatable(nsims,n), times=100)
#Unit: milliseconds
#                  expr      min       lq   median       uq      max neval
# preallocate(nsims, n) 428.4022 432.3249 434.2910 436.4806 489.2061   100
#   datatable(nsims, n) 238.9006 242.3517 244.1229 246.5998 303.6133   100

- dnlbrky

---

你应该在同一台电脑上运行你的解决方案和其他解决方案（因为有些电脑比其他电脑更快）。我发现你的解决方案比使用“microbenchmark”的循环要慢一点；而使用“replicate(nsims, mean(rnorm(1:n)))”则稍微快一些。在这两种情况下，差别大约是3%。 - Frank

感谢建议，@Frank，已经看到编辑后的答案。在我的系统上，data.table比@Eduardo_Leoni的循环解决方案提供了约1.8倍的加速。如果我漏掉了什么，请告诉我，因为我绝对不是基准测试的专家... - dnlbrky

非常有帮助。使用这个技巧，我将计算时间从250秒减少到60秒。 - hubert_farnsworth

Answer 4

拉丁超立方抽样是一种易于应用且对结果产生重大影响的方法。基本上，您可以从均匀分布中获取拉丁超立方样本（例如使用lhs包中的randomLHS()函数），然后使用qnorm(uniformsample)将其转换为所需的分布。

- eli