使用 hcubature 集成到 R 中

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，但正确答案应该是 0.1348。谁能帮帮我？非常感谢。

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。

使用 hcubature 集成到 R 中

Integration in R using hcubature

integration

r