[R] bivariate vector numerical integration with infinite range

[baptiste auguie]

Thanks, adaptIntegrate() seems perfectly suited, I'll just need to
figure a transformation rule for the infinite limit. The suggestion of
x->1/x does not seem to work here because it also transforms 0 into
-infinity. I think exp(pi* sinh(x)) could be a better choice,
according to Numerical Recipes.

Thanks,

baptiste


On 21 September 2010 14:26, Hans W Borchers <hwborchers at googlemail.com> wrote:
> 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
>>
>>
>
> ______________________________________________
> 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.
>