rnumerical-integration

'Strange' behavior of integrate function with truncdist in R

Mathematically the following is impossible for

library(truncdist)
q = function(x, L, R ) dtrunc(x, "exp", rate=0.1, a=L,b=R) 
integrate(q, L=2, R=3, lower  =0, upper = 27 )
integrate(q, L=2, R=3, lower  =0, upper = 29 )
integrate(q, L=2, R=3, lower  =27, upper = 29 )
integrate(q, L=2, R=3, lower  =0, upper = 30 )

We found the first integral giving a positive number, the second one evaluating to zero by adding the third interval which integrates itself to zero. Is this an issue in integrate or truncdist?

We can use the following to find more such issues

z=numeric()
for(i in 1:50){
  z[i]=integrate(q, L=2, R=3, lower  =0, upper = i)$value
}

What do I need to do to find the correct integrals (which all are 1 when integrating from 0 to i>=3)?

Solution

  • From help("integrate"):

    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.

    You found an example of this:

    curve(q(x, 2, 3), from = -1, to = 30)

    resulting plot

    You shouldn't integrate distribution density functions numerically. Use the cumulative distribution function:

    diff(ptrunc(c(0, 29), "exp", rate = 0.1, a = 2, b = 3))
#[1] 1