5

I was doing some integration into a loop using integrate and I come up with an error I can't understand neither get rid of. Here is a MWE I could extract:

u_min = 0.06911363
u_max = 1.011011 
m = 0.06990648
s = 0.001092265
integrate(f = function(v){pnorm(v, mean = m, sd = s, lower.tail =  FALSE)}, u_min, u_max)

this returns an error "the integrale is probably divergent" which is obviously false. I tried to modify the parameters a little bit and got this working for example:

u_min <- 0.07
u_max <- 1.1
m <- 0.0699
s <- 0.00109
integrate(f = function(v){pnorm(v, mean = m, sd = s, lower.tail =  FALSE)}, u_min, u_max)

I tried to have a look into the integrate function with debugbut it's a wrapper of a Ccode. Also I'm not a specialist of quadrature techniques. I saw this SO post but couldn't make anything from it.

thanks

Community
  • 1
  • 1
ClementWalter
  • 4,814
  • 1
  • 32
  • 54
  • From `?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`. Your function is basically always 0. – nicola Dec 09 '15 at 11:20
  • @nicola So what would one do to not get an error returned? How would I be able to work around this? – Plinth May 09 '16 at 03:23
  • @Plinth you have to narrow the interval around values where the function is not zero, see nicola's answer below – ClementWalter May 09 '16 at 08:26

2 Answers2

5

The default tolerance of .Machine$double.eps^0.25 (= 0.0001220703) needs to be lowered. Try, for example, this:

f <- function(v) pnorm(v, mean = m, sd = s, lower.tail =  FALSE)
integrate(f, u_min, u_max, rel.tol = 1e-15)

## 0.0009421867 with absolute error < 1.1e-17
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
2

I'd use this work-around:

integrate(f = function(v){pnorm(v, mean = m, sd = s, lower.tail =  FALSE)}, 
      max(u_min,m-10*s),min(u_max,m+10*s))$value  + (u_min-m+10*s)*(u_min<m+10*s)

What I've done:

  • pnorm with lower.tail=FALSE is basically zero when very far on the right from the mean. So there is no point in "stretching" the right limit of the integral. So, when u_max > m+10*s, you just integrate to m + 10*s. You can of course change the 10 factor to add precision;
  • on the other hand, on the left pnorm is basically always 1; so you can enhance the left limit and the missing part is just u_min - m+10*s. Same logic as above.
nicola
  • 24,005
  • 3
  • 35
  • 56