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)
 
 

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.

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 Dec. 5, 2023, 3:05 a.m.