Dear list,
[cross-posting from Stack Overflow where this question has remained
unanswered for two weeks]
I'd like to perform a numerical integration in one dimension,
I = int_a^b f(x) dx
where the integrand f: x in IR -> f(x) in IR^p is vector-valued.
integrate() only allows scalar integrands, thus I would need to call
it many (p=200 typically) times, which sounds suboptimal. The cubature
package seems well suited, as illustrated below,
library(cubature)
Nmax <- 1e3
tolerance <- 1e-4
integrand <- function(x, a=0.01) c(exp(-x^2/a^2), cos(x))
adaptIntegrate(integrand, -1, 1, tolerance, 2, max=Nmax)$integral
[1] 0.01772454 1.68294197
However, adaptIntegrate appears to perform quite poorly when compared
to integrate. Consider the following example (one-dimensional
integrand),
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, tolerance, 1, max=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 was using many more function evaluations for a
similar precision. Both methods apparently use Gauss-Kronrod
quadrature, though ?integrate adds a "Wynn's Epsilon algorithm". Could
that explain the large timing difference?
I'm open to suggestions of alternative ways of dealing with
vector-valued integrands.
Thanks.
baptiste
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
