2011-10-28 22 views
6

Mi piacerebbe eseguire un'integrazione numerica in una dimensione, dove l'integrando è valore vettoriale. integrate() consente solo integram scalari, quindi dovrei chiamarlo più volte. Il pacchetto cubature sembra molto adatto, ma sembra non funzionare abbastanza bene per gli integrali 1D. Si consideri il seguente esempio (integrando valori scalari e integrazione 1D),performance di adaptIntegrate vs. integra

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 

qualche modo, adaptIntegrate sembra utilizzare molte valutazioni più funzionali per una precisione simile. Entrambi i metodi usano apparentemente la quadratura di Gauss-Kronrod (caso 1D: regola di quadratura gaussiana di 15 punti), sebbene ?integrate aggiunga un "algoritmo di Epsilon di Wynn". Questo spiegherebbe la grande differenza temporale?

Sono aperto a suggerimenti di modi alternativi di trattare con integrandi valori vettoriali quali

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 

Grazie.

+0

Non mi dispiace; cosa c'è di sbagliato nel confronto one-to-one che ho per integrando con valori scalari? – baptiste

+0

Ho eseguito il test con 'fDim = 2' (ultimo esempio, 345 anche le valutazioni), il confronto è solo un caso di chiamare' integrare' due volte, 'str (lapply (c (integrand1, integrand2), integra, -1,1 , rel.tol = tolerance, subdivisions = Nmax)) 'dà 10 + 1 = 11 valutazioni. Il mio punto di vista, sì, 'adaptIntegrate' è rivolto all'integrazione multidimensionale e facoltativamente agli integrandi con valori vettoriali, ma il caso dell'integrazione unidimensionale è molto meno efficiente di chiamare' integrarsi 'ripetutamente, ma un ampio margine (~ 30 volte Qui). – baptiste

+0

Hai visto questo pacchetto: http://cran.r-project.org/web/packages/R2Cuba/ –

risposta

2

C'è anche R2Cuba pacchetto CRAN che implementa diversi algoritmi di integrazione multidimensionali:

ho provato a verificare questa con la funzione esempio, in un caso come semplice che non poteva ottenere tutti gli algoritmi di lavoro (sebbene Non ci ho provato molto duramente), e pochi metodi che ho avuto modo di lavorare sono stati molto più lenti di adaptIntegrate con le impostazioni predefinite, ma forse nella tua vera applicazione questo pacchetto potrebbe valere la pena di provarci.