Search code examples
rintegral

Multiple integration with functions of variables in the limits

I need to numerically integrate the following:

enter image description here

I tried to use cubature and pracma but they don't seem to support functional integration limits. I found a attempt to use cubature by:

library(cubature)

integrand <- function(arg) {
   x <- arg[1]
   y <- arg[2]
   z <- arg[3]
   w <- arg[4]
   v<- arg[5]
   ff <- dnorm(x, 10,2)*dnorm(y, 10,2)*dnorm(z, 10,2)*dnorm(w, 10,2)* dnorm(v, 10,2)* (x+y+z+w+v<=52)
   return(ff)
}

R <- cuhre(f = integrand, 
           lowerLimit=c(0,0,0,0,0), 
           upperLimit=c(20,20,20,20,20), 
           relTol = 1e-5, absTol= 1e-5)

But the author doesn't guarantee that it's correct to do it.

Is there a way to numerically integrate multiple integrals with functional limits in R?

Solution

  • The domain of integration is the canonical simplex scaled by the factor 42. To evaluate an integral on a simplex, use the SimplicialCubature package:

    integrand <- function(arg) {
  x <- arg[1]
  y <- arg[2]
  z <- arg[3]
  w <- arg[4]
  v <- arg[5]
  dnorm(x, 10, 2) * 
    dnorm(y, 10, 2) * 
    dnorm(z, 10, 2) * 
    dnorm(w, 10, 2) * 
    dnorm(v, 10, 2)
}

library(SimplicialCubature)
Simplex <- 42 * CanonicalSimplex(5)

    Here is the command to run: 

    adaptIntegrateSimplex(integrand, S = Simplex)
# $integral
# [1] 0.03252553
# 
# $estAbsError
# [1] 0.3248119
# 
# $functionEvaluations
# [1] 9792
# 
# $returnCode
# [1] 1
# 
# $message
# [1] "error: maxEvals exceeded - too many function evaluations"

    The algorithm has reached the maximal number of function evaluations and the estimated absolute error is 0.3248119, while the estimated value of the integral is 0.03252553. This is a big error.

    We can increase the maximum number of function evaluations allowed. Taking 1e6, the computation is a bit slow and we get:

    adaptIntegrateSimplex(integrand, S = Simplex, maxEvals = 1e6)
# $integral
# [1] 0.03682535
# 
# $estAbsError
# [1] 0.001004083
# 
# $functionEvaluations
# [1] 999811
# 
# $returnCode
# [1] 1
# 
# $message
# [1] "error: maxEvals exceeded - too many function evaluations"

    The estimated error has decreased to 0.001004083, quite better.

    Note that we can approximate this integral by using simulations, because this integral is the measure of the simplex under a multivariate normal distribution:

    library(mvtnorm)
Sigma <- 2^2 * diag(5)
Mean <- rep(10, 5)
set.seed(666)
sims <- rmvnorm(1e6, mean = Mean, sigma = Sigma)

f <- function(X){ # test whether 0 < x < 42, 0 < x + y < 42, 0 < x + y + z < 42, ...
  all(X > 0 & cumsum(X) < 42)
}

mean(apply(sims, 1, f))
# 0.037083