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.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With