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

In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

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,
       ...)

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.

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.

Author(s)

Martin Maechler

See Also

Standard R's qnorm.

Examples

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))))
if(doIncr) .mpfr_erange_set(value = myERng)

log2(abs(.mpfr_erange()))# 62 62 if the above increase worked
(lrgOK <- all(log2(abs(.mpfr_erange())) >= 62))
                                        # currently(2023-01) 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  233   bits  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)

Rmpfr documentation built on Aug. 8, 2023, 5:14 p.m.