G.Gamma: Generalized Gamma Probability Distribution
In ggamma: Generalized Gamma Probability Distribution

Description Usage Arguments Details References Examples

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

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)

`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 ≤ 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.

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} Γ(k)}

Cumulative density function

F(x) = \frac{γ(k, (x/a)^b)}{Γ(k)}

The above function can be written in terms of a Gamma(α, β). 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.

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

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);