CODEWITHSUNDEEP
r
Tou Mou
Tou Mou

Reputation: 1274

How to compute double integral in r?

Suppose that we have the following density :

  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)

I'm wanting to compute the follwing integral :

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

Unfortunately , it provides this error :

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

I don't know how to resolve this. Thank you for help in advance !

Upvotes: 2

Views: 942

Answers (2)

Rui Barradas
Rui Barradas

Reputation: 76460

With package cubature, functions hcubature and pcubature the integrand would have to be changed a bit. The integrators from that package accept integrand functions of one variable only, that can be a vector in a multidimensional real space. In this case, R2. The values of x and y would have to be assigned in the integrand or change to become x[1] and x[2] in its expression.

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)

Upvotes: 4

Allan Cameron
Allan Cameron

Reputation: 173888

If you need to do a double integral, you could just integrate twice:

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

This gives an answer that is pleasingly close to 1, as expected.

Created on 2020-09-05 by the reprex package (v0.3.0)

Upvotes: 2

Related Questions