pqnormAsymp {DPQmpfr}

Asymptotic Approximations of Extreme Tail 'pnorm()' and 'qnorm()'

Description

These functions provide the first terms of asymptotic series approximations to pnorm()'s (extreme) tail, from Abramawitz and Stegun's 26.2.13 (p.932), or qnorm() where the approximations have been derived via iterative plugin using Abramowitz and Stegun's formula.

Usage

pnormAsymp(x, k, lower.tail = FALSE, log.p = FALSE)
qnormAsymp(p, lp = .DT_Clog(p, lower.tail = lower.tail, log.p = log.p),
           order, M_2PI =,
           lower.tail = TRUE, log.p = missing(p))

Arguments

x

positive (at least non-negative) numeric vector.

k

integer \ge 0 indicating how many terms the approximation should use; currently k \le 5.

p

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

lp

numeric (vector) of log(1-p) values; if not specified, computed from p, depending on lower.tail and log.p.

order

an integer in \{0,1,\dots,5\}, specifying the approximation order.

M_2PI

the number 2\pi in the same precision as p or lp, i.e., numeric or of class "mpfr".

lower.tail

logical; if true, probabilities are P[X \le x], otherwise upper tail probabilities, P[X > x].

log.p

logical; if TRUE (default for qnormAsymp !!), probabilities p are given as \log(p) in argument p or \log{(1-p)} in lp.

Details

see both help pages pnormAsymp and qnormAsymp from our package DPQ.

Value

vector/array/mpfr like first argument x or p or lp, respectively.

Author(s)

Martin Maechler

See Also

pnorm. The same functions “numeric-only” are in my DPQ package with more extensive documentation.

Examples

require("Rmpfr") # (in strong dependencies of this pkg {DPQmpfr})
x <- seq(1/64, 10, by=1/64)
xm  <- mpfr( x, 96)
"TODO"

## More extreme tails: ----------------------------------------------
##
## 1. pnormAsymp() ---------------------
lx <- c((2:10)*2, 25, (3:9)*10, (1:9)*100, (1:8)*1000, (2:7)*5000)
lxm <- mpfr(lx, 256)
Px <- pnorm(lxm, lower.tail = FALSE, log.p=TRUE)
PxA <- sapplyMpfr(setNames(0:5, paste("k =",0:5)),
                  pnormAsymp, x=lxm, lower.tail = FALSE, log.p=TRUE)
if(interactive())
  roundMpfr(PxA, 40)
# rel.errors :
relE <- asNumeric(1 - PxA/Px)
options(width = 99) -> oop # (nicely printing the matrices)
cbind(lx, relE)
matplot(lx, abs(relE), type="b", cex = 1/2, log="xy", pch=as.character(0:5),
        axes=FALSE,
        main = "|relE( <pnormAsymp(lx, k=*, lower.tail=FALSE, log.p=TRUE) )|")
sfsmisc::eaxis(1, sub10=2); sfsmisc::eaxis(2)
legend("bottom", paste("k =", 0:5), col=1:6, lty=1:5,
       pch = as.character(0:5), pt.cex=1/2, bty="n")
## NB: rel.Errors go down to  7e-59 ==> need precision of  -log2(7e-59) ~ 193.2 bits

## 2. qnormAsymp() ---------------------
QPx <- sapplyMpfr(setNames(0:5, paste("k =",0:5)),
                  function(k) qnormAsymp(Px, order=k, lower.tail = FALSE, log.p=TRUE))
(relE.q <- asNumeric(QPx/lx - 1))
         # note how consistent the signs are (!) <==> have upper/lower bounds

matplot(-asNumeric(Px), abs(relE.q), type="b", cex = 1/2, log="xy", pch=as.character(0:5),
        xlab = quote(-Px), axes=FALSE,
        main = "|relE( <qnormAsymp(Px, k=*, lower.tail=FALSE, log.p=TRUE) )|")
sfsmisc::eaxis(1, sub10=2); sfsmisc::eaxis(2)
legend("bottom", paste("k =", 0:5), col=1:6, lty=1:5,
       pch = as.character(0:5), pt.cex=1/2, bty="n")

options(oop) # {revert to previous state}