A gamma variate generator with shape parameter less than unity

The gamma distribution is one of the most important probability distributions in statistics. It has applications in a broad range of fields in natural sciences, engineering, and social sciences. The probability density is defined by

where $\alpha$ ($>0$) is the shape parameter, $\lambda$ ($>0$) the scale parameter, and $\Gamma(x)$ the gamma function.

Generating a gamma variate, a random number drawn from a gamma distribution, is an important and fundamental problem in statistical computing. Owing to its importance, various numerical algorithms have been developed since the early days of computer age [1, 3, 4, 8, 7, 9, 10, 6].

The nature of the gamma distribution is controlled by the shape parameter $\alpha$. Accordingly, we need to use different gamma variate generators. When $1<\alpha$, the distribution function increases from zero at $x=0$, has a maximum at $x=(\alpha-1)\lambda$, and then decays as $x$ further increases. To generate a gamma variate with $\alpha>1$, probably Marsaglia & Tsang [8]’s algorithm is one of the best, as detailed in literature [8, 10, 6]. When $\alpha=1$, the distribution is reduced to the exponential distribution, and it is easy to generate a variate.

When $0<\alpha<1$, which is the range of our interest, the gamma distribution monotonically decreases from positive infinity at $x=0$ to zero at $x\rightarrow\infty$. To generate a gamma variate for $0<\alpha<1$, various acceptance-rejection methods have been developed [1, 3, 4, 7]. Combining power-law and exponential distributions, Ahrens & Dieter [1] proposed a piecewise rejection method. Best [3] extended the Ahrens & Dieter [1] method, by adjusting a switching point and introducing a squeeze technique. Devroye [4] developed a different rejection method, based on the exponential power distribution. As of 2024, it is employed by a Python package, NumPy [2]. Using a generalized exponential distribution [5], Kundu & Gupta [7] developed a rejection method and its piecewise extensions. Finally, Tanizaki [9] developed a ratio-of-uniforms method, which works either for $0<\alpha<1$ and for $1<\alpha$.

In this note, we propose efficient gamma generators for $0<\alpha<1$, advancing an earlier study [7]. In Section 2, we develop an acceptance-rejection algorithm, using the generalized exponential distribution. In Section 3, we construct fractional functions for the squeeze technique. In Section 4, we discuss piecewise envelope functions. In Section 5, we evaluate the performance of the proposed and previous methods. Section 6 contains discussions and the summary.

We construct a base algorithm to generate a random variate drawn from the gamma distribution with shape $0<\alpha<1$. We focus on the $\lambda=1$ case, because we can obtain the results for $\lambda\neq 1$, by multiplying the outputs by $\lambda$.

In line of Kundu & Gupta [7], we employ the generalized exponential (GE) distribution [5],

where $\alpha$ ($>0$) and $\lambda$ ($>0$) are the shape and scale parameters. Note that Refs. [5] and [7] use different conventions for the scale parameter, and we employ Gupta & Kundu [5]’s. The cumulative distribution function (CDF) of the GE distribution is

Hereafter we consider the GE distribution with $\lambda=1$. We generate a uniform random variate $U_{1}$ in the $(0,1)$ range:

Equating $U_{1}$ with (4) with $\lambda=1$, we find that a random variate $x$ that follows the generalized exponential distribution $f_{\rm GE(\alpha,1)}$ is drawn by

To help our discussion, we use a parameter $\beta$, a function $g(x)$, and a rejection function $R_{1}(x)$.

The function $g(x)$ is continuous at $x=0$. When $x\geq 0$, as $g(0)=1$ and $g^{\prime}(x)=(1+x-e^{x})/(x^{2}e^{x})\leq 0$, we find $0<g(x)\leq 1$. Since $0<\beta<1$, we find

All the equal signs are met when $x=0$.

Comparing the GE (3) and gamma distributions (2), we rewrite (2) as follows.

Here, the relation $\alpha\Gamma(\alpha)=\Gamma(\alpha+1)$ is used. Eq. (10) suggests that we can apply a rejection method to the GE-distributed variate $x$, to obtain a gamma-distributed variate. Using another uniform variate $U_{2}\sim U(0,1)$, we evaluate the acceptance condition:

or an equivalent inequality:

If either (11) or (12) is satisfied, then we take this number $x$. If it is not satisfied, we discard the number, and then we go back to the beginning. The final output follows a gamma distribution with shape $\alpha$. These procedures are summarized in Algorithm 1 in Table 1. We need the total density $1/\Gamma(\alpha+1)$ to generate one gamma distribution.

We improve the speed of the algorithm by the so-called squeeze method. We avoid the power and exponential operations in the acceptance condition (11), by using the following inequalities for $x\geq 0$.

To show this, we start from several inequalities:

We can prove them in the following way.

For the left inequality of (13), we differentiate the rejection function

With help from the Taylor series $e^{x}=\sum_{n=0}^{\infty}\dfrac{x^{n}}{n!}$, we find

In the right-hand side, all the terms in the numerator and in the denominator are positive, but the numerator terms contain additional factor of $1/2,1/3,\cdots$. This indicates that the ratio of the numerator to the denominator is less than $1/2$. Together with the $x=0$ case, we find

This suggests, together with (14),

For the right inequality of (13), with help from (8), (14) and (15), we obtain

Thus (24) and (25) prove (13).

In practice, we utilize the inequalities (13) or equivalent logical expressions, before evaluating $R_{1}(x)$. Translating $\beta\rightarrow(1-\alpha)$ again, we show the final algorithm in Algorithm 2 (Table 2).

Next, we split the envelope function for the rejection method into two parts. Across the switching point at $x=s$ ($>0$), we consider the GE distribution in the left and an exponential tail in the right. With help from the CDF of the GE distribution (Eq. (4)), we rewrite the gamma distribution (Eq. (2)) as follows.

Here, $S_{L}$ and $S_{R}$ are the densities of the left and right parts,

and $R_{2}(x)$ is the second rejection function

which satisfies $0<R_{2}(x)<1$ for $x>s$. Equation (28) tells us that we can apply the rejection method to the GE distribution with $R_{1}(x)$ in the left part, and that we can apply the rejection method to the exponential density with $R_{2}(x)$ in the right part.

Using a uniform variate $U_{1}$, when $U_{1}\leq S_{L}/(S_{L}+S_{R})$, we proceed to the GE part. We generate the GE-distributed number in the range $[0,s]$. Using a uniform variate $U_{2}$, the GE variate $x$ can be drawn by

Then we evaluate the acceptance condition, as discussed in Section 3.

In the right part, we obtain an exponential variate from $U_{2}$,

and then we apply the rejection method. We can similarly apply the squeeze technique to (30). For $x\geq 1$, from (16) and (17), we obtain a squeeze relation

Finally, we choose $s$. First, we look for $s$ that minimizes the total density

The solution can be found by solving $(1-e^{-s})^{\alpha-1}+(\alpha-1){s^{\alpha-2}}-{s^{\alpha-1}}=0$ by a root finder. In this case, we find an approximation by trial and error,

This gives a near-minimum $S(\alpha,s^{*})$ within an error of $<1.5\times 10^{-6}$. Another choice is $s=1$, in analogy with Ahrens & Dieter [1]. This simplifies several parameters in the algorithm.

These procedures are shown in Algorithm 3 (Table 3). The first three lines initialize coefficients, and then the lines inside the loop generate the variate. We reuse $(S_{L}+S_{R})/S_{L}\cdot U_{1}$ or $(S_{L}+S_{R})/S_{R}\cdot(U_{1}-S_{L}/(S_{L}+S_{R}))$ as a uniform variate, in order to reduce the total number of random variates. When $s=1$, we can further simplify the code, which is a trivial task.

We have carried out several benchmark tests. We wrote C codes for the eight methods: 1) the Ahrens & Dieter [1] method, 2) the Best [3] method, 3) the Devroye [4] method, 4) the Kundu & Gupta [7] method (Algorithm 3 in their paper), 5) Algorithm 1 (Table 1), 6) Algorithm 2 (Table 2), 7) Algorithm 3 (Table 3) with $s=s^{*}$, and 8) Algorithm 3 with $s=1$. In the last case, we have simplified the algorithm, as mentioned in Sec. 4. We use a uniform random generator (gsl_rng_uniform) in the GNU Scientific Library. We use Intel oneAPI Compiler (icx) v2024.1 with the -lgsl -O0 -lm option and the clang compiler (clang) v14.0 with the -lgsl -lm option on AMD Ryzen 5955 processor on Ubuntu Linux 24.04.

Figure 1 shows the normalized number of random variates to generate a gamma distribution, as a function of $\alpha$. The curves indicate theoretical predictions. Algorithms 1 & 2 require $1/\Gamma(\alpha+1)$, the same number as the most efficient method to date, the Devroye [4] method (not shown). Algorithm 3’s number is given by (34). The circles indicate our Monte Carlo results for eleven shape parameters, $\alpha=0.01,0.1,0.2,0.3,\cdots 0.9,0.99$. In each case, $10^{8}$ random numbers are generated. They are in excellent agreement with the theories. The proposed methods, Algorithm 3 in particular, need fewer random numbers than most of the other algorithms. The $s=s^{*}$ case is better than the $s=1$ case by $\lesssim 0.8$%.

Next we compare the speed of the algorithms. Considering practical applications, two settings are considered [9]. In the first setting, we repeatedly draw one random variate $10^{8}$ times for all 11 parameters. We measure the elapsed time in seconds three times per method, and then take the average elapsed time. Figure 2(a) shows the results. The Devroye [4] method and Algorithm 2 are the fastest. They outperform the well established methods [1, 3]. Although the two compilers give different results, from Algorithms 1 and 2, we see that the squeeze technique is effective. Note that the Kundu & Gupta [7] method does not use the squeeze technique.

In the second setting, we draw $10^{8}$ random variates all at once, and then repeat this procedure for the 11 parameters. Some algorithms initialize internal parameters, before generating many random variates. For example, in Algorithm 3, the first three lines in Table 3 are used only once per parameter, before generating many variates at once. These lines are repeated each time we draw a variate in the first setting. Figure 2(b) shows the elapsed times. For reference, results with high-level compiler optimizations (the -O2 option) are also displayed. Now the piecewise algorithms such as the Best [3] method and Algorithm 3 become faster than in Fig. 2(a), suggesting that they are fast, once the parameters are initialized. Together with the Best [3] method and the Devroye [4] method, Algorithms 2–3 are close and promising. In Algorithm 3, the $s=1$ cases are always faster than the $s=s^{*}$ cases.

We have developed the rejection algorithms to generate a gamma variate with the shape parameter $0<\alpha<1$. In line of Kundu & Gupta [7], we have employed the GE distribution [5]. The GE distribution is useful, because it resembles the gamma distribution, and because its CDF has a simple form. Kundu & Gupta [7] employed the GE distribution with the scale parameter $2$, that is, $f_{\rm GE}(x;\alpha,2)$, while we employ the scale parameter $1$, $f_{\rm GE}(x;\alpha,1)$, to obtain better acceptance efficiency. Therefore, the parameters and the rejection function in their study are different from ours. Despite a different strategy, our acceptance efficiency is as good as the best one [4].

To improve the performance, we have constructed the squeeze functions in (13). As demonstrated, this makes the code faster. One may notice that the left-hand side (the lower bound) can be negative for $x>4/\beta$, but logically there is no problem. In addition, since $F_{\rm GE}(4/\beta;\alpha,1)\geq F_{\rm GE}(4;\alpha,1)=(1-e^{-4})^{\alpha}\gtrsim 0.982$, more than 98% of the random numbers are found for $x<4/\beta$, where the lower bound condition is useful.

We have considered piecewise extensions of our algorithms, with an exponential tail [1, 3, 7]. This makes the acceptance efficiency even better. The squeeze functions are presented in (33). The left inequality (the lower bound) is identical to one in the Best [3] method, while the right one (the upper bound) is a new addition. This addition actually makes the Best [3] method few percents faster (not shown). As for Algorithm 3, despite the near-optimum density $S(\alpha,s^{*})<S(\alpha,1)$, the $s=1$ case is always faster than the $s=s^{*}$ case, because $s=1$ simplifies the algorithm. Therefore $s=1$ is our choice. In the second setting (Fig. 2(b)), Algorithm 3 with $s=1$ is always faster than Algorithm 2, owing to its good acceptance rate. It is promising when drawing many variates at once. In the first setting (Fig. 2(a)), in agreement with Tanizaki [9], several methods become slower due to the parameter initialization. This is also the case for Algorithm 3. In contrast, Algorithm 2 remains fast in the first and second settings, because of its simplicity.

Considering several aspects, we conclude that Algorithm 2 is one of the best two generators, together with the Devroye [4] method. These two often outperform the others, including the well established methods [1, 3]. With Intel compiler, Algorithm 2 is the best, and the Devroye [4] method is the best with clang. The proposed methods are simple enough, and easy to implement, as shown in Tables 1–3. In addition to the compiler-dependence, there can be a CPU-dependence (and possibly a GPU-dependence). Considering various runtime environments, we want to have as many good algorithms as possible, and ours are certainly one of them.

In summary, we have developed gamma variate generators for $0<\alpha<1$, using the GE distribution, the squeeze technique, and the piecewise envelope functions. The proposed methods are excellent in acceptance efficiency. The numerical tests suggest that Algorithm 2 is one of the best two generators for $0<\alpha<1$. Algorithm 3 with $s=1$ can be a good option when drawing many variates at once. We hope that the proposed methods are useful in practical applications.

Program codes in C language are available from the author upon request.