r - 集成与Adapt不一致

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

我已经使用adaptIntegrate包中的cubature实现了此操作，如下所示，针对源单元A，目标单元B和简化的分散核，其中，当点间距离> 1.25时，分散为1，否则为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时，我得到一个不同的积分，该目标单元格C的距离相同，但位于单元格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

最佳答案

这样做的原因似乎是adaptIntegrate的工作方式，因为显然，您唯一需要更改的就是集成的顺序。可能由于单独的近似积分而导致不相同的结果(请参阅第一个响应here)，但这似乎更像是一个错误。

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

和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