Double Integral in R

How to compute double integral in r?

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

How to calculate the double integration in R

Here is how I would do it (for one version of beta and X). Note that a double integral is just two single integrals nested into one another. A parameter n defines the number of random samples I am using for estimating the expectation in the integral.

beta <- function(t){
    return(t*t*t-1.6*t*t+0.76*t+1)
}

myX <- function(a,t){
    pt <- c(1,t,t*t,t*t*t)
    return(sum(a*pt))
}

## computes the expectation by averaging over n samples
myE <- function(n,s,t){
    samp <- sapply(seq(n),function(x){
          a <- rnorm(4)
          myX(a,s)*myX(a,t)})
    return(mean(samp,na.rm=T))
}

## funtion inside the first integral
myIntegrand1 <- function(s,t,n){
   return(beta(s)*myE(n,s,t))
}

## function inside the second integral
myIntegrand2 <- function(t,n){
    v <- integrate(myIntegrand1,0,1,t=t,n=n)
    return(beta(t)*v$value)
}

## computes sigma
mySig <- function(n){
  v <- integrate(myIntegrand2,0,1,n=n)
  return( 0.25*v$value)
}
## tests various values of n (number of samples drawn to compute the expectation)
sapply(seq(3),function(x)
             c("100"=mySig(100),"1000"=mySig(1000),"10000"=mySig(10000)))
## output shows you the level of precision you may expect:
##       [,1]     [,2]     [,3]
##  100  48.61876 47.85445 58.2094
## 1000  52.95681 50.61860 50.61702
## 10000 54.88292 53.02073 54.48635

How to double integrate with +inf?

Use the cubature package. It evaluates multiple integrals, allowing infinite bounds.

Here is an example:

library(cubature)
f <- function(x) exp(-x[1]-x[2])
pcubature(f, c(0,0), c(Inf,Inf))$integral
# 1

library(pracma)
f2 <- function(x,y) f(c(x,y))
integral2(f2, 0, 1000, 0, 1000, vectorized = FALSE)$Q
# 1

Double Integral implementation in R

First I needed to convince myself that Mathematica was correct. Yeah, I suppose that stems from my distrust of authority, but it was not difficult. It required realizing that the integral of unity from upper to lower was simply the difference of upper minus lower, so it could be reduced to a single variable problem and solved with R's integrate function:

integrate( function(x){sin(x)-0.5},lower=pi/6,upper=5*pi/6)
0.6848533 with absolute error < 7.6e-15

So that led me to try the same strategy in the 'mosaicCalc' framework:

> one = makeFun(1 ~ y + x)
> bx.y = antiD(one(y = y, x = x) ~ y)
> bx.yx = antiD(bx.y(y = sin(x)-1/2, x = x) ~ x)
> bx.yx(x = 5*pi/6) - bx.yx(x = pi/6)
[1] 0.6848533

That got the right answer but it didn't seem to properly represent and preserve the "functional cascade" (to invent a term I've never heard used before). I wanted the limits to appear in a manner that reflected a more general functional set of calls, so eventually came up with this which seems satisfactory:

one = makeFun(1 ~ y + x)
bx.y = antiD(one(y = y, x = x) ~ y)
bx.yx = antiD(bx.y(x=x, y = sin(x)) - bx.y(x=x,y=1/2) ~ x)
bx.yx(x = 5*pi/6) - bx.yx(x = pi/6)
[1] 0.6848533

This does support more complex functional integrands. I don't have Mathematica set up to do that integration with anything but with the mosaicCalc setup above I get 1.284286 for x^2 as an integrand function. You might want to check.

In working this problem it did seem to me that the order of integration was reversed. In my hazy memory of these problems from 40, no 50 years ago, it seemed I had always used the dx calculation as the inner one, but I do realize that is arbitrary. At any rate the roles you assigned to the x and y values in the second anti-derivative didn't seem appropriate. You were getting the result of a 2D integration of unity from lower limits of x=pi/6 and y=0.5 to upper limits of x=5*pi/6, y=1)

library(cubature) # package capable of 2D-integration with fixed limits
adaptIntegrate(function(x){1}, lower=c(a=pi/6, b=0.5), upper=c(a=5*pi/6,b=1 ))
$integral
[1] 1.047198

$error
[1] 0

$functionEvaluations
[1] 17

$returnCode
[1] 0

Why is my double integral is R not working

Your variables are all out of wack. This should do what you want

InnerFunc <- function(x, y) { x + (y^2) }
InnerIntegral <- function(y) { sapply(y, 
  function(z) { integrate(InnerFunc, 0, 2, y=z)$value }) }
integrate(InnerIntegral, 0, 1)
# 2.666667 with absolute error < 3e-14

Notice how InnerFunc is now a function of both x and y. Also notice how we apply over the vector passed to InnerIntegral via the y parameter (rather than ignoring the z value passed in). We also pass in the current value of y to InnerFunc.

Although you typed

InnerFunc <- function(x, y) { x + (y^2) }

the math you drew looks should really result in

InnerFunc <- function(x, y) { x * (y^2) }

so I'm not sure what you were really after.

---

Related Topics