Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

performance of adaptIntegrate vs. integrate

I'd like to perform a numerical integration in one dimension, where the integrand is vector-valued. integrate() only allows scalar integrands, thus I would need to call it several times. The cubature package seems well suited, but it seems to perform quite poorly for 1D integrals. Consider the following example (scalar-valued integrand and 1D integration),

library(cubature)
integrand <- function(x, a=0.01) exp(-x^2/a^2)*cos(x)
Nmax <- 1e3
tolerance <- 1e-4

# using cubature's adaptIntegrate
time1 <- system.time(replicate(1e3, {
  a <<- adaptIntegrate(integrand, -1, 1, tol=tolerance, fDim=1, maxEval=Nmax)
}) )

# using integrate
time2 <- system.time(replicate(1e3, {
  b <<- integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax)
}) )

time1
user  system elapsed 
  2.398   0.004   2.403 
time2
user  system elapsed 
  0.204   0.004   0.208 

a$integral
> [1] 0.0177241
b$value
> [1] 0.0177241

a$functionEvaluations
> [1] 345
b$subdivisions
> [1] 10

Somehow, adaptIntegrate seems to be using many more function evaluations for a similar precision. Both methods apparently use Gauss-Kronrod quadrature (1D case: 15-point Gaussian quadrature rule), though ?integrate adds a "Wynn's Epsilon algorithm". Would that explain the large timing difference?

I'm open to suggestions of alternative ways of dealing with vector-valued integrands such as

integrand <- function(x, a = 0.01) c(exp(-x^2/a^2), cos(x))
adaptIntegrate(integrand, -1, 1, tol=tolerance, fDim=2, maxEval=Nmax)
$integral
[1] 0.01772454 1.68294197

$error
[1] 2.034608e-08 1.868441e-14

$functionEvaluations
[1] 345

Thanks.

like image 487
baptiste Avatar asked Oct 28 '11 03:10

baptiste


1 Answers

There is also R2Cuba package in CRAN which implements several multidimensional integration algorithms:

I tried to test this with your example function, and in such a simple case I couldn't get all the algorithms to work (although I didn't try really hard), and few methods I did get to work were considerably slower than adaptIntegrate with default setting, but perhaps in your true application this package could be worth of trying.

like image 176
Jouni Helske Avatar answered Nov 18 '22 08:11

Jouni Helske