In this post, I would like to discuss how to generate Gamma distributed random variables. Gamma random variate has a number of applications. One of the most important application is to generate Dirichlet distributed random vectors, which plays a key role in topic modeling and other Bayesian algorithms.

A good starting point is a book by Kroese et al. [1] where detailed discussion about how to generate a number of different random distributed variables. In the book (Section 4.2.6), they list the following methods for Gamma distribution:

- Marsaglia and Tsang’s method [2]
- Ahrens and Dieter’s method [3]
- Cheng and Feast’s method [4]
- Best’s method [5]

Here, we focus on Marsaglia and Tsang’s method, which is used in GSL Library and Matlab “gamrnd” command (you can check this by typing “open gamrnd”).

## Gamma Distribution

As there are at least two forms of Gamma distribution, we focus the following formalism of PDF:

\[

\begin{align}

f(x; \alpha; \beta) = \frac{\beta^{\alpha} x^{\alpha – 1} e^{-\beta x} }{\Gamma(\alpha)}, \,\,\,\, x \geq 0

\end{align}

\]where \( \alpha > 0 \) is called the shape parameter and \( \beta > 0 \) is called the scale parameter.

## Marsaglia and Tsang’s Method

The algorithm works as follows for \( \mathbf{X} \sim \mbox{Gamma}(\alpha, 1) \) for \( \alpha \geq 1 \):

- Set \( d = \alpha – 1/3 \) and \( c = 1/\sqrt{9d} \).
- Generate \( Z \sim \mbox{N}(0,1) \) and \( U \sim \mbox{U}(0,1)\) independently.
- If \( Z > -1/c \) and \( \log U < \frac{1}{2} Z^{2} + d - d V + d \times ln V \), where \( V = (1+cZ)^{3}\), return \( \mathbf{X} = d \times V\); otherwise, go back to Step 2.

Two notes:

- This method can be easily extended to the cases where \( 1 > \alpha > 0 \). We just generate \( \mathbf{X} \sim \mbox{Gamma}(\alpha + 1, \beta)\), then \( \mathbf{X}^{\prime} = \mathbf{X} \times U^{1/\alpha}\) where \(U \sim (0,1)\). Thus, \( \mathbf{X}^{\prime} \sim \mbox{Gamma}(\alpha, \beta)\). See details in Page 9 of [2].
- For \( \beta \neq 1 \), we firstly generate \( \mathbf{X} \sim \mbox{Gamma}(\alpha, 1)\), then \( \mathbf{X}/\beta \sim \mbox{Gamma}(\alpha, \beta)\).

## Implementations

Here, we discuss some implementations details. We start from the original proposed C code from [2]:

1
2
3
4
5
6
7
8
| float rgama(float a) {
float d,c,x,v,u;
d5a-1./3.; c51./sqrt(9.*d);
for(;;) do {x5RNOR; v51.1c*x;} while(v,50.);
v5v*v*v; u5UNI;
if( u,1.-.0331*(x*x)*(x*x) ) return (d*v);
if( log(u),0.5*x*x1d*(1.-v1log(v)) ) return (d*v); }
} |

float rgama(float a) {
float d,c,x,v,u;
d5a-1./3.; c51./sqrt(9.*d);
for(;;) do {x5RNOR; v51.1c*x;} while(v,50.);
v5v*v*v; u5UNI;
if( u,1.-.0331*(x*x)*(x*x) ) return (d*v);
if( log(u),0.5*x*x1d*(1.-v1log(v)) ) return (d*v); }
}

which is slightly different from the algorithm proposed above. Again, this does not handle \( 1 > \alpha > 0 \) and \( \beta \neq 1 \).

This is the Matlab sample code from [1]:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| function x=gamrand(alpha,lambda)
% Gamma(alpha,lambda) generator using Marsaglia and Tsang method
% Algorithm 4.33
if alpha>1
d=alpha-1/3; c=1/sqrt(9*d); flag=1;
while flag
Z=randn;
if Z>-1/c
V=(1+c*Z)^3; U=rand;
flag=log(U)>(0.5*Z^2+d-d*V+d*log(V));
end
end
x=d*V/lambda;
else
x=gamrand(alpha+1,lambda);
x=x*rand^(1/alpha);
end |

function x=gamrand(alpha,lambda)
% Gamma(alpha,lambda) generator using Marsaglia and Tsang method
% Algorithm 4.33
if alpha>1
d=alpha-1/3; c=1/sqrt(9*d); flag=1;
while flag
Z=randn;
if Z>-1/c
V=(1+c*Z)^3; U=rand;
flag=log(U)>(0.5*Z^2+d-d*V+d*log(V));
end
end
x=d*V/lambda;
else
x=gamrand(alpha+1,lambda);
x=x*rand^(1/alpha);
end

This code does everything.

This is the code snippet extracted from GSL library:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
| double gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
{
/* assume a > 0 */
if (a < 1)
{
double u = gsl_rng_uniform_pos (r);
return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
}
{
double x, v, u;
double d = a - 1.0 / 3.0;
double c = (1.0 / 3.0) / sqrt (d);
while (1)
{
do
{
x = gsl_ran_gaussian_ziggurat (r, 1.0);
v = 1.0 + c * x;
}
while (v <= 0);
v = v * v * v;
u = gsl_rng_uniform_pos (r);
if (u < 1 - 0.0331 * x * x * x * x)
break;
if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
break;
}
return b * d * v;
}
} |

double gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
{
/* assume a > 0 */
if (a < 1)
{
double u = gsl_rng_uniform_pos (r);
return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
}
{
double x, v, u;
double d = a - 1.0 / 3.0;
double c = (1.0 / 3.0) / sqrt (d);
while (1)
{
do
{
x = gsl_ran_gaussian_ziggurat (r, 1.0);
v = 1.0 + c * x;
}
while (v <= 0);
v = v * v * v;
u = gsl_rng_uniform_pos (r);
if (u < 1 - 0.0331 * x * x * x * x)
break;
if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
break;
}
return b * d * v;
}
}

Again, this does everything. Note this code is based on a slightly different PDF: \( {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b}\).

**Reference**

[1] D.P. Kroese, T. Taimre, Z.I. Botev. **Handbook of Monte Carlo Methods**. John Wiley & Sons, 2011.

[2] G. Marsaglia and W. Tsang. **A simple method for generating gamma variables**. ACM Transactions on Mathematical Software, 26(3):363-372, 2000.

[3] J. H. Ahrens and U. Dieter. **Computer methods for sampling from gamma, beta, Poisson, and binomial distributions**. Computing, 12(3):223-246, 1974.

[4] R. C. H. Cheng and G. M. Feast. **Some simple gamma variate generators**. Computing, 28(3):290-295, 1979.

[5] D. J. Best. **A note on gamma variate generators with shape parameters less than unity**. Computing, 30(2):185-188, 1983.