qpoisR: Pure R Implementation of R's qpois() with Tuning Parameters

Description

A pure R implementation, including many tuning parameter arguments, of R's own Mathlib C code algorithm, but with more flexibility.

It is using Vectorize(qpoisR1, *) where the hidden qpoisR1 works for numbers (aka ‘scalar’, length one) arguments only, the same as the C code.

Usage

qpoisR(p, lambda, lower.tail = TRUE, log.p = FALSE,
       yLarge = 4096, # was hard wired to 1e5
       incF = 1/64,   # was hard wired to .001
       iShrink = 8,   # was hard wired to 100
       relTol = 1e-15,# was hard wired to 1e-15
       pfEps.n = 8,   # was hard wired to 64: "fuzz to ensure left continuity"
       pfEps.L = 2,   # was hard wired to 64:   "   "   ..
       fpf = 4, # *MUST* be >= 1 (did not exist previously)
       trace = 0)

Arguments

p, lambda, lower.tail, log.p

qpois() standard argument, see its help page.

yLarge

a positive number; in R up to 2021, was internally hardwired to yLarge = 1e5: Uses more careful search for $y \ge y_L$, where $y$ is the initial approximate result, derived from a Cornish-Fisher expansiion.

incF

a positive “increment factor” (originally hardwired to 0.001), used only when y >= yLarge; defines the initial increment in the search algorithm as incr <- floor(incF * y).

iShrink

a positive increment shrinking factor, used only when y >= yLarge to define the new increment from the old one as incr <- max(1, floor(incr/iShrink)) where the LHS was hardired original to (incr/100).

relTol

originally hard wired to 1e-15, defines the convergence tolerance for the search iterations when y >= yLarge; the iterations stop when (new) incr <= y * relTol.

pfEps.n, pfEps.L

positive factors defining “fuzz to ensure left continuity”, both originally hardwired to 64. originally, the fuzz adjustment was

p <- p * (1 - 64 *.Machine$double.eps)

Now, pfEps.L is used if(log.p) is true and pfEps.n is used otherwise ("n"ormal case), and the adjustments also depend on lower.tail, and also on fpf :

fpf

a number larger than 1, together with pfEps.n determines the fuzz-adjustment to p in the case (lower=tail=FALSE, log.p=FALSE): with e <- pfEps.n * .Machine$double.eps, the adjustment p <- p * (1 + e) is made iff 1 - p > fpf*e.

trace

logical (or integer) specifying if (and how much) output should be produced from the algorithm.

Value

a numeric vector like p recycled to the common lengths of p and lambda.

Details

The defaults and exact meaning of the algorithmic tuning arguments from yLarge to fpf were experimentally determined are subject to change.

Examples

Run this code

# NOT RUN {
x <- 10*(15:25)
Pp <- ppois(x, lambda = 100, lower.tail = FALSE)  # no cancellation
qPp <- qpois(Pp, lambda = 100, lower.tail=FALSE)
table(x == qPp) # all TRUE ?
## future: if(getRversion() >= "4.2") stopifnot(x == qPp) # R-devel
qpRp <- qpoisR(Pp, lambda = 100, lower.tail=FALSE)
all.equal(x, qpRp, tol = 0)
stopifnot(all.equal(x, qpRp, tol = 1e-15))
# }