使用 hcubature 在 R 中集成答案

【问题标题】：Integration in R using hcubature使用 hcubature 在 R 中集成
【发布时间】：2020-04-20 15:21:00
【问题描述】：

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。有人可以帮我吗？非常感谢。

【问题讨论】：

rho的值是多少？
对不起，rho
尝试同一个包中的另一个函数：cuhre(f1, 1L, c(-Inf, -Inf), c(Inf, Inf))。这给了[1] 0.1347824。
但是Matlab给出0.1348，R使用积分函数也给出0.1348。有什么问题？非常感谢。
如何在R中选择合适的函数进行集成？不同的功能给出不同的答案。

标签： r integration

【解决方案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。

【讨论】：