[R] bivariate vector numerical integration with infinite range

baptiste auguie baptiste.auguie at googlemail.com
Tue Sep 21 22:16:50 CEST 2010


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?

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

On 21 September 2010 17:11, Hans W Borchers <hwborchers at googlemail.com> wrote:
> baptiste auguie <baptiste.auguie <at> googlemail.com> writes:
>
>>
>> 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.
>
> Yes, that's one way.
> But you can also split the integral in two parts, one from 0 to 1 and then
> from 1 to Inf. The first one is a finite hypercube and the second can be
> transformed with x --> 1/x into [0, 1].
>
> I usually prefer the second approach for higher-dimensional applications
> as the Jacobian appears to be simpler. In the literature you will find
> discussions on how far out the finite hypercube should reach for lowering
> the absolute error.
>
> Hans Werner
>
>> Thanks,
>>
>> 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.
>



More information about the R-help mailing list