DPQ (version 0.6-2)

qgammaAppr: Compute (Approximate) Quantiles of the Gamma Distribution

Description

Compute approximations to the quantile (i.e., inverse cumulative) function of the Gamma distribution.

Usage

qgammaAppr(p, shape, lower.tail = TRUE, log.p = FALSE,
           tol = 5e-07)
qgamma.R  (p, alpha, scale = 1, lower.tail = TRUE, log.p = FALSE,
           EPS1 = 0.01, EPS2 = 5e-07, epsN = 1e-15, maxit = 1000,
           pMin = 1e-100, pMax = (1 - 1e-14),
           verbose = getOption("verbose"))
qgammaApprKG(p, shape, lower.tail = TRUE, log.p = FALSE)
 
qgammaApprSmallP(p, shape, lower.tail = TRUE, log.p = FALSE)
 
 
.qgammaApprBnd(a, logEps = -36.043653389117156)

Value

numeric

Arguments

p

numeric vector (possibly log tranformed) probabilities.

shape, alpha, a

shape parameter, non-negative.

scale

scale parameter, non-negative, see qgamma.

lower.tail, log.p

logical, see, e.g., qgamma(); must have length 1.

tol

tolerance of maximal approximation error.

EPS1

small positive number. ...

EPS2

smaller positive number. ...

epsN

even smaller positive number, typically circa 4.5\(\epsilon_c\).

maxit

maximal number of iterations for qgamma.R().

pMin, pMax

boundaries for p. in qgamma.R(p), where p. is defined to be in \([0, 1]\), i.e., not on the log scale even when log.p = TRUE.

verbose

logical indicating if the algorithm should produce “monitoring” information.

logEps

\( = \log(\epsilon_c)\) := log(.Machine$double.eps), see the approximation .qgammaApprBnd(), below.

Author

Martin Maechler

Details

qgammaApprSmallP(p, a) should be a good approximation in the following situation when both p and shape \(= \alpha =: a\) are small :

If we look at Abramowitz&Stegun \(gamma*(a,x) = x^-a * P(a,x)\) and its series \(g*(a,x) = 1/gamma(a) * (1/a - 1/(a+1) * x + ...)\),

then the first order approximation \(P(a,x) = x^a * g*(a,x) ~= x^a/gamma(a+1)\) and hence its inverse \(x = qgamma(p, a) ~= (p * gamma(a+1)) ^ (1/a)\) should be good as soon as \(1/a >> 1/(a+1) * x\)

<==> x << (a+1)/a = (1 + 1/a)

<==> x < eps *(a+1)/a

<==> log(x) < log(eps) + log( (a+1)/a ) = log(eps) + log((a+1)/a) ~ -36 - log(a) where log(x) ~= log(p * gamma(a+1)) / a = (log(p) + lgamma1p(a))/a

such that the above

<==> (log(p) + lgamma1p(a))/a < log(eps) + log((a+1)/a)

<==> log(p) + lgamma1p(a) < a*(-log(a)+ log(eps) + log1p(a))

<==> log(p) < a*(-log(a)+ log(eps) + log1p(a)) - lgamma1p(a) =: bnd(a)

Note that qgammaApprSmallP() indeed also builds on lgamma1p().

.qgammaApprBnd(a) provides this bound \(bnd(a)\); it is simply a*(logEps + log1p(a) - log(a)) - lgamma1p(a), where logEps is \(\log(\epsilon)\) = log(eps) where eps <- .Machine$double.eps, i.e. typically (always?) logEps\(= \log \epsilon = -52 * \log(2) = -36.04365\).

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.

See Also

qgamma for R's Gamma distribution functions.

Examples

  ## TODO :  Move some of the curve()s from ../tests/qgamma-ex.R !!