qnormRR Documentation

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

Description

Compute's R level implementations of R's qnorm() as implemented in C code (in R's ‘Rmathlib’).

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", "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; "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, it simply the result of Vectorize(qnormR1, ...).

Value

a numeric vector like the input q.

Author(s)

Martin Maechler

References

For the asymptotic approximations in later versions, see the reference(s) on qnormAsymp's help page.

See Also

qnorm

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
lp <- - 4^(1:30) # effect of  'trace = *' :
qp0 <- qnormR(lp, lower.tail=FALSE, log.p=TRUE, trace=1, version="4.0.x")
qp1 <- qnormR(lp, lower.tail=FALSE, log.p=TRUE, trace=1, version="2020")
qp2 <- qnormR(lp, lower.tail=FALSE, log.p=TRUE, trace=1, version="2022-08")

matplot(-lp, cbind(qp0, qp1, qp2), log="xy", type="l", lwd=2, axes=FALSE,
        main = "qnormR(lp, log.p=TRUE, lower=F, version = * )")
sfsmisc::eaxis(1, nintLog=15, sub=2); sfsmisc::eaxis(2)
lines(-lp, sqrt(-2*lp), col=4)
leg <- as.expression(c(paste("version=", ver.qn), quote(sqrt(-2*lp))))
legend("top", leg, bty='n', col=1:4, 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) ,]