时间戳

首页 > 解决方案 > 使用 hcubature 在 R 中集成

r - 使用 hcubature 在 R 中集成

问题描述

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

标签： rintegration

解决方案

这个问题对于任何数值积分方法都是常见的：如果你没有在正确的地方评估函数，它就会看起来是恒定的。您的关节密度以 附近为中心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。

推荐阅读