如何计算r中的二重积分

Question

假设我们有以下密度：

bvtnorm <- function(x, y, mu_x = 10, mu_y = 5, sigma_x = 3, sigma_y = 7, rho = 0.4) {

function(x, y) 
1 / (2 * pi * sigma_x * sigma_y * sqrt(1 - rho ^ 2)) * 
exp(- 1 / (2 * (1 - rho ^ 2)) * (((x - mu_x) / sigma_x) ^ 2 + 
((y - mu_y) / sigma_y) ^ 2 - 2 * rho * (x - mu_x) * (y - mu_y) / 
(sigma_x * sigma_y)))
}

f2 <- bvtnorm(x, y)

我想计算以下积分：

integral_1=1-adaptIntegrate(f2, lowerLimit = c(-Inf,0), upperLimit = c(+Inf,+Inf))

不幸的是，它提供了这个错误：

Error in f(tan(x), ...) : argument "y" is missing, with no default

我不知道如何解决这个问题。提前感谢您的帮助！

Answer 1

对于包cubature、函数hcubature和pcubature，被积函数必须稍微更改。该包中的积分器只接受一个变量的被积函数，该被积函数可以是多维实空间中的向量。在这种情况下，R2。x和y的值必须在被积函数中赋值，或者在其表达式中变为x[1]和x[2]。

bvtnorm <- function(x, mu_x = 10, mu_y = 5, sigma_x = 3, sigma_y = 7, rho = 0.4) {

y <- x[2]
x <- x[1]
1 / (2 * pi * sigma_x * sigma_y * sqrt(1 - rho ^ 2)) * 
exp(- 1 / (2 * (1 - rho ^ 2)) * (((x - mu_x) / sigma_x) ^ 2 + 
((y - mu_y) / sigma_y) ^ 2 - 2 * rho * (x - mu_x) * (y - mu_y) / 
(sigma_x * sigma_y)))
}
library(cubature)
eps <- .Machine$double.eps^0.5
hcubature(bvtnorm, lowerLimit = c(-Inf, 0), upperLimit = c(+Inf,+Inf), tol = eps)
pcubature(bvtnorm, lowerLimit = c(-Inf, 0), upperLimit = c(+Inf,+Inf), tol = eps)

Answer 2

如果你需要做二重积分，你只需要integrate两次：

bvtnorm <- function(y, mu_x = 10, mu_y = 5, sigma_x = 3, sigma_y = 7, rho = 0.4) {

function(x) 
1 / (2 * pi * sigma_x * sigma_y * sqrt(1 - rho ^ 2)) * 
exp(- 1 / (2 * (1 - rho ^ 2)) * 
(((x - mu_x) / sigma_x) ^ 2 + 
((y - mu_y) / sigma_y) ^ 2 - 2 * rho * (x - mu_x) * (y - mu_y) / 
(sigma_x * sigma_y)))
}

f3 <- function(y)
{
f2 <- bvtnorm(y = y)
integrate(f2, lower = -Inf, upper = Inf)$value
}

integrate(Vectorize(f3), -Inf, Inf)
#> 1.000027 with absolute error < 1.8e-05

正如预期的那样，这给出了一个令人愉快的接近1的答案。

^{由reprex包(v0.3.0(于2020-09-05创建}