2

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?

qqwqert007
  • 125
  • 1
  • 1
  • 8
  • Increase the `subdivisions`, decrease the `rel.tol` and if you have enough memory you'll get enough precision. – Khashaa Jan 06 '15 at 14:23
  • 1
    See the note in the help file: `Like all numerical integration routines, these evaluate the function on a finite set of points. If the function is approximately constant (in particular, zero) over nearly all its range it is possible that the result and error estimate may be seriously wrong.` That is, all of the interesting features of the function fall within a hole in the integration grid when the range is large. The second issue is probably because the function is an even function and it doesn't check that you include identical infinite limits. – James Jan 06 '15 at 14:25
  • @James, post as answer? (There may also be a duplicate somewhere -- I'll see if I can find it) – Ben Bolker Jan 06 '15 at 14:35
  • not quite identical but closely related: http://stackoverflow.com/questions/25585316/r-integrate-returns-wrong-solution-is-using-wrong-quadrature-points ; http://r.789695.n4.nabble.com/puzzle-with-integrate-over-infinite-range-td2548525.html ; http://stackoverflow.com/questions/27297339/problems-with-integrate-returning-integral-of-0 – Ben Bolker Jan 06 '15 at 14:38
  • @BenBolker I think the problem is merely (at least the last part I found suprising), that `integrate` seems to change `-Inf` to `Inf` when used as upper limit. But I'm on thin ice here so please correct me if I'm wrong. – J.R. Jan 06 '15 at 15:03

1 Answers1

1

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

J.R.
  • 3,838
  • 1
  • 21
  • 25
  • 2
    you can see why this happens if you look at the code for `integrate`. It checks whether the upper bound is *finite* (`is.finite()`), not whether it is equal to `+Inf`. This seems like a bug(let), might be worth raising for discussion on `r-devel@r-project.org` ... it does otherwise respect the mathematical definition (e.g. `integrate(function(x) x, lower=2,upper=1)` is -1.5) ... – Ben Bolker Jan 06 '15 at 15:23
  • You are of course right! Yep just looked at the function (should probably have done it at first), a quick fix to make it behave exactly as expected. – J.R. Jan 06 '15 at 15:38
  • yes, but you will need a lot of test cases to make sure you got all the possibilities right (e.g. all combinations of (-Inf,+Inf,-1e100,1e100,0,-1,1) at both ends) – Ben Bolker Jan 06 '15 at 15:55