How to compute double integral in r?

Question

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 !

The function `f2` produced by `bvtnorm` takes two parameters, x and y, but `adaptIntegrate` expects a function of a single variable, so it can't be used for a double integral in this way — Allan Cameron, Sep 05 '20 at 18:41
[double integral in R](https://stackoverflow.com/q/43189512/4752675) may help. — G5W, Sep 05 '20 at 18:58

score 4 · Answer 1 · answered Sep 05 '20 at 19:29

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)

This is nice Rui. I don't know if the lower integration limits are intentionally set to `c(-Inf, 0)`, in the question, but if they are changed to `c(-Inf, -Inf)` this returns an integral of exactly 1, so I think this is the better answer. — Allan Cameron, Sep 05 '20 at 19:55
@AllanCameron That is due to a default `rel.tol = .Machine$double.eps^0.25` (0.0001220703), try my `rel.tol = eps` (0.0000000149) and the integral is 1. — Rui Barradas, Sep 05 '20 at 20:00
@RuiBarradas , i will use your answer because it requires a minor modification compared to the first one ! Thank you both for your valuable help ! — Tou Mou, Sep 05 '20 at 20:23

score 2 · Accepted Answer · answered Sep 05 '20 at 19:11

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)}

Yes this help a lot ! I already tried to use "rmutil" package but it's very slow when bounds are equal to Inf. In the opposite , your solutions seems usefull for me. Thank you a lot again ! — Tou Mou, Sep 05 '20 at 20:19

How to compute double integral in r?

2 Answers2