DPQ (version 0.5-3)

qnormR: 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"))

Value

a numeric vector like the input q.

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.

Author

Martin Maechler

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

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