rprecisionnumerical-integration

Precision issue of integrate function in R


When using the integrate() function in R, some precision issues make me confused. For example, I want to integrate this function:

p <- function(x)exp(-x^2)

I choose the interval [-1e3,1e3], and it gives the right answer:

> integrate(p, -1e3, 1e3)
1.772454 with absolute error < 7.8e-09

But if I expand the interval to [-1e9,1e9], it's wrong:

> integrate(p, -1e9, 1e9)
0 with absolute error < 0

In order not to tune the interval, I try to use Inf:

> integrate(p, -Inf, Inf)
1.772454 with absolute error < 4.3e-06

It works, but soon I find there are some other problems :

> integrate(p, -Inf, -Inf)
1.772454 with absolute error < 4.3e-06

In this case it can't give the right answer 0.

I wonder why these precision problem will happen and are there any methods to avoid such problem?


Solution

  • You are advised to use Inf,-Inf. Also it seems, interestingly, that integrate() is converting -Inf to Inf when used as upper limit. :

    > integrate(p, -Inf, -1*.Machine$double.xmax)
    0 with absolute error < 0
    > integrate(p, -Inf, -2*.Machine$double.xmax)
    1.772454 with absolute error < 4.3e-06
    

    This didn't result in what we would expect. Lets try to split the integral in two, first:

    > integrate(p, 0, Inf)
    0.8862269 with absolute error < 2.2e-06
    > integrate(p, 0, 1e100)
    0 with absolute error < 0
    > integrate(p, 0, 1e2)
    0.8862269 with absolute error < 9.9e-11
    

    This seems perfectly consistent with the advice of using Inf,-Inf, but notice what happens, when switching signs:

    > integrate(p, 0, -Inf)
    0.8862269 with absolute error < 2.2e-06
    > integrate(p, 0, -1e100)
    0 with absolute error < 0
    > integrate(p, 0, -1e2)
    -0.8862269 with absolute error < 9.9e-11
    

    At last you can always change tolerance, subdivisions, but this will not help you in this case.

    EDIT: As @Ben pointed out, it's a minor bug because integrate checks whether the upper limit is finite and not if it's Inf (or -Inf).