[R] bivariate vector numerical integration with infinite range

[baptiste auguie]

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