r - 与 adaptIntegrate 的集成不一致

Question

我正在估计一个物种在网格景观中扩散的概率，给定一个具有最大扩散距离的扩散核(距离的函数)。我正在尝试计算等式中描述的区域到区域的扩散概率。 8 this (open access) paper .这涉及四重积分，分别针对源单元和目标单元中源点和目标点的每种可能组合评估分散函数的值。

我已经使用 cubature 包中的 adaptIntegrate 实现了这一点，如下所示，对于源单元 A、目标单元 B 和分散为 1 的简化分散内核当点间距离 > 1.25 时，否则为 0。这如下图所示，其中无法到达单元格 B 的红色区域，因为单元格 A 中的任何点都不在 1.25 的距离内。

library(cubature)
f <- function(xmin, xmax, ymin, ymax) {
  adaptIntegrate(function(x) {
    r <- sqrt((x[3] - x[1])^2 + (x[4] - x[2])^2)
    ifelse(r > 1.25, 0, 1)
  }, 
  lowerLimit=c(-0.5, -0.5, xmin, ymin), 
  upperLimit=c(0.5, 0.5, xmax, ymax), 
  maxEval=1e5)
}

f(xmin=1.5, xmax=2.5, ymin=-0.5, ymax=0.5)

# $integral
# [1] 0.01949567
# 
# $error
# [1] 0.001225998
# 
# $functionEvaluations
# [1] 100035
# 
# $returnCode
# [1] 0

当考虑目标单元格 C 时，我得到了不同的积分，它与单元格 A 的距离相同，但位于单元格 A 的上方而不是右侧。

f(xmin=-0.5, xmax=0.5, ymin=1.5, ymax=2.5)

# $integral
# [1] 0.01016105
# 
# $error
# [1] 0.0241325
# 
# $functionEvaluations
# [1] 100035
# 
# $returnCode
# [1] 0

为什么这些积分如此不同(0.01949567 与 0.01016105)？我编码不正确吗？更改公差和最大评估次数似乎没有太大区别。或者，是否有更好的方法来编写此类问题的解决方案？

我意识到有关一般方法的问题可能更适合 stats.stackexchange.com，但我已在此处发布，因为我怀疑我可能忽略了编码本身的某些内容。

---

编辑:对于 A -> B 情况，嵌套的 integrate 返回类似于第一个 adaptIntegrate 解决方案的解决方案。对于 A -> C 情况，它返回 Error in integrate(function(ky) { : the integral is probably divergent。

g <- function(Bx, By, Ax, Ay) {
  r <- sqrt((Ax - Bx)^2 + (Ay - By)^2)
  ifelse(r > 1.25, 0, 1)
}

integrate(function(Ay) {
  sapply(Ay, function(Ay) {
    integrate(function(Ax) {
      sapply(Ax, function(Ax) {
        integrate(function(By) {
          sapply(By, function(By) {
            integrate(function(Bx) g(Bx, By, Ax, Ay), 1.5, 2.5)$value # Bx
          })
        }, -0.5, 0.5)$value # By
      })
    }, -0.5, 0.5)$value # Ax
  })
}, -0.5, 0.5)$value # Ay

# [1] 0.019593

Answer 1

原因似乎是 adaptIntegrate 的工作方式，因为很明显，您唯一需要更改的是集成顺序。不相同的结果可能是因为单独的近似积分(参见第一个响应 here )，但这似乎更像是一个错误。

这里是计算f(xmin=1.5, xmax=2.5, ymin=-0.5, ymax=0.5)时r的值

和f(xmin=-0.5, xmax=0.5, ymin=1.5, ymax=2.5)

所以函数内部一定有什么东西在发生，因为值的范围有很大的不同。

对此的一种替代方法是蒙特卡洛积分，这在这种情况下很好，因为您的点是均匀分布的。

MCI <- function(Ax, Ay, Bx, By, N, r) {
  d <- sapply(list(Ax, Ay, Bx, By), function(l) runif(N, l[1], l[2]))
  sum(sqrt((d[, 1] - d[, 3])^2 + (d[, 2] - d[, 4])^2) <= r) / N
}

set.seed(123)
MCI(c(-0.5, 0.5), c(-0.5, 0.5), c(1.5, 2.5), c(-0.5, 0.5), 100000, 1.25)
# [1] 0.0194
MCI(c(-0.5, 0.5), c(-0.5, 0.5), c(-0.5, 0.5), c(1.5, 2.5), 100000, 1.25)
# [1] 0.01929