Numerical problems with qnorm

Question

I'm having a numeric issue using qnorm(psn()). The problem is numeric.

Firstly, the Skew-Normal CDF round the result, since psn(9) is not 1:

 library(sn)
 psn(9)
#[1] 1

then

 qnorm(psn(9))
#[1] Inf

And see that:

 qnorm(.9999999999999999)
#[1] 8.209536
 qnorm(.99999999999999999)
#[1] Inf

note that 8.209536 is not that big, so this rounding is very imprecise.

Then, my final problem is the calculation of qnorm(psn()), that is part of my Copula density. Any hint on how can I avoid these numerical problems?

Answer 1

(This is not a resolution to your dilemma, more of an explanation of why I think you're seeing this and perhaps not likely to find an easy solution.)

I think this is getting into the realm where normal floating-point precision isn't going to work for you. For instance, doing the inverse of your function:

options(digits=22)
pnorm(8.209536)
# [1] 0.99999999999999989
pnorm(8.209536) - 1
# [1] -1.1102230246251565e-16

which is very close to

.Machine$double.eps
# [1] 2.2204460492503131e-16

which, according to ?.Machine, is

double.eps: the smallest positive floating-point number 'x' such that
          '1 + x != 1'.  It equals 'double.base ^ ulp.digits' if either
          'double.base' is 2 or 'double.rounding' is 0; otherwise, it
          is '(double.base ^ double.ulp.digits) / 2'.  Normally
          '2.220446e-16'.

It might be possible to translate what you need into higher-precision using auxiliary packages like gmp or Rmpfr. (I don't know if they support qnorm-like operations.)