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?
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).