我想在一维中执行数值积分,其中被积分数是矢量值integrate()仅允许标量被整数,因此我需要多次调用它。 cubature包似乎很适合,但对于一维积分似乎效果很差。考虑以下示例(标量值被积和一维积分),

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

某种程度上,adaptIntegrate似乎正在使用更多的函数求值以达到类似的精度。尽管?integrate添加了“Wynn's Epsilon算法”,但两种方法显然都使用高斯-克朗罗德正交(1D情况:15点高斯正交规则)。那可以解释时间上的巨大差异吗?

我乐于接受向量值积分对象的替代方法的建议,例如
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

谢谢。

最佳答案

CRAN中还有R2Cuba软件包,它实现了几种多维集成算法:

我尝试使用您的示例函数对此进行测试,在这种简单情况下,我无法使所有算法都起作用(尽管我并未尽力而为),而且我上手的几种方法都比adaptIntegrate慢得多使用默认设置,但也许在您的真实应用程序中,此程序包值得尝试。

09-18 23:08