qnormI: Gaussian / Normal Quantiles `qnorm()` via Inversion

Description

Compute Gaussian or Normal Quantiles qnorm(p, *) via inversion of our “mpfr-ified” arbitrary accurate pnorm(), using our unirootR() root finder.

Usage


qnormI(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE,
       trace = 0, verbose = as.logical(trace),
       tol,
       useMpfr = any(prec > 53),
       give.full = FALSE,
       ...)

Value

If give.full is true, return a list, say r, of

unirootR(.) results, with length(r) == length(p).

Otherwise, return a “numeric vector” like p, e.g., of

class "mpfr" when p is.

Arguments

p: vector of probabilities.
mean: vector of means.
sd: vector of standard deviations.
log.p: logical; if TRUE, probabilities p are given as log(p).
lower.tail: logical; if TRUE (default), probabilities are \(P[X \le x]\) otherwise, \(P[X > x]\).

trace: integer passed to unirootR(). If positive, information about a search interval extension will be printed to the console.
verbose: logical indicating if progress details should be printed to the console.
tol: optionally the desired accuracy (convergence tolerance); if missing or not finite, it is computed as \(2^-{pr+2}\) where the precision \(pr\) is basically max(getPrec(p+mean+sd)).
useMpfr: logical indicating if mpfr arithmetic should be used.

give.full: logical indicating if the full result of unirootR() should be returned (when applicable).
...: optional further arguments passed to unirootR() such as maxiter, verbDigits, check.conv, warn.no.convergence, and epsC.

Author

Martin Maechler

Examples

Run this code

doX <- Rmpfr:::doExtras() # slow parts only if(doX)
cat("doExtras: ", doX, "\n")
p  <- (0:32)/32
lp <- -c(1000, 500, 200, 100, 50, 20:1, 2^-(1:8))
if(doX) {
  tol1 <- 2.3e-16
  tolM <- 1e-20
  tolRIlog <- 4e-14
} else { # use one more than a third of the points:
   ip <- c(TRUE,FALSE, rep_len(c(TRUE,FALSE,FALSE), length(p)-2L))
   p <-  p[ip]
  lp <- lp[ip]
  tol1 <- 1e-9
  tolM <- 1e-12
  tolRIlog <- 25*tolM
}

f.all.eq <- function(a,b)
  sub("^Mean relative difference:", '', format(all.equal(a, b, tol=0)))
for(logp  in c(FALSE,TRUE)) {
  pp <- if(logp) lp else p
  mp <- mpfr(pp, precBits = if(doX) 80 else 64) # precBits = 128 gave "the same" as 80
  for(l.tail in c(FALSE,TRUE)) {
      qn <- qnorm (pp, lower.tail = l.tail, log.p = logp)
     qnI <- qnormI(pp, lower.tail = l.tail, log.p = logp, tol = tol1)
     qnM <- qnormI(mp, lower.tail = l.tail, log.p = logp, tol = tolM)
     cat(sprintf("Accuracy of qnorm(*, lower.t=%-5s, log.p=%-5s): %s || qnI: %s\n",
                 l.tail, logp, f.all.eq(qnM, qn ),
                               f.all.eq(qnM, qnI)))
     stopifnot(exprs = {
        all.equal(qn,  qnI, tol = if(logp) tolRIlog else 4*tol1)
        all.equal(qnM, qnI, tol = tol1)
     })
  }
}

## useMpfr, using mpfr()  :
if(doX) {
  p2 <- 2^-c(1:27, 5*(6:20), 20*(6:15))
  e2 <- 88
} else {
  p2 <- 2^-c(1:2, 7, 77, 177, 307)
  e2 <- 60
}
system.time( pn2 <- pnorm(qnormI(mpfr(p2, e2))) ) # 4.1 or 0.68
           all.equal(p2, pn2, tol = 0) # 5.48e-29 // 5.2e-18
2^-e2
stopifnot(all.equal(p2, pn2, tol = 6 * 2^-e2)) # '4 *' needed


## Boundary -- from limits in mpfr maximal exponent range!
## 1) Use maximal ranges:
(old_eranges <- .mpfr_erange()) # typically -/+ 2^30
(myERng <- (1-2^-52) * .mpfr_erange(c("min.emin","max.emax")))
(doIncr <- !isTRUE(all.equal(unname(myERng), unname(old_eranges)))) # ==>
## TRUE only if long is 64-bit, i.e., *not* on Windows
if(doIncr) .mpfr_erange_set(value = myERng)

log2(abs(.mpfr_erange()))# 62 62 if(doIncr) i.e. not on Windows
(lrgOK <- all(log2(abs(.mpfr_erange())) >= 62)) # FALSE on Windows
## The largest quantile for which our mpfr-ized qnorm() does *NOT* underflow :
cM <- if(doX) { "2528468770.343293436810768159197281514373932815851856314908753969469064"
      } else    "2528468770.34329343681"
##               1 3 5 7 9  1 3 5 7 9 1 3 5 7 9 1 3 5 7 9 1 3 5 7 9 1 3 5 7 9 1 3 5 7 9 1 3
##                       10         20        30        40        50        60        70
(qM <- mpfr(cM))
(pM <- pnorm(-qM)) # precision   if(doX) 233 else 70  bits of precision ;
## |--> 0 on Windows {limited erange}; otherwise and if(doX) :
## 7.64890682545699845135633468495894619457903458325606933043966616334460003e-1388255822130839040
log(pM) # 233 bits: -3196577161300663205.8575919621115614148120323933633827052786873078552904

if(lrgOK) withAutoprint({
  
  try( qnormI(pM) ) ## Error: lower < upper not fulfilled (evt. TODO)
  ## but this works
  print(qnI <- qnormI(log(pM), log.p=TRUE)) #  -2528468770.343293436
  all.equal(-qM, qnI, tol = 0) # << show how close; seen  1.084202e-19
  stopifnot( all.equal(-qM, qnI, tol = 1e-18) )
})

if(FALSE) # this (*SLOW*) gives 21 x the *same* (wrong) result --- FIXME!
  qnormI(log(pM) * (2:22), log.p=TRUE)






if(doX) ## Show how bad it is (currently ca. 220 iterations, and then *wrong*)
 str(qnormI(round(log(pM)), log.p=TRUE, trace=1, give.full = TRUE))
if(requireNamespace("DPQ"))
  new("mpfr", as(DPQ::qnormR(pM, trace=1), "mpfr")) # as(*, "mpfr") also works for +/- Inf
  # qnormR1(p=         0, m=0, s=1, l.t.= 1, log= 0): q = -0.5
  #    somewhat close to 0 or 1: r := sqrt(-lp) =  1.7879e+09
  #    r > 5, using rational form R_3(t), for t=1.787897e+09  -- that is *not* accurate
  # [1] -94658744.369295865460462720............

## reset to previous status if needed
if(doIncr) .mpfr_erange_set( , old_eranges)