calculating double integrals in R quickly

Question

I'm looking for a solution for a double integral that is faster than

integrate(function(y) { 
   sapply(y, function(y) {
     integrate(function(x) myfun(x,y), llim, ulim)$value
   })
 }, llim, ulim)

with eg

myfun <- function(x,y) cos(x+y)
llim <- -0.5
ulim <- 0.5

I found an old paper that referred to a FORTRAN program called quad2d, but I couldn't find anything else but some help pages for matlab for the rest. So I'm looking for a C or FORTRAN library that can do double integrals quick (i.e. without the sapply loop), and that can be called from R. All ideas are very much appreciated, as long as they're all GPL compatible.

If the solution involves calling other functions from the libraries that are already shipped with R, I'd love to hear from them as well.

Attending

Have you considered: pracma::dblquad, pracma:simpson2d, or the functions in the cubature and R2Cuba packages? – Depurative

Answer 1

The cubature package does 2D (and N-D) integration using an adaptive algorithm. It should outperform simpler approaches for most integrands.

Decern

Thx for the pointer, I'm running the tests this weekend and see which one is making me most happy. – Attending

For future reference, the cubature approach would be: adaptIntegrate(function(x) cos(x[1] + x[2]), c(llim, llim), c(ulim, ulim)). – Precognition

Hi, I am also struggling with a double integral with integrate which I find quite slow but I don't know understand how you specify the limits of integration with adaptIntegrate I want to integrate between 35 and u for the first integral and between 35 and 95 for the second, with u moving between 35 and 95. How would you write it ? Thanks – Maintopmast

Answer 2

The pracma package that Joshua pointed out contains a version of quad2d.

quad2d(myfun, llim, ulim, llim, ulim)

This gives the same answer, within numerical tolerance, as your loop, using the example function.

By default, with your example function, quad2d is slower than the loop. If you drop n down, you can make it faster, but I guess it depends upon how smooth your function is, and how much accuracy you are willing to sacrifice for speed.

Creech

Thx for the example as well. I'll test it on the real functions, as they're a bit more complex. I need accuracy the most, but speed really is an issue. I'll keep you updated. – Attending