在R中使用hcubature进行集成

Question

rho <- 0.2
f1 <- function(X) { 
  x <- X[1] 
  y <- X[2] 
  (pnorm(x-10)-pnorm(y-10))^2*(exp(-((x-10)^2-2*rho*(x-10)*(y-10)+(y-10)^2)/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2))) 
} 
library(cubature) 
round(hcubature(f1, c(-Inf, -Inf), c(Inf, Inf), tol = 1e-12)$integral, 6) 

这给出了0.1348，但正确答案应该是0。有人能帮帮我吗？非常感谢。

Answer 1

这个问题在任何数值积分方法中都是常见的:如果你不在正确的地方计算函数，它将看起来是恒定的。关节密度集中在c(10, 10)附近，而hcubature函数主要在c(0, 0)附近求值，因此它看起来是恒定的，接近0。您可以看到如下所示：

rho <- 0.2

# First, plot the density function

fn <- function(x, y) (pnorm(x-10)-pnorm(y-10))^2*(exp(-((x-10)^2-2*rho*(x-10)*(y-10)+(y-10)^2)/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2)))

# Record the points used by hcubature

pts <- matrix(nrow = 0, ncol = 2)
f1 <- function(X) { 
  pts <<- rbind(pts, X)
  fn(X[1], X[2])
}
library(cubature)
round(hcubature(f1, c(-Inf, -Inf), c(Inf, Inf), tol = 1e-12)$integral, 6)
#> [1] 0

# Draw the points on the plot
plot(pts, type = "p")

# Now show contours of the function
x <- y <- seq(min(pts), max(pts), length.out = 100)
z <- outer(x, y, fn)
contour(x, y, z, col = "red", add = TRUE)

​

​

如果你把函数放在接近0的位置，那就没问题了。使用以下命令运行上面的代码

fn <- function(x, y) (pnorm(x)-pnorm(y))^2*(exp(-((x)^2-2*rho*(x)*(y)+(y)^2)/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2)))

并将结果打印为0.134783。