adaptIntegrate와의 통합 불일치

Question

최대 분산 거리를 가진 분산 커널 (거리의 함수)을 고려할 때 격자가있는 풍경을 따라 종의 분산 확률을 추정합니다. 나는 면적 - 면적 분산 확률을 식에 기술 된대로 계산하려고 시도하고있다. 8 of this (open access) paper. 이것은 소스 및 표적 세포의 각각의 가능한 소스 및 타겟 포인트의 조합에 대한 분산 함수의 값을 평가하는 4 중 적분을 포함한다.

I는 다음과 같이 셀 B 타겟팅 소스 셀 A를 들어, cubature 패키지 adaptIntegrate 이것을 구현하며 한 분산이 1 인 단순 분산 커널시 간 점 거리 별도로> 1.25 0. 이 그림은 아래에 그래픽으로 표시되어 있습니다. 셀 A의 점이 1.25 내에 있지 않기 때문에 셀 B의 빨간색 영역에 도달 할 수 없습니다.

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 

는 I 얻을 다른 적분 표적 세포를 고려하면, 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에 더 적합하다는 것을 알았지 만 코딩 자체가 간과되는 부분이있을 수 있으므로 여기에 게시했습니다.

---

편집

: integrate 중첩 된 A -> B 경우를 들어 는 첫 번째 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)의 메인테이너 내가 위의 질문에보고하고,이 버그가 될 가능성이 있음으로 Cubature C 라이브러리 같은 결과를 얻을 수 있음을 나에게 알렸다; 오히려, h- 적응 형 큐브 루틴 (R 패키지가 인터페이스 인)은 경우에 따라 Cubature의 p- 적응 루틴보다 정확하지 않습니다. 이는 적절한 지역에서 doubles the number of sampling points입니다.

Naras도 내 질문에 제시된 두 가지 경우에 대해 일관된 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)