R help - Sep 2010 - bivariate vector numerical integration with infinite range

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

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

baptiste auguie <baptiste.auguie <at> googlemail.com> writes:
> 
> 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,
The function adaptIntegral() in the 'cubature' package integrates 
multi-valued functions over n-dimensional finite hypercubes, as do the
functions in 'R2Cuba'.
If the hypercube is (partly) infinite, a transformation such as x --> 1/x
per infinite axis (as in NR) has to be applied.

For two dimensions, another approach could be to apply the integrate()
function twice, because this 1-dimensional integration function can handle
infinite intervals.
Hint: First inegrate over the finite dimension(s).

Example: Integrate sin(x)/exp(y) for 0 <= x <= pi, 0 <= y <= Inf
(value: 2)
----
    f1 <- function(y) 1/exp(y)
    f2 <- function(x) sin(x) * integrate(f1, 0, Inf)$value
    integrate(f2, 0, pi)
    # 2 with absolute error < 2.2e-14
----
Note that the absolute error is not correct, as the first integration has
its own error term. You have to do your own error estimation.

Hans Werner
> 
> [...]
> 
> 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
> 
>

Got it, thanks!

baptiste

On 21 September 2010 22:38, Hans W Borchers <hwborchers at googlemail.com>
wrote:
> baptiste auguie <baptiste.auguie <at> googlemail.com> writes:
>
>>
>> Thanks. I am having trouble getting adaptIntegrate to work with a
>> multivalued integrand though, and cannot find a working example.
>> Anyone had better luck with it?
>
> The function to be integrated needs a vector as input:
>
> ? ?f <- function(x) {
> ? ? ? ?res <- 1 / (sqrt(x[1])*(1+x[1]))
> ? ? ? ?c(res, res/2, 2*res)
> ? ?}
>
> ? ?adaptIntegrate(f, lowerLimit=c(0.1, 0), upperLimit=c(10, 1), fDim = 3)
> ? ?$integral
> ? ?[1] 1.9164832 0.9582416 3.8329665
>
> ? ?$error
> ? ?[1] 1.265252e-05 6.326261e-06 2.530504e-05
>
> ? ?$functionEvaluations
> ? ?[1] 323
>
> ? $returnCode
> ? [1] 0
>
> Hans Werner
>
>> library(cubature)
>> >
>> > f <- function(x, y) {
>> + ? res <- 1 / (sqrt(x)*(1+x))
>> + ? c(res, res/2, 2*res)
>> + }
>> >
>> > adaptIntegrate(f, lowerLimit=c(0.1, 0), upperLimit=c(10, 1), fDim
= 3)
>> [1] "adaptIntegrate: Error in evaluation function f(x) for
x="
>> ? ? ? ? ? ?res
>> [1,] 0.07355275 0.03677638 0.1471055
>> [2,] 0.94280904 0.47140452 1.8856181
>> Error in adaptIntegrate(f, lowerLimit = c(0.1, 0), upperLimit = c(10,
?:
>> ?adaptIntegrate: Result f(x) is not numeric or has wrong dimension
>>
>> Best,
>>
>> 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.
>