Dear list,
I'm seeking some advice regarding a particular numerical integration I
wish to perform.
The integrand f takes two real arguments x and y and returns a vector
of constant length N. The range of integration is [0, infty) for x and
[a,b] (finite) for y. Since the integrand has values in R^N I did not
find a built-in function to perform numerical quadrature, so I wrote
my own after some inspiration from a post in R-help,
library(statmod)
## performs 2D numerical integration
## using Gauss-Legendre quadrature
## with N points for x and y
vAverage <- function(f, a1,b1, a2,b2, N=5, ...){
GL <- gauss.quad(N)
nodes <- GL$nodes
weights <- GL$weights
C2 <- (b2 - a2) / 2
D2 <- (b2 + a2) / 2
y <- nodes*C2 + D2
C1 <- (b1 - a1) / 2
D1 <- (b1 + a1) / 2
x <- nodes*C1 + D1
value <- 0
for (ii in seq_along(x)){
tmp <- 0
for (jj in seq_along(y)){
tmp <- tmp + C1 * weights[jj] * f(x[jj], y[ii], ...)
}
value <- value + C2 * weights[ii] * tmp
}
value
}
## test function, the result is pi for y=1
f <- function(x, y) {
res <- 1 / (sqrt(x)*(1+x))
c(res, res/2, 2*res)
}
## Transformation rule from Numerical Recipes
## to deal with the [0, infty) range of x
mixedrule <- function(x, y, f, ...)
{
t <- exp(pi*sinh(x))
dtdx <- t*(pi*cosh(x))
f(t, y, ...)*dtdx
}
vAverage(mixedrule, -4, 4, 0.0, 1, 20, f) - c(pi, pi/2, 2*pi)
## -3.535056e-06 -1.767528e-06 -7.070112e-06
So it seems to work. I wonder though if I may have missed an easier
(and more reliable) way to perform such integration using base
functions or an add-on package that I may have overlooked.
Best regards,
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.
