如何计算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 package (v0.3.0) 于 2020-09-05 创建}