TAOCP 3.4.1 Exercise 16
Let f(x)=\frac{x^{a-1}e^{-x}}{\Gamma(a)}, \qquad x>0,\qquad 0<a\le1, be the gamma density of order $a$.
Section 3.4.1: Numerical Distributions
Exercise 16. ▶ [HM25] (J. H. Ahrens.) Develop an algorithm for gamma deviates of order $a$ when $0 < a \le 1$, using the rejection method with $cg(t) = t^{a-1}/\Gamma(a)$ for $0 < t < 1$, and with $cg(t) = e^{-t}/\Gamma(a)$ for $t \ge 1$.
Verified: yes
Solve time: 4m35s
Solution
Let
$$ f(x)=\frac{x^{a-1}e^{-x}}{\Gamma(a)}, \qquad x>0,\qquad 0<a\le1, $$
be the gamma density of order $a$.
The exercise prescribes the comparison function
$$ cg(x)= \begin{cases} \dfrac{x^{a-1}}{\Gamma(a)}, & 0<x<1,\[1ex] \dfrac{e^{-x}}{\Gamma(a)}, & x\ge1. \end{cases} $$
Since
$$ f(x)\le cg(x) $$
for all $x>0$, rejection sampling applies.
We first compute the normalizing constant $c$:
$$ \begin{aligned} c &=\int_0^\infty cg(x),dx \ &=\frac1{\Gamma(a)} \left( \int_0^1 x^{a-1},dx + \int_1^\infty e^{-x},dx \right). \end{aligned} $$
Now
$$ \int_0^1 x^{a-1},dx=\frac1a, \qquad \int_1^\infty e^{-x},dx=e^{-1}, $$
hence
$$ c=\frac{1/a+e^{-1}}{\Gamma(a)}. $$
Therefore the normalized proposal density is
$$ g(x)=\frac{cg(x)}{c}. $$
The total mass of the two branches is
$$ p_L = \frac{\int_0^1 cg(x),dx}{c} = \frac{1/a}{1/a+e^{-1}} = \frac1{1+a e^{-1}}, $$
and
$$ p_R = \frac{\int_1^\infty cg(x),dx}{c} = \frac{e^{-1}}{1/a+e^{-1}} = \frac{a e^{-1}}{1+a e^{-1}}. $$
We now derive the conditional proposal densities.
Left branch
Conditioned on $0<x<1$,
$$ g_L(x) = \frac{cg(x)}{\int_0^1 cg(t),dt} = \frac{x^{a-1}/\Gamma(a)}{(1/a)/\Gamma(a)} = a x^{a-1}, \qquad 0<x<1. $$
Its distribution function is
$$ G_L(x)=\int_0^x a t^{a-1},dt=x^a. $$
Hence, if $U\sim\mathrm{Uniform}(0,1)$,
$$ X=U^{1/a} $$
has density $g_L$.
The rejection probability is determined from
$$ \frac{f(x)}{cg(x)} = \frac{x^{a-1}e^{-x}/\Gamma(a)} {x^{a-1}/\Gamma(a)} = e^{-x}, \qquad 0<x<1. $$
Thus the left-branch acceptance probability is $e^{-x}$.
Right branch
Conditioned on $x\ge1$,
$$ g_R(x) = \frac{cg(x)}{\int_1^\infty cg(t),dt} = \frac{e^{-x}/\Gamma(a)}{e^{-1}/\Gamma(a)} = e^{1-x}, \qquad x\ge1. $$
Its distribution function is
$$ G_R(x) = \int_1^x e^{1-t},dt = 1-e^{1-x}. $$
Hence, if $U\sim\mathrm{Uniform}(0,1)$,
$$ X=1-\ln(1-U) $$
has density $g_R$.
The rejection probability is
$$ \frac{f(x)}{cg(x)} = \frac{x^{a-1}e^{-x}/\Gamma(a)} {e^{-x}/\Gamma(a)} = x^{a-1}, \qquad x\ge1. $$
Since $0<a\le1$, we have $a-1\le0$, therefore $x^{a-1}\le1$ for $x\ge1$, as required.
Algorithm G
This yields the following rejection algorithm for gamma deviates of order $a$, $0<a\le1$.
G1. [Choose branch.]
Generate $U_0\sim\mathrm{Uniform}(0,1)$.
If
$$ U_0\le p_L=\frac1{1+a e^{-1}}, $$
go to step G2. Otherwise go to step G3.
G2. [Left proposal.]
Generate independent $U_1,U_2\sim\mathrm{Uniform}(0,1)$.
Set
$$ X\leftarrow U_1^{1/a}. $$
If
$$ U_2\le e^{-X}, $$
accept $X$. Otherwise return to G1.
G3. [Right proposal.]
Generate independent $U_1,U_2\sim\mathrm{Uniform}(0,1)$.
Set
$$ X\leftarrow 1-\ln(1-U_1). $$
If
$$ U_2\le X^{a-1}, $$
accept $X$. Otherwise return to G1.
The accepted values have density
$$ f(x)=\frac{x^{a-1}e^{-x}}{\Gamma(a)}. $$
The expected number of proposals per accepted value is
$$ c=\frac{1/a+e^{-1}}{\Gamma(a)}. $$
$\square$