Gamma distributions

1 All moments

The moment generating function of the Gamma distribution is \[\Ex[\exp(\rv{x} s)]=\left(1-{\frac {s}{\lambda }}\right)^{-\alpha }{\text{ for }}s<\lambda\] which gives us expressions for all moments fairly easily: \[\begin{aligned} \Ex[\rv{x}]&=\left.\frac{\dd}{\dd s}\Ex[\exp(\rv{x} t)]\right|_{s=0}\\ &=\left.\frac{\dd}{\dd s}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha }\right|_{s=0}\\ &=\left.\frac{\alpha}{\lambda}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha -1}\right|_{s=0}\\ &=\alpha/\lambda\\ \Ex[\rv{x}^2] &=\left.\frac{\dd}{\dd s}\frac{\alpha}{\lambda}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha -1}\right|_{s=0}\\ &=\left.\frac{\alpha^2+\alpha}{\lambda^2}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha -2}\right|_{s=0}\\ &=\frac{\alpha(\alpha+1)}{\lambda^2}\\ \var[\rv{x}] &=\frac{\alpha}{\lambda^2} \quad\text{(handy later)}\\ \var[\rv{x}^2] &=\frac{\Gamma (\alpha) \Gamma(\alpha+4)-\Gamma^2 (\alpha+2)}{\Gamma^2 (\alpha)\lambda ^4 }\quad\text{(handy later)}\\ \Ex[\rv{x}^3] &=\left.\frac{\dd}{\dd s}\frac{\alpha^2+\alpha}{\lambda^2}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha -2}\right|_{s=0}\\ &=\frac{\alpha(\alpha+1)(\alpha+2)}{\lambda^3}\left(1-{\frac {s}{\lambda }}\right)^{-\alpha -2}\\ &=\frac{\alpha(\alpha+1)(\alpha+2)}{\lambda^3}\\ \Ex[\rv{x}^4] &=\frac{\alpha(\alpha+1)(\alpha+2)(\alpha+3)}{\lambda^4}\\ &\dots\\ \Ex[\rv{x}^n] &=\frac{\langle \alpha \rangle_{n}}{\lambda^n}\\ \end{aligned} \]

Here \(\langle \alpha \rangle_{n}:=\frac{\Gamma(\alpha+n)}{\Gamma(\alpha)}\) is the rising factorial.

2 As exponential family

The gamma distribution is a two-parameter exponential family with natural parameters \([\theta_1, \theta_2]=[\alpha − 1,−\lambda]^\top\), and natural statistics \([x,\ln(x)]^\top.\) In particular, it is a natural exponential family with quadratic variance function (Morris 1983).

The log-partition function \(A(\theta)\) is given by: \[ A(\theta)=-\left(\theta_1+1\right) \log \left(-\theta_2\right)+\log \Gamma\left(\theta_1+1\right) \] So in canonical form the density is \(f\left(x ; \theta_1, \theta_2\right)=\frac{1}{x} \exp \left(\theta_1 \log (x)+\theta_2 x-A(\theta)\right).\)

In terms of the natural parameters, the mean and variance are \[ \begin{aligned} \Ex[\rv{x}]&=\frac{\theta_1+1}{-\theta_2}\\ \var[\rv{x}]&=\frac{\theta_1+1}{\theta_2^2}. \end{aligned} \]

3 Linear combinations of Gammas

Is the Gamma family closed under addition? For fixed scale/rate parameters, yes. See Gamma-Dirichlet algebra.

If we are summing gamma variates which differ in the rate parameter \(\lambda\), the result is not in general Gamma distributed. It can be compactly expressed in terms of Gamma densities, for what that is worth. If we wish to know the distribution of the sum of a set of arbitrary Gamma random variables with different rate parameters, we can try to approximate the sum with moment matching (Mathai 1982; Moschopoulos 1985) if that feels worthwhile.

Note that multiplying a gamma RV by a scalar changes the rate, so gamma variates are not closed under affine combination as Gaussians ones are. The moral is that we cannot assume the convenient algebra of linear combination additivity as in the Gaussian process community. What kind of algebra do we get?

4 Gamma-Dirichlet algebra

There are various operations which give us similar conveniences, however. For those, we need to also be aware of the Dirichlet and Beta distributions. Here are some useful properties, drawn from, or extrapolated from, Dufresne (1998), Lin (2016), and Pérez-Abreu and Stelzer (2014) which use those properties.

First, we fix some notation. From here on, all variables denoted \(\rv{x}_{a}\) (for some \(a>0\), with or without superscript) have a \(\operatorname{Gamma}(a, 1)\) distribution; all variables denoted \(\rv{b}_{a, b}\) (for some \(a, b>0\), with or without superscript) have a \(\operatorname{Beta}_{a, b}\) distribution. In all expressions, the variables \(\rv{x}_{a_{1}}, \rv{x}_{a_{2}}', \ldots, \rv{b}_{a_{1}, b_{1}}, \rv{b}_{a_{2}, b_{2}}', \ldots\) are independent unless I say otherwise.

5 Conjugate prior for

Fink (1997) summarises Miller (1980), which extends Damsleth (1975):

Suppose that data \(x_1, \ldots, x_n\) are independent and identically distributed from a gamma process where both the shape, \(\alpha\), and the reciprocal scale, \(\beta\), parameters are unknown. The likelihood function, \(L\left(\alpha, \beta \mid x_1, \ldots x_n\right)\), proportional to the parameters, is \[ \begin{aligned} & L\left(\alpha, \beta \mid x_1, \ldots, x_n\right) \propto\left\{\begin{array}{cc} \frac{P^{\alpha-1} \exp (-\beta S)}{\left(\Gamma(\alpha) \beta^{-\alpha}\right)^n} & \text { where } \quad x_i>0, i=1 \ldots n \\ 0 & \text { otherwise } \end{array}\right. \\ & \text { where } \quad S=\sum_{i=1}^n x_i \quad P=\prod_{i=1}^n x_i . \end{aligned} \]

The sufficient statistics are \(n\), the number of data points, \(P\), the product of the data, and \(S\), the sum of the data. The factors of [the equation] proportional to parameters \(\alpha\) and \(\beta\) make up the kernel of the conjugate prior, \(\pi(\cdot, \cdot)\). We specify the conjugate prior with hyperparameters \(p, q, r, s>0\) \[ \pi(\alpha, \beta \mid p, q, r, s)=\left\{\begin{array}{cc} \frac{1}{K} \frac{p^{\alpha-1} \exp (-\beta q)}{\Gamma(\alpha)^r \beta^{-\alpha s}} & \text { where } \quad \alpha, \beta>0 \\ 0 & \text { otherwise. } \end{array}\right. \]

[…] The posterior joint distribution of \(\alpha\) and \(\beta\) is specified by the hyperparameters \[ p^{\prime}=p P \quad q^{\prime}=q+S \quad r^{\prime}=r+n \quad s^{\prime}=s+n . \]

From this we find the predictive, \[ \begin{aligned} &L\left(x|\alpha, \beta \right) \pi(\alpha, \beta \mid p, q, r, s) \mathrm{d}\alpha \mathrm{d}\beta \\ &\propto \iint \frac{x^{\alpha-1} \exp (-\beta x)}{\Gamma(\alpha) \beta^{-\alpha}} \frac{p^{\alpha-1} \exp (-\beta q)}{\Gamma(\alpha)^r \beta^{-\alpha s}}\mathrm{d}\alpha \mathrm{d}\beta \\ % &\propto \iint \frac{P^{\alpha-1} p^{\alpha-1} \exp -\beta S \exp (-\beta q)}{\Gamma(\alpha)^r \beta^{-\alpha s}\left(\Gamma(\alpha) \beta^{-\alpha}\right)^n}\mathrm{d}\alpha \mathrm{d}\beta\\ &\propto \iint \frac{(x p)^{\alpha-1} \exp -\beta (x+q)}{\Gamma(\alpha)^{r+1} \beta^{-\alpha (x+n)}}\mathrm{d}\alpha \mathrm{d}\beta \end{aligned} \] The predictive is reasonably intractable e.g. because there is nothing obviously nice to do with the \(1/\Gamma(\alpha)\) term; even estimating the mean requires approximations.

Maybe if you really want to have a good conjugate prior relation for a non-negative variate we should consider something a little (but only a little) less messy, such as [inverse Gaussian] (./inverse_gaussian_distribution.qmd) or lognormal distributions.

6 Generalized Gamma Convolution

As noted under divisible distributions, the class of Generalized Gamma Convolution (GGC) is a construction that represents some startling (to me) processes as a certain type of generalization of Gamma distributions. This family includes Pareto (Thorin 1977b) and Lognormal (Thorin 1977b) distributions. Those Thorin papers introduced the idea originally; possibly it is easy to start from one of the textbooks or overviews (Bondesson 2012; James, Roynette, and Yor 2008; Steutel and van Harn 2003; Barndorff-Nielsen, Maejima, and Sato 2006).

AFAICT this allows us to prove lots of nice things about such distributions. It is less easy to get implementable computational methods this way.

The GGC convolves a Gamma distribution with some measure and makes a new divisible distribution. James, Roynette, and Yor (2008):

we say that a positive r.v. \(\Gamma\) is a generalized gamma convolution \((\mathrm{GGC})\) - … if there exists a positive Radon measure \(\mu\) on \(] 0, \infty[\) such that: \[ \begin{aligned} E\left[e^{-\lambda \Gamma}\right]=& \exp \left\{-\int_{0}^{\infty}\left(1-e^{-\lambda x}\right) \frac{d x}{x} \int_{0}^{\infty} e^{-x z} \mu(d z)\right\} \\ =& \exp \left\{-\int_{0}^{\infty} \log \left(1+\frac{\lambda}{z}\right) \mu(d z)\right\} \\ \text { with: } & \int_{] 0,1]}|\log x| \mu(d x)<\infty \text { and } \int_{[1, \infty[} \frac{\mu(d x)}{x}<\infty. \end{aligned} \]

Barndorff-Nielsen, Maejima, and Sato (2006) and Pérez-Abreu and Stelzer (2014) generalize the GGC to vector- and matrix-valued distributions.

7 Parameter estimation

The method of moments is obvious. The Maximum likelihood version is surprisingly fiddly and has no closed form, but a low-bias closed-form approximation is given by Ye and Chen (2017). Wikipedia’s summary:

The estimate for the shape \(k\) is

\[ \hat{k} = \frac{N \sum_{i=1}^N x_i}{N \sum_{i=1}^N x_i \ln x_i - \sum_{i=1}^N x_i \sum_{i=1}^N \ln x_i} \]

and the estimate for the scale \(\theta\) is

\[ \hat{\theta} = \frac{1}{N^2} \left(N \sum_{i=1}^N x_i \ln x_i - \sum_{i=1}^N x_i \sum_{i=1}^N \ln x_i\right) \]

Using the sample mean of \(x\), the sample mean of \(\ln x\), and the sample mean of the product \(x\ln x\) simplifies the expressions to:

\[ \hat{k} = \bar{x} / \hat{\theta} \] \[ \hat{\theta} = \overline{x\ln{x}} - \bar{x} \overline{\ln{x}}. \]

If the rate parameterization is used, the estimate of \(\hat{\beta} = 1/\hat{\theta}\).

These estimators are not strictly maximum likelihood estimators, but are instead referred to as mixed type log-moment estimators. They have however similar efficiency as the maximum likelihood estimators.

Although these estimators are consistent, they have a small bias. A bias-corrected variant of the estimator for the scale \(\theta\) is

\[ \tilde{\theta} = \frac{N}{N - 1} \hat{\theta} \]

A bias correction for the shape parameter \(k\) is given as (Louzada, Ramos, and Ramos 2019)

\[ \tilde{k} = \hat{k} - \frac{1}{N} \left(3 \hat{k} - \frac{2}{3} \left(\frac{\hat{k}}{1 + \hat{k}}\right) - \frac{4}{5} \frac{\hat{k}}{(1 + \hat{k})^2} \right) \]

A Bayesian update is given with a similar form in Louzada and Ramos (2018) but they have a slightly different parameterisation, so I will need to come back to that when I have time to translate that.

8 Simulating Gamma variates

9 References