[R] bivariate vector numerical integration with infinite range

[David Winsemius]

Baptiste;

You should see if this meets your requirements:

help(adaptIntegrate, package=cubature)

(I got errors when I ran the code and NaN's when I looked at the  
output of test function, "f".)

 > vAverage(mixedrule, -4, 4, 0.0, 1, 20, f) - c(pi, pi/2, 2*pi)
Error: object 'mixedrule' not found

 >  f(-4,0)
[1] NaN NaN NaN
Warning message:
In sqrt(x) : NaNs produced
 >  f(-3.9,0)
[1] NaN NaN NaN
Warning message:
In sqrt(x) : NaNs produced
 >  f(-3.9,0.1)
[1] NaN NaN NaN
Warning message:
In sqrt(x) : NaNs produced
-- 
David
On Sep 21, 2010, at 4:11 AM, baptiste auguie wrote:

> 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 at 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.

David Winsemius, MD
West Hartford, CT