Gamma distribution

In probability theory and statistics, the gamma distribution is a continuous probability distribution. For integer values of the parameter k it is also known as the Erlang distribution.

Probability density function
The probability density function of the gamma distribution can be expressed in terms of the gamma function:


 * $$ f(x;k,\theta) = x^{k-1} \frac{e^{-x/\theta}}{\theta^k \, \Gamma(k)}

\ \mathrm{for}\ x > 0 \,\!$$

where $$k > 0$$ is the shape parameter and $$\theta > 0$$ is the scale parameter of the gamma distribution. (NOTE: this parameterization is what is used in the infobox and the plots.)

Alternatively, the gamma distribution can be parameterized in terms of a shape parameter $$\alpha = k$$ and an inverse scale parameter $$\beta = 1/\theta$$, called a rate parameter:


 * $$ g(x;\alpha,\beta) = x^{\alpha-1} \frac{\beta^{\alpha} \, e^{-\beta\,x} }{\Gamma(\alpha)}  \ \mathrm{for}\ x > 0 \,\!$$

Both parameterizations are common because they are convenient to use in certain situations and fields.

Properties
The cumulative distribution function can be expressed in terms of the incomplete gamma function,


 * $$ F(x;k,\theta) = \int_0^x f(u;k,\theta)\,du

= \frac{\gamma(k, x/\theta)}{\Gamma(k)} \,\!$$

The information entropy is given by:


 * $$S=k\theta+(1-k)\ln(\theta)+\ln(\Gamma(k))+(1-k)\psi(k)\,$$

where $$\psi(k)$$ is the polygamma function.

If $$X_i \sim \mathrm{Gamma}(\alpha_i, \beta)$$ for $$i=1, 2, \cdots, N$$ and $$\bar{\alpha} = \sum_{k=1}^N \alpha_i$$ then



\left[ Y = \sum_{i=1}^N X_i \right] \sim \mathrm{Gamma} \left( \bar{\alpha}, \beta \right) $$

provided all $$X_i$$ are independent. The gamma distribution exhibits infinite divisibility.

If $$X \sim \operatorname{Gamma}(k, \theta)$$, then $$\frac X \theta \sim \operatorname{Gamma}(k, 1)$$. Or, more generally, for any $$t > 0$$ it holds that $$tX \sim \operatorname{Gamma} (k, t \theta)$$. That is the meaning of &theta; (or &beta;) being the scale parameter.

Parameter estimation
The likelihood function is


 * $$L=\prod_{i=1}^N f(x_i;k,\theta)$$

from which we calculate the log-likelihood function


 * $$\ell=(k-1)\sum_{i=1}^N\ln(x_i)-\sum x_i/\theta-Nk\ln(\theta)-N\ln\Gamma(k)$$

Finding the maximum with respect to $$\theta$$ by taking the derivative and setting it equal to zero yields the maximum likelihood estimate of the $$\theta$$ parameter:


 * $$\theta=\frac{1}{kN}\sum_{i=1}^N x_i$$

Substituting this into the log-likelihood function gives:


 * $$\ell=(k-1)\sum_{i=1}^N\ln(x_i)-Nk-Nk\ln\left(\frac{\sum x_i}{kN}\right)-N\ln\Gamma(k)$$

Finding the maximum with respect to $$k$$ by taking the derivative and setting it equal to zero yields:


 * $$\ln(k)-\psi(k)=\ln\left(\frac{1}{N}\sum_{i=1}^N x_i\right)-\frac{1}{N}\sum_{i=1}^N\ln(x_i)$$

where $$\psi(k)=\frac{\Gamma'(k)}{\Gamma(k)}$$ is the digamma function.

There is no closed-form solution for $$k$$. The function is numerically very well behaved, so if a numerical solution is desired, it can be found using Newton's method. An initial value of $$k$$ can be found either using the method of moments, or using the approximation:


 * $$\ln(k)-\psi(k) \approx \frac{1}{k}\left(\frac{1}{2} + \frac{1}{12k+2}\right)$$

If we let $$s = \ln\left(\frac{1}{N}\sum_{i=1}^N x_i\right)-\frac{1}{N}\sum_{i=1}^N\ln(x_i),$$ then $$k$$ is approximately


 * $$k \approx \frac{3-s+\sqrt{(s-3)^2 + 24s}}{12s}$$

which is within 1.5% of the correct value.

Generating Gamma random variables
Given the scaling property above, it is enough to generate Gamma variables with $$\beta = 1$$ as we can later convert to any value of &beta; with simple division.

Using the fact that if $$X \sim \operatorname{Gamma}(1, 1)$$, then also $$X \sim \operatorname {Exp} (1)$$, and the method of generating exponential variables, we conclude that if U is uniformly distributed on (0, 1 ], then $$-\ln U \sim \operatorname{Gamma} (1, 1)$$. Now, using the "&alpha;-addition" property of Gamma distribution, we expand this result:


 * $$\sum _{k=1} ^n {-\ln U_k} \sim \operatorname{Gamma} (n, 1)$$,

where $$U_k$$ are all uniformly distributed on (0, 1 ] and independent.

All that is left now is to generate a variable distributed as $$\operatorname{Gamma} (\delta, 1)$$ for $$0 < \delta < 1$$ and apply the "&alpha;-addition" property once more. This is the most difficult part, however.

We provide an algorithm without proof. It is an instance of the acceptance-rejection method:


 * 1) Let m be 1.
 * 2) Generate $$V_{2m - 1}$$ and $$V_{2m}$$ &mdash; independent uniformly distributed on (0, 1 ] variables.
 * 3) If $$V_{2m - 1} \le v_0$$, where $$v_0 = \frac e {e + \delta}$$, then go to step 4, else go to step 5.
 * 4) Let $$\xi_m = \left( \frac {V_{2m - 1}} {v_0} \right) ^{\frac 1 \delta}, \ \eta_m = V_{2m} \xi _m^ {\delta - 1}$$.  Go to step 6.
 * 5) Let $$\xi_m = 1 - \ln {\frac {V_{2m - 1} - v_0} {1 - v_0}}, \ \eta_m = V_{2m} e^{-\xi_m}$$.
 * 6) If $$\eta_m > \xi_m^{\delta - 1} e^{-\xi_m}$$, then increment m and go to step 2.
 * 7) Assume $$\xi = \xi_m$$ to be the realization of $$\operatorname {Gamma} (\delta, 1)$$.

Now, to summarize,


 * $$\frac 1 \beta \left( \xi - \sum _{k=1} ^{[\alpha]} {\ln U_k} \right) \sim \operatorname{Gamma}(\alpha, \beta)$$ ,

where $$[\alpha]$$ is the integral part of &alpha;, &xi; has been generating using the algorithm above with $$\delta = \{\alpha\}$$ (the fractional part of &alpha;), $$U_k$$ and $$V_l$$ are distributed as explained above and are all independent.

Related distributions

 * $$X \sim \mathrm{Exponential}(\theta)$$ is an exponential distribution if $$X \sim \mathrm{Gamma}(1, \theta)$$.
 * $$cX \sim \mathrm{Gamma}(k, c\theta)$$ if $$X \sim \mathrm{Gamma}(k, \theta)$$ for any c > 0.
 * $$Y \sim \mathrm{Gamma}(N, \theta)$$ is a gamma distribution if $$Y = X_1 + \cdots + X_N$$ and if the $$X_i \sim \mathrm{Exponential}(\theta)$$ are all independent and share the same parameter $$\theta$$.
 * $$X \sim \chi^2(\nu)$$ is a chi-square distribution if $$X \sim \mathrm{Gamma}(k=\nu/2, \theta = 2)$$.
 * If $$k$$ is an integer, the gamma distribution is an Erlang distribution (so named in honor of A. K. Erlang) and is the probability distribution of the waiting time until the $$k$$-th "arrival" in a one-dimensional Poisson process with intensity $$1/\theta$$.
 * $$X \sim \mathrm{Gamma}(k, \theta)$$ then $$Y \sim \mathrm{InvGamma}(k, \theta^{-1})$$ if $$Y = 1/X$$, where $$\mathrm{InvGamma}$$ is the inverse-gamma distribution.
 * $$Y = X_1/(X_1+X_2) \sim \mathrm{Beta}$$ is a beta distribution if $$X_1 \sim \mathrm{Gamma}$$< and $$X_2 \sim \mathrm{Gamma}$$ and are also independent.
 * $$Y \sim \mathrm{Maxwell}(\beta)$$ is a Maxwell-Boltzmann distribution if $$X \sim \mathrm{Gamma}(\alpha = 3/2, \beta)$$.
 * $$Y \sim N(\mu = \alpha \beta, \sigma^2 = \alpha \beta^2)$$ is a normal distribution as $$Y = \lim_{\alpha \to \infty} X$$ where $$X \sim \mathrm{Gamma}(\alpha, \beta)$$.
 * The real vector $$(X_1/S,\ldots,X_n/S)\sim \operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_n)$$ follows a Dirichlet distribution if $$X_i\sim\operatorname{Gamma}(\alpha_i,1)$$ are independent, and $$S=X_1+\cdots+X_n$$.