0

I am doing numerical integration in R with the function integrate. I am integrating a function fun1 and have discovered some nonsmooth behavior. Changing the integral's lower bound slightly produces a big jump in the calculated value of the integral; compare aaa and ccc below. aaa is approximately correct, as can be verified using online numerical integral calculators. However, ccc is far from the truth.

For illustrative purposes, I split the integration interval into two and sum the two parts. The result is very different from integrating over a single interval. This should not be the case, as algebraically the two are equal.

fun1 <- function(x, d){
 min( 0, -(exp(x)-0.1*abs(exp(x)-d)-0.05*(exp(x)-d))^(-1) ) * dnorm(x=x,mean=0,sd=0.05)
}
aaa <- integrate(f=fun1, d=0.8932789, lower=-3.159041, upper=40       )$value; aaa
bbb <- integrate(f=fun1, d=0.8932789, lower=-3.159041, upper=-3.157379)$value; bbb
ccc <- integrate(f=fun1, d=0.8932789, lower=-3.157379, upper=40       )$value; ccc
aaa-(bbb+ccc) # This should be (but is not) close to zero because algebraically, aaa-(bbb+ccc)=0

Questions: Why is this happening? How do I fix this? I.e. how do I get a correct value for ccc?

Richard Hardy
  • 375
  • 6
  • 20
  • 2
    Plot your function in the desired interval. Next, read the `Note` section of `?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` – nicola Aug 27 '20 at 14:28
  • 1
    https://stackoverflow.com/questions/20665689/integrate-a-very-peaked-function-in-r ; https://stackoverflow.com/questions/62122505/integrate-gives-totally-wrong-number ; https://stackoverflow.com/questions/26023899/integrate-in-r-gives-terribly-wrong-answer?rq=1 – Ben Bolker Aug 27 '20 at 14:34
  • 2
    Function `cubature::hcubature` (see [my answer](https://stackoverflow.com/a/62122769/8245406) to the question in @BenBolker [second link](https://stackoverflow.com/questions/63617774/nonsmooth-behavior-of-numerical-integral-in-r-function-integrate#comment112497359_63617774)) gives `aaa-(bbb+ccc)` equal to `[1] 1.994734e-07`. – Rui Barradas Aug 27 '20 at 16:31
  • @RuiBarradas, thank you! Currently it looks like `hcubature` is doing better than `integrate`. I have tried out quite a bit of tuning tricks to make `integrate` work properly, but `hcubature` gets a more accurate answer with way less programming effort. The only drawback for now is that it is about 10x slower than `integrate` for my data. – Richard Hardy Aug 28 '20 at 13:20

0 Answers0