qgammaAppr: Compute (Approximate) Quantiles of the Gamma Distribution

In DPQ: Density, Probability, Quantile ('DPQ') Computations

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)

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.

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.

Value

numeric

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.

See Also

qgamma for R's Gamma distribution functions.

Examples

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

DPQ documentation built on Aug. 21, 2026, 3 p.m.