與adaptIntegrate的積分不一致性

Question

我估計了一個物種在網格化景觀中的散佈概率，給定了散佈核心（距離函數）和最大散佈距離。我試圖計算如公式1中所述的面積 - 面積擴散概率。 8的this (open access) paper。這涉及四倍積分，分別評估源和目標單元中源點和目標點的每種可能組合的擴散函數的值。

我已經與adaptIntegrate從cubature封裝中實現這一點，如下，源小區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的權

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 vs 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)

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

當這樣一定是有什麼功能裏面發生了，因爲值的範圍顯着不同是r值。

對此的一種替代方法是蒙特卡羅積分，在這種情況下，這是很好的，因爲你的點是均勻分佈的。

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 

Answer 2

一般距離的措施是(x1-x2)^2+(y1-y2)^2。你能解釋一下爲什麼在構造r時從y中減去x的值？考慮修改後的版本：

f <- function(xmin, xmax, ymin, ymax) { 
    adaptIntegrate(function(x) { 
    r <- sqrt((x[4] - x[3])^2 + (x[2] - x[1])^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.01016105 

$error 
[1] 0.0241325 

$functionEvaluations 
[1] 100035 

$returnCode 
[1] 0 
#--------- 
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 

Answer 3

將R cubature包（NARAS）的維護者已通知我，該數值積分C庫給出了相同的結果，因爲我在這個問題上面的報告，並認爲這是不太可能是一個錯誤;相反，h-自適應容積式例程（R包是一個接口）在某些情況下比Cubature的p-自適應例程（在合適的區域中爲doubles the number of sampling points）要精確。

Naras還提供了以下julia代碼，該代碼證明了我的問題中呈現的兩種情況（返回值的元素是估計的積分，然後是估計的絕對誤差）的一致的pcubature解決方案。

using Cubature 

# integrand 
f = x -> ifelse(sqrt((x[3] - x[1])^2 + (x[4] - x[2])^2) > 1.25, 0, 1) 

# A to B case 
pcubature(f, [-0.5, -0.5, 1.5, -0.5], [0.5, 0.5, 2.5, 0.5], abstol=1e-5)  
# (0.019593408732917292,3.5592555263398717e-6) 

# A to C case 
pcubature(f, [-0.5, -0.5, -0.5, 1.5], [0.5, 0.5, 0.5, 2.5], abstol=1e-5) 
# (0.019593408732918302,3.559255527241928e-6)