Search code examples
rfloating-pointprecisionnormal-distribution

Adding floating point precision to qnorm/pnorm?

I would be interested to increase the floating point limit for when calculating qnorm/pnorm from their current level, for example:

x <- pnorm(10) # 1
qnorm(x) # Inf
qnorm(.9999999999999999444) # The highst limit I've found that still return a <<Inf number

Is that (under a reasonable amount of time) possible to do? If so, how?

Solution

  • If the argument is way in the upper tail, you should be able to get better precision by calculating 1-p. Like this:

    > x = pnorm(10, lower.tail=F)
> qnorm(x, lower.tail=F)
10

    I would expect (though I don't know for sure) that the pnorm() function is referring to a C or Fortran routine that is stuck on whatever floating point size the hardware supports. Probably better to rearrange your problem so the precision isn't needed.

    Then, if you're dealing with really really big z-values, you can use log.p=T:

    > qnorm(pnorm(100, low=F, log=T), low=F, log=T)
100

    Sorry this isn't exactly what you're looking for. But I think it will be more scalable -- pnorm hits 1 so rapidly at high z-values (it is e^(-x^2), after all) that even if you add more bits they will run out fast.