使用snowfall在R中进行并行蒙特卡罗模拟

Question

我尝试比较多达数千个估计的beta分布。每个beta分布由两个形状参数alpha和beta来描述。现在，我从每个分布中抽取10万个样本。最终结果是要得到每个样本抽取中概率最高的分布顺序。我的第一个方法是使用lapply生成一个N * NDRAWS数字值的矩阵，当N超过10000时会消耗太多内存（10000 * 100000 * 8字节）。因此，我决定采用顺序方法对每个单独的抽样进行排序，然后累加所有抽样的顺序，得到最终的顺序，如下面的示例所示：

set.seed(12345)
N=100
NDRAWS=100000
df <- data.frame(alpha=sample(1:20, N, replace=T), beta=sample(1:200, N, replace=T))  

vec    <- vector(mode = "integer", length = N )

for(i in 1:NDRAWS){
  # order probabilities after a single draw for every theta
  pos <- order(rbeta(N, shape1=df$alpha, shape2=df$beta) )

  # sum up winning positions for every theta
  vec[pos] <- vec[pos] + 1:N
}

# order thetas
ord <- order(-vec)

df[ord,]

这只消耗了N * 4字节的内存，因为没有巨大的矩阵，而是一个长度为N的单个向量。我的问题是，如何利用我的4个CPU核心，而不仅仅使用一个核心，通过使用snowfall（或任何其他多核包）来加速此操作？

# parallelize using snowfall pckg
library(snowfall)
sfInit( parallel=TRUE, cpus=4, type="SOCK")
sfLapply( 1:NDRAWS, function(x) ?????? )
sfStop()

非常感谢您的帮助！

- Christian Borck

---

问题是如何将for循环转换成sfLapply，以便最终排序的数据框与我的示例中的相同？ - Christian Borck

这让我感到困惑：vec[pos] <- vec[pos] + 1:N。你试图将一个向量放入一个值中？ - PascalVKooten

我知道这可能有点令人困惑，也许有更直接的方法。我的做法是为每个1:N绘制添加单个有序绘制的位置插槽。 - Christian Borck

啊，好吧。我现在的看法是，如果你要将结果相加到每次迭代中的一个单独项（该项不在内存中共享），那么它就不是顺序的了。 - PascalVKooten

感谢您的努力和指出foreach和doParallel。还有您在StackOverflow上给新手的提示;)我会在未来的讨论中记住这些。 - Christian Borck

Answer 1

这可以使用与随机森林或自助法并行化相同的方式进行并行化。只需在每个工作线程上执行顺序代码，但每个线程使用更少的迭代次数。这比将for循环的每个迭代拆分为单独的并行任务要高效得多。

下面是使用foreach包和doParallel后端的完整示例：

set.seed(12345)
N=100
NDRAWS=100000
df <- data.frame(alpha=sample(1:20, N, replace=T),
                 beta=sample(1:200, N, replace=T))
library(doParallel)
nworkers <- detectCores()
cl <- makePSOCKcluster(nworkers)
clusterSetRNGStream(cl, c(1,2,3,4,5,6,7))
registerDoParallel(cl)

vec <- foreach(ndraws=rep(ceiling(NDRAWS/nworkers), nworkers),
               .combine='+') %dopar% {
  v <- integer(N)
  for(i in 1:ndraws) {
    pos <- order(rbeta(N, shape1=df$alpha, shape2=df$beta) )
    v[pos] <- v[pos] + 1:N
  }
  v
}
ord <- order(-vec)
df[ord,]

请注意，这与顺序版本给出不同的结果，因为工作进程生成了不同的随机数。我使用了并行程序包提供的并行随机数支持，因为这是一个良好的实践。

- Steve Weston

---

谢谢你，史蒂夫。这正是我一直在寻找的。 正如你所提到的，由于随机数发生器的影响，结果略有不同，但这是一个很好的权衡，考虑到速度的提高。 你知道这是否在Linux系统上也适用吗？核心会自动检测吗？ - Christian Borck

是的，我在Linux上编写并测试了它，但它也应该可以在Windows和Mac OS X上运行。我使用了detectCores函数来确定核心数。 - Steve Weston

Answer 2

好的，功能已经实现了。但我不确定每次迭代时你会返回什么。

也许可以试试这个？

myFunc <- function(xx, N) {
  pos <- order(rbeta(N, shape1=df$alpha, shape2=df$beta) )
  vec[pos] + 1:N
}

使用 doParallel 可以让你添加结果：

require(doParallel)
registerDoParallel(cores=4)
foreach(i=1:NDRAWS, .combine='+') %dopar% myFunc(i, N)

- PascalVKooten

---

嗨Dualinity，感谢您的快速回复。我不小心发布了一个示例函数“sum(rnorm(x))”。实际上，我想将我的for循环代码示例转换为可以在并行处理中处理而不需要构建巨大矩阵的东西。 - Christian Borck

是的，我的错，抱歉。它刚刚在我面前爆炸了。 - PascalVKooten

还有一个并行的 for 循环，这听起来不像是你正在寻找的吗？ - PascalVKooten

@ChristianBorck 通常，如果您想补充这篇帖子，方法是发布自己的答案，在其中引用我的帮助并提供完整的解决方案。或者，如果您仍需要帮助，则可以提出进一步的问题。如果这个答案有所帮助，您可以在拥有足够声望时投票支持它或回答它。欢迎来到StackOverflow。 - PascalVKooten