时间戳

首页 > 解决方案 > 加快蒙特卡洛模拟？

r - 加快蒙特卡洛模拟？

问题描述

我目前正在进行研究，通过使用蒙特卡洛模拟来预测葛根（一种侵入性藤蔓）在五年内将在俄克拉荷马州传播的地方。我创建了一个带有存在点的栅格并将其加载到 R 中。

对于每个蒙特卡罗模拟（每年），我运行 6000 次迭代以提供准确的结果。但是，由于“for”循环，这需要很长时间。第一年通常在 3 天内完成运行，但第二年已经运行了 3 周以上，仍未完成。

有没有办法加快这个过程？

每一年都建立在前一年的基础上。我在前两年提供了以下代码：

OK.rast16<-raster("OK_rast20161.tif")
p.a16<-as.matrix(OK.rast16)
table(p.a16)

# Set the random number seed so results can be reproduced if needed
set.seed(10)

drow.pa16 <- 133.7197873 # distance of grid cell (meters) in n-s direction
trow.pa16 <- length(p.a16[,1]) # total number of rows

dcol.pa16 <- 133.7197873 # distance of grid cell (meters) in e-w direction
tcol.pa16 <- length(p.a16[1,]) # total number of rows

#####Year 1 of infection in 2016#####

kudzu_sim1_16 <- matrix(0,trow.pa16,tcol.pa16)
for(m in 1:6000)
{
  OK.kudzu_16 <- p.a16 # initialize matrix of annual dispersal
  for(i in 1:trow.pa16)
{
  for(j in 1:tcol.pa16)
{
  if(!is.na(p.a16[i,j]) & p.a16[i,j] == 1)
  {
    for(k in 1:trow.pa16)
    {
      for(l in 1:tcol.pa16)
      {
        if(!is.na(OK.kudzu_16[k,l]) & OK.kudzu_16[k,l] == 0)
        {
          distcalc <- sqrt((abs(i-k)*drow.pa16)^2+(abs(j-l)*dcol.pa16)^2)
          prob <- exp(0.0369599-0.00474401*distcalc)
          if(prob>runif(1)) {OK.kudzu_16[k,l] <- 1}
        }
      }
    }
   }
  }
 }
 kudzu_sim1_16 <- OK.kudzu_16+kudzu_sim1_16
}

#####Year 2 of infection in 2016####

kudzu_sim2_16 <- matrix(0,trow.pa16,tcol.pa16)
for(m in 1:6000)
{
  OK.kudzu1_16 <- OK.kudzu_16 # initialize matrix of annual dispersal
   for(i in 1:trow.pa16)
  {
    for(j in 1:tcol.pa16)
   {
     if(!is.na(OK.kudzu_16[i,j]) & OK.kudzu_16[i,j] == 1)
  {
    for(k in 1:trow.pa16)
    {
      for(l in 1:tcol.pa16)
      {
        if(!is.na(OK.kudzu1_16[k,l]) & OK.kudzu1_16[k,l] == 0)
        {
          distcalc <- sqrt((abs(i-k)*drow.pa16)^2+(abs(j-l)*dcol.pa16)^2)
          prob <- exp(0.0369599-0.00474401*distcalc)
          if(prob>runif(1)) {OK.kudzu1_16[k,l] <- 1}
        }
      }
    }
  }
 }
}
kudzu_sim2_16 <- OK.kudzu1_16+kudzu_sim2_16
}

这是要加载以启动代码的栅格：

葛根好

标签： rfor-loopiterationsimulationmontecarlo

解决方案

原本是一条评论，但很快就超过了长度限制：

1) 133mx133m 是一个非常小的网格尺寸，用于对整个州这么大的东西进行空间模拟。它可能有助于找到一种方法使分辨率成为模拟的参数（而不是硬连线数字），简化代码以使其在更大的分辨率下运行良好，然后提高分辨率。该raster函数有一个名为的可选参数res，可用于控制分辨率。

2）虽然矢量化肯定会有所帮助，但不太可能将运行数周且没有输出的算法转换为在合理时间内运行的算法。也许您需要从根本上重新考虑您的算法。您似乎正在将每个网格单元与每个其他网格单元进行比较。这在我看来在生物学上并不合理。葛根在当地传播。为什么明年俄克拉荷马州狭长地带的一个给定的 133m x 133m 小区会发生什么取决于 Texoma 湖上另一个小区的当前状态？如果您的模拟具有任何生物现实主义，exp(0.0369599-0.00474401*distcalc)那么对于两个这样的细胞应该可以忽略不计，但您的代码不会忽略它。在某种意义上本地化的算法可能会更好。

3）您的矩阵中有很多条目对应于俄克拉荷马州以外的点。除非您的模型旨在查看 kudzu 如何在德克萨斯州的大部分地区扩散，否则这些可能与您的程序无关。如果是这样，一个根本不同的数据结构（例如位置列表）可能更可取，它只包含俄克拉荷马州的点条目。或者可能不是。只是想一想。

4）要获得更详细的帮助，如果您解释您的算法实际上是什么（而不仅仅是它打算做什么），将会有所帮助。快速阅读您的代码并不是很明显。

推荐阅读