qnormR: Pure R version of R's 'qnorm()' with Diagnostics and Tuning...

In DPQ: Density, Probability, Quantile ('DPQ') Computations

Pure R version of R's qnorm() with Diagnostics and Tuning Parameters

Description

Computes R level implementations of R's qnorm() as implemented in C code (in R's ‘Rmathlib’), historically and present.

Usage

qnormR1(p, mu = 0, sd = 1, lower.tail = TRUE, log.p = FALSE, trace = 0, version = )
qnormR (p, mu = 0, sd = 1, lower.tail = TRUE, log.p = FALSE, trace = 0,
        version = c("4.0.x", "1.0.x", "1.0_noN", "2020-10-17", "2022-08-04"))

Arguments

p	probability p, 1-p, or \log(p), \log(1-p), depending on lower.tail and log.p.
mu	mean of the normal distribution.
sd	standard deviation of the normal distribution.
lower.tail, log.p	logical, see, e.g., qnorm().
trace	logical or integer; if positive or TRUE, diagnostic output is printed to the console during the computations.
version	a character string specifying which version or variant is used. The current default, "4.0.x" is the one used in R versions up to 4.0.x. The two "1.0*" versions are as used up to R 1.0.1, based on Algorithm AS 111, improved by a branch for extreme tails by Wichura, and a final Newton step which is only sensible when log.p=FALSE. That final stepped is skipped for version = "1.0_noN", “noN” := “no Newton”. "2020-10-17" is the one committed to the R development sources on 2020-10-17, which prevents the worst for very large |p| when log.p=TRUE. "2022-08-04" uses very accurate asymptotic formulas found on that date and provides full double precision accuracy also for extreme tails.

Details

For qnormR1(p, ..), p must be of length one, whereas qnormR(p, m, s, ..) works vectorized in p, mu, and sd. In the DPQ package source, qnormR is simply the result of Vectorize(qnormR1, ...).

Value

a numeric vector like the input q.

Author(s)

Martin Maechler

References

For versions "1.0.x" and "1.0_noN":
Beasley, J.D. and Springer, S.G. (1977) Algorithm AS 111: The Percentage Points of the Normal Distribution. JRSS C (Appied Statistics) 26, 118–121; \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2307/2346889")}.

For the asymptotic approximations used in versions newer than "4.0.x", i.e., "2020-10-17" and later, see the reference(s) on qnormAsymp's help page.

See Also

qnorm, qnormAsymp.

Examples

qR <- curve(qnormR, n = 2^11)
abline(h=0, v=0:1, lty=3, col=adjustcolor(1, 1/2))
with(qR, all.equal(y, qnorm(x), tol=0)) # currently shows TRUE
with(qR, all.equal(pnorm(y), x, tol=0)) # currently: mean rel. diff.: 2e-16
stopifnot(with(qR, all.equal(pnorm(y), x, tol = 1e-14)))

(ver.qn <- eval(formals(qnormR)$version)) # the possible versions
(doExtras <- DPQ:::doExtras()) # TRUE e.g. if interactive()
lp <- - 4^(1:30) # effect of  'trace = *' :
qpAll <- sapply(ver.qn, function (V)
    qnormR(lp, log.p=TRUE, trace=doExtras, version = V))
head(qpAll) # the "1.0" versions underflow quickly ..

cAdj <- adjustcolor(palette(), 1/2)
matplot(-lp, -qpAll, log="xy", type="l", lwd=3, col=cAdj, axes=FALSE,
        main = "- qnormR(lp, log.p=TRUE, version = * )")
sfsmisc::eaxis(1, nintLog=15, sub=2); sfsmisc::eaxis(2)
lines(-lp, sqrt(-2*lp), col=cAdj[ncol(qpAll)+1])
leg <- as.expression(c(paste("version=", ver.qn), quote(sqrt(-2 %.% lp))))
matlines(-lp, -qpAll[,2:3], lwd=6, col=cAdj[2:3])
legend("top", leg, bty='n', col=cAdj, lty=1:3, lwd=2)

## Showing why/where R's qnorm() was poor up to 2020: log.p=TRUE extreme tail
##% MM: more TODO? --> ~/R/MM/NUMERICS/dpq-functions/qnorm-extreme-bad.R
qs <- 2^seq(0, 155, by=1/8)
lp <- pnorm(qs, lower.tail=FALSE, log.p=TRUE)
## The inverse of pnorm() fails BADLY for extreme tails:
## this is identical to qnorm(..) in R <= 4.0.x:
qp <- qnormR(lp, lower.tail=FALSE, log.p=TRUE, version="4.0.x")
## asymptotically correct approximation :
qpA <- sqrt(- 2* lp)
##^
col2 <- c("black", adjustcolor(2, 0.6))
col3 <- c(col2, adjustcolor(4, 0.6))
## instead of going toward infinity, it converges at  9.834030e+07 :
matplot(-lp, cbind(qs, qp, qpA), type="l", log="xy", lwd = c(1,1,3), col=col3,
        main = "Poorness of qnorm(lp, lower.tail=FALSE, log.p=TRUE)",
        ylab = "qnorm(lp, ..)", axes=FALSE)
sfsmisc::eaxis(1); sfsmisc::eaxis(2)
legend("top", c("truth", "qnorm(.) = qnormR(., \"4.0.x\")", "asymp. approx"),
       lwd=c(1,1,3), lty=1:3, col=col3, bty="n")

rM <- cbind(lp, qs, 1 - cbind(relE.qnorm=qp, relE.approx=qpA)/qs)
rM[ which(1:nrow(rM) %% 20 == 1) ,]

DPQ documentation built on Nov. 3, 2023, 5:07 p.m.