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

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

Arguments

probability \(p\), \(1-p\), or \(\log(p)\), \(\log(1-p)\), depending on lower.tail and log.p.

mean of the normal distribution.

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.

Value

a numeric vector like the input q.

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

Examples

Run this code

# NOT RUN {
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)))

## Showing why/where R's qnorm() was poor up to 2020: log.p=TRUE extreme tail
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 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) ,]
# }