[R] two dimensional integration

[Berend Hasselman]

On 16-02-2013, at 18:01, julia cafnik <julia.cafnik at gmail.com> wrote:

> Dear R-users,
> 
> I'm wondering how to calculate this double integral in R:
> 
> int_a^b  int_c^y g(x, y) dx dy
> 
> where g(x,y) = exp(- alpha (y - x)) * b
> 

A very similar question was asked about nine years ago: http://tolstoy.newcastle.edu.au/R/help/04/10/5831.html
The final message in the thread gave a workable answer.

In your case define function g (leave out the multiply by b since you can always do that outside the integral).

g <- function(x,y,alpha=1) exp(-alpha*(y-x))

Assume a=0, b=1 and c=0.

Then following the final message in the thread mentioned above there are two ways of getting an answer:

if function g is fully vectorized:

integrate( function(y) {
    sapply(y, function(y) { 
        integrate(function(x) g(x,y),0,y)$value
    })
},0,1)

and if g is not vectorized and only takes scalars as input

integrate(function(y) {
  sapply(y, function(y) {
    integrate(function(x) {
      sapply(x, function(x) g(x,y)) },0,y)$value
  })
},0,1)

Both of the approaches result in

0.3678794 with absolute error < 4.1e-15

which corresponds to the exact analytical answer 1/exp(1) (if I didn't make a mistake)

You can also do this with package cubature with Gabor's approach (from the mentioned thread) like this

library(cubature)              
h <- function(z) g(z[1],z[2])*(z[1]<z[2])
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4) 

but it's also very slow with higher tolerances.

> adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)
$integral
[1] 0.3678723
$error
[1] 3.677682e-05
$functionEvaluations
[1] 268413
$returnCode
[1] 0


Berend