DPQ (version 0.5-3)

qnormAppr: Approximations to 'qnorm()', i.e., \(z_\alpha\)

Description

Approximations to the standard normal (aka “Gaussian”) quantiles, i.e., the inverse of the normal cumulative probability function.

The qnormUappr*() are relatively simple approximations from Abramowitz and Stegun, computed by Hastings(1955): qnormUappr() is the 4-coefficient approximation to (the upper tail) standard normal quantiles, qnorm(), used in some qbeta() computations.

qnormUappr6() is the “traditional” 6-coefficient approximation to qnorm(), see in ‘Details’.

Usage

qnormAppr(p)
qnormUappr(p, lp = .DT_Clog(p, lower.tail=lower.tail, log.p=log.p),
           lower.tail = FALSE, log.p = FALSE,
           tLarge = 1e10)
qnormUappr6(p,
            lp = .DT_Clog(p, lower.tail=lower.tail, log.p=log.p),
               # ~= log(1-p) -- independent of lower.tail, log.p
            lower.tail=FALSE, log.p=FALSE,
            tLarge = 1e10)

Value

numeric vector of (approximate) normal quantiles corresponding to probabilities p

Arguments

p

numeric vector of probabilities, possibly transformed, depending on log.p. Does not need to be specified, if lp is instead.

lp

log(1 - p*), assuming \(p*\) is the lower.tail=TRUE, log.p=FALSE version of p. If passed as argument, it can be much more accurate than when computed from p by default.

lower.tail

logical; if TRUE (not the default here!), probabilities are \(P[X \le x]\), otherwise (by default) upper tail probabilities, \(P[X > x]\).

log.p

logical; if TRUE, probabilities \(p\) are given as \(\log(p)\) in argument p.

tLarge

a large number \(t0\); if \(t >= t0\), where \(t := sqrt(-2 * lp)\), the result will be \(= t\).

Author

Martin Maechler

Details

qnormAppr(p) uses the simple 4 coefficient rational approximation to qnorm(p), provided by Abramowitz and Stegun (26.2.22), p.933, to be used only for \(p > 1/2\) and typically qbeta() computations, e.g., qbeta.R.
The relative error of this approximation is quite asymmetric: It is mainly < 0.

qnormUappr(p) uses the same rational approximation directly for the Upper tail where it is relatively good, and for the lower tail via “swapping the tails”, so it is good there as well.

qnormUappr6(p, *) uses the 6 coefficient rational approximation to qnorm(p, *), from Abramowitz and Stegun (26.2.23), again mostly useful in the outer tails.

References

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.

Hastings jr., Cecil (1955) Approximations for Digital Computers. Princeton Univ. Press.

See Also

qnorm (in base R package stats), and importantly, qnormR and qnormAsymp() in our DPQ.

Examples

pp <- c(.001, .005, .01, .05, (1:9)/10, .95, .99, .995, .999)
z_p <- qnorm(pp)
(R <- cbind(pp, z_p, qA = qnormAppr(pp),
            qUA = qnormUappr(pp, lower.tail= TRUE),
            qA6 = qnormUappr6(pp, lower.tail=TRUE)))
## Errors, absolute and relative:
cbind(pp, (relE <- cbind(
               errA  = z_p - R[,"qA" ],
               errUA = z_p - R[,"qUA"],
               rE.A  = 1 - R[,"qA" ]/z_p,
               rE.UA = 1 - R[,"qUA"]/z_p,
               rE.A6 = 1 - R[,"qA6"]/z_p)))

lp <- -c(1000, 500, 200, 100, 50, 20:10, seq(9.75, 0, by = -1/8))
qnormUappr(lp=lp) # 'p' need not be specified if 'lp' is

curve(qnorm(x, lower.tail=FALSE), n=1001)
curve(qnormUappr(x), add=TRUE,    n=1001, col = adjustcolor("red", 1/2))

## Error curve:
curve(qnormUappr(x) - qnorm(x, lower.tail=FALSE), n=1001,
      main = "Absolute Error of  qnormUappr(x)")
abline(h=0, v=1/2, lty=2, col="gray")

curve(qnormUappr(x) / qnorm(x, lower.tail=FALSE) - 1, n=1001,
      main = "Relative Error of  qnormUappr(x)")
 abline(h=0, v=1/2, lty=2, col="gray")

curve(qnormUappr(lp=x) / qnorm(x, log.p=TRUE) - 1, -200, -1, n=1001,
      main = "Relative Error of  qnormUappr(lp=x)"); mtext(" & qnormUappr6()  [log.p scale]", col=2)
curve(qnormUappr6(lp=x) / qnorm(x, log.p=TRUE) - 1, add=TRUE, col=2, n=1001)
abline(h=0, v=1/2, lty=2, col="gray")

curve(qnormUappr(lp=x) / qnorm(x, log.p=TRUE) - 1,
      -2000, -.1, ylim = c(-2e-4, 1e-4), n=1001,
      main = "Relative Error of  qnormUappr(lp=x)"); mtext(" & qnormUappr6()  [log.p scale]", col=2)
curve(qnormUappr6(lp=x) / qnorm(x, log.p=TRUE) - 1, add=TRUE, col=2, n=1001)
abline(h=0, v=1/2, lty=2, col="gray")

## zoom out much more --> see jump (to be closer to truth) at ~ x = -650'000:
curve(qnormUappr(lp=x) / qnorm(x, log.p=TRUE) - 1,
      -800000, -.1, ylim = 2e-5*c(-1,1), n=1001,
      main = "Relative Error of  qnormUappr(lp=x)  [log.p scale]")
abline(h=0, v=1/2, lty=2, col="gray")
(t_ <- pretty(sqrt(-2*c(-830000, 0))))
axis(3, at=-t_^2/2, labels=t_, mgp=c(1.3, .4,0))