ggamma (version 1.0.2)

G.Gamma: Generalized Gamma Probability Distribution

Description

Fast implementation of density, distribution function, quantile function and random generation for the Generalized Gamma probability distribution.

Usage

dggamma(x, a, b, k, log = F)
pggamma(q, a, b, k, lower.tail = TRUE, log.p = FALSE)
qggamma(p, a, b, k, lower.tail = TRUE, log.p = FALSE)
rggamma(n, a, b, k)

Arguments

x, q

vector of quantiles.

a, b, k

Parameters of the distribution, all of which must be positive.

log, log.p

logical; if TRUE, probabilities p are given as log(p).

lower.tail

logical; if TRUE (default), probabilities are \(P[X \le x]\) otherwise, \(P[X > x]\).

p

vector of probabilities.

n

number of observations. If length(n) > 1, the length is taken to be the number required.

Details

The generalized gamma distribution proposed by Stacy (1962) has parameters \(a, d, p\), but here we adopt the reparametrization $$ a = a $$ $$ b = p $$ $$ k = \frac{d}{p} $$ as is used by the R package *flexsurv*.

Probability density function $$ f(x) = \frac{b x^{bk-1} \exp[-(x/a)^b]}{a^{bk} \Gamma(k)} $$

Cumulative density function $$ F(x) = \frac{\gamma(k, (x/a)^b)}{\Gamma(k)} $$

The above function can be written in terms of a \(Gamma(\alpha, \beta)\). Let \(T \sim Gamma(k, 1)\) and its cumulative distribution be denoted as \(F_T(t)\), then the cumulative density function of the generalized gamma distribution can be written as $$ F(x) = F_T( (x/a)^b ) $$ which allows us to write the quantile function of the generalized gamma in terms of the gamma one (\(Q_T(u)\) is the quantile function of \(T\)) $$ Q(u) = (Q_T(u) \cdot a)^{1/b} $$ from which random numbers can be drawn.

References

Stacy, E. W. (1962). A generalization of the gamma distribution. The Annals of mathematical statistics, 33(3), 1187-1192.

Examples

x = seq(0.001, 5, length=1000);
plot(x, dggamma(x, 3, 1.8, 0.5), col=2, type="l", lwd=4, ylim=c(0, 1));
lines(x, pggamma(x, 3, 1.8, 0.5), col=4, type="l", lwd=4, ylim=c(0, 1));
legend("right", c("PDF", "CDF"), col=c(2, 4), lwd=4);

r = rgamma(n = 100, 2, 2);
lik = function(params) -sum(dggamma(r, params[1], params[2], params[3], log=TRUE));
optPar = optim(lik, par=c(1, 1, 1), method="L-BFGS", lower=0.00001, upper=Inf)$par;
x = seq(0.001, 5, length=1000);
plot(x, dgamma(x, 2, 2), type="l", col=2, lwd=4, ylim=c(0, 1));
lines(x, dggamma(x, optPar[1], optPar[2], optPar[3]), col=4, lwd=4);
legend("topright", c("Gamma(shape=2, rate=2)", "MLE Gen. Gamma"), col=c(2, 4), lwd=4);