DPQ (version 0.5-9)

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)
 
 

Value

numeric

Arguments

p

numeric vector (possibly log tranformed) probabilities.

shape, alpha

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

small positive number. ...

epsN

small positive number. ...

maxit

maximal number of iterations. ...

pMin, pMax

boundaries for p. ...

verbose

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

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