【问题标题】:Grid-based Monte Carlo simulation using mc2d使用 mc2d 进行基于网格的蒙特卡罗模拟
【发布时间】:2018-01-17 17:37:55
【问题描述】:

我有一个简单的公式来计算四个不同网格之和的“CL”:

CL = A + B + C + D

我正在尝试使用 mc2d 包根据每个输入网格中的不确定性来表达 CL 结果中的不确定性。

每个输入网格中的不确定性基于均匀分布表示,不确定性范围为:

对于光栅 A +/- 6

对于光栅 B,+/- 10

对于光栅 C,+/- 12

对于光栅 D +/- 5

我想运行蒙特卡罗模拟,根据从这些不确定性范围内随机选择的每个网格单元的值对输入网格进行多次求和(例如 n = 1000)

下面的代码失败了:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
    cannot use this function

我不确定如何解决此错误。任何建议将不胜感激。

library(raster) 
library(rgdal) 
library(mc2d) #load package

# Input grid files
A <- raster(matrix(c(20, 35, 40, 60), nrow=2))
B <- raster(matrix(c(6, 10, 13, 14), nrow=2))
C <- raster(matrix(c(6, 8, 12, 14), nrow=2))
D <- raster(matrix(c(35, 40, 50, 60), nrow=2))

# combine the RasterLayer objects into a RasterStack 
s <- stack(A, B, C, D)

# Uncertainty distance for each raster used to establish the 
# min/max for the uniform distributions in the function just below
uA <- 6
uB <- 10
uC <- 12
uD <- 5

# Function for generating CL
ndunc(1000)
fun <- function(x) {
  # Convert each raster to mcnode object
  dA <- mcdata(x[1], type="0")
  dB <- mcdata(x[2], type="0")
  dC <- mcdata(x[3], type="0")
  dD <- mcdata(x[4], type="0")

  # Define uniform distirbutions for each raster
  stA <- mcstoc(runif, type="U", min=dA-uA, max=dA+uA, rtrunc=TRUE, linf=0)
  stB <- mcstoc(runif, type="U", min=dA-uB, max=dA+uB, rtrunc=TRUE, linf=0)
  stC <- mcstoc(runif, type="U", min=dA-uC, max=dA+uC, rtrunc=TRUE, linf=0)
  stD <- mcstoc(runif, type="U", min=dA-uD, max=dA+uD, rtrunc=TRUE, linf=0)

  # Apply Monte Carlo to this equation
  CL <- stA + stB + stC + stD
  # Extract the min, mean, and max CL from the 1000 iterations
  c(min(CL), mean(CL), max(CL))
}

# Need to create rasters for min(CL), mean(CL), and max(CL)
x <- calc(s, fun) 
names(x) = c('min', 'mean', 'max') 
plot(x)

【问题讨论】:

  • 原始问题已被修改为包含在其他地方提供的答案。代码现在可以正常工作了。

标签: r montecarlo


【解决方案1】:

函数中的错误是由 rtrunc=TRUE 中的拼写错误(缺少等号)引起的

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2019-05-09
    • 1970-01-01
    • 1970-01-01
    • 2023-01-24
    • 1970-01-01
    • 2013-09-24
    相关资源
    最近更新 更多