R: Approximations to 'qnorm()', i.e., z

qnormAppr

R Documentation

Approximations to 'qnorm()', i.e., z_α

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)

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 ≤ 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 t_0; if t ≥ t_0, where t := √{-2 \cdot lp}, the result will be = t.

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.

Value

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

Author(s)

Martin Maechler

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.

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