[R] bivariate vector numerical integration with infinite range

[baptiste auguie]

Hi,

thanks for the pointer to cubature (which i had probably dismissed too
quickly). Your tests with "f" should not work: the domain of f(x,.) is
restricted to positive reals, but this domain of integration is then
transformed in mixedrule() to map the semi-infinite range to a more
reasonable domain (for example [-4,4]).

Thanks,

baptiste




On 21 September 2010 14:16, David Winsemius <dwinsemius at comcast.net> wrote:
> 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
>
>