Function integration in R returns an error
I have a function defined as
tail <- function(x) {585.1961*x^-2.592484}
When I tried to get the integral from 5000 to Inf, it returns an error:
> integrate(tail, 5000, Inf)
Error in integrate(tail, 5000, Inf) : the integral is probably divergent
However, the integration from 1000 to Inf and 1000 to 5000 worked fine:
> integrate(tail, 1000, Inf)
0.006134318 with absolute error < 2.5e-05
> integrate(tail, 1000, 5000)
0.005661634 with absolute error < 4.9e-09
Isn't it integrate(tail, 5000, Inf)
just equal integrate(tail, 1000, Inf) - integrate(tail, 1000, 5000)
? Why did this lead to a diverging integral?
+3
source to share
1 answer
Your default ( .Machine$double.eps^0.25
) is too large, so we change it:
> tail <- function(x) {585.1961*x^-2.592484}
> integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^0.5 )
0.0004727982 with absolute error < 1.5e-09
and this is about the same as:
> integrate(tail, 1000, Inf)$val-integrate(tail, 1000, 5000)$val
[1] 0.0004726847
Of course, you can just set your tolerance .Machine$double.eps
, but this comes at the cost of time:
> library(microbenchmark)
> a<-function(x) for(i in 1:50) integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^x )
> microbenchmark(a(0.5),a(0.7),a(1))
Unit: milliseconds
expr min lq median uq max neval
a(0.5) 10.44027 10.97920 11.12981 11.40529 19.70019 100
a(0.7) 13.02904 13.69813 13.95942 14.89460 23.02422 100
a(1) 15.14433 15.96499 16.12595 16.38194 26.27847 100
eg. a time increase of around 50 pct.
+2
source to share