2011-10-28 19 views
6

Me gustaría realizar una integración numérica en una dimensión, donde el integrando tiene un valor vectorial de. integrate() solo permite integrales escalares, por lo que necesitaría llamarlo varias veces. El paquete cubature parece adecuado, pero parece funcionar bastante mal para integrales 1D. Consideremos el siguiente ejemplo (integrando con valores escalares y la integración 1D),rendimiento de adaptIntegrate vs. integrar

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 

alguna manera, adaptIntegrate parece estar usando muchas evaluaciones más función para una precisión similar. Ambos métodos aparentemente usan la cuadratura Gauss-Kronrod (caso 1D: regla de cuadratura gaussiana de 15 puntos), aunque ?integrate agrega un "algoritmo Epsilon de Wynn". ¿Eso explicaría la gran diferencia de tiempo?

Estoy abierto a sugerencias de formas alternativas de tratar integrandos con valores vectoriales, tales como

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 

Gracias.

+0

No lo siento; ¿Qué pasa con la comparación uno a uno que doy para el integrando de valor escalar? – baptiste

+0

Hice la prueba con 'fDim = 2' (último ejemplo, 345 evaluaciones también), la comparación es simplemente un caso de invocar' integrate' dos veces, 'str (lapply (c (integrand1, integrand2), integrate, -1,1 , rel.tol = tolerancia, subdivisiones = Nmax)) 'da 10 + 1 = 11 evaluaciones. Mi punto es, sí, 'adapteIntegrate' se dirige a la integración multidimensional, y opcionalmente a integrandos con valores vectoriales, pero el caso de la integración unidimensional es mucho menos eficiente que llamar' integra' repetidamente, pero un margen grande (~ 30 veces aquí). – baptiste

+0

¿Has visto este paquete: http://cran.r-project.org/web/packages/R2Cuba/ –

Respuesta

2

También hay R2Cuba paquete en CRAN que implementa varios algoritmos de integración multidimensional:

he intentado poner a prueba esta función con su ejemplo, y en tal caso sencillo que no podía conseguir todos los algoritmos para trabajar (aunque No lo intenté muy duro), y pocos métodos que obtuve para trabajar fueron considerablemente más lentos que adaptIntegrate con la configuración predeterminada, pero tal vez en su verdadera aplicación este paquete podría valer la pena intentarlo.