Gamma Variates By Numerical Inversion

Why Numerically Invert a Gamma CDF?

The gamma distribution is used to represent positive quantities whose variability and skewness can be adjusted through its parameters. In simulation, a standard approach is to draw $U$ uniformly from $(0,1)$, then find the value whose gamma cumulative distribution function equals $U$. For many shape parameters, that inverse has no elementary expression, so it must be found numerically.

Numerical inversion is useful when the required gamma shape is not an integer and a special-purpose generator is unavailable. The essential task is a one-dimensional root search: find $x$ such that $F_X(x)=U$. Correctly distinguishing the gamma shape from its scale, and controlling the error in that search, are central to obtaining a valid simulated variate.

Definition: A gamma random variable with shape $\alpha>0$ and scale $\theta>0$, written $X\sim\operatorname{Gamma}(\alpha,\theta)$, has density $$ f_X(x)=\frac{x^{\alpha-1}e^{-x/\theta}}{\Gamma(\alpha)\theta^\alpha},\qquad x>0, $$ and cumulative distribution function $$ F_X(x)=\frac{1}{\Gamma(\alpha)}\int_0^{x/\theta} t^{\alpha-1}e^{-t}\,dt,\qquad x\ge 0. $$ Here $\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}e^{-t}\,dt$ is the gamma function. The scale parameter $\theta$ multiplies the natural unit of the distribution; it is not the rate parameter $1/\theta$.

Inverse Transform for a Gamma Variate

Let $U\sim\operatorname{Uniform}(0,1)$. Since the gamma CDF is continuous and strictly increasing for $x>0$, there is a unique positive solution $x$ to $F_X(x)=U$. Finding that solution produces a gamma variate.

Theorem: If $U\sim\operatorname{Uniform}(0,1)$ and $F_X$ is a continuous, strictly increasing CDF, then $X=F_X^{-1}(U)$ has CDF $F_X$. For a gamma distribution, numerical inversion approximates this inverse by solving $F_X(x)-U=0$ on $x>0$.

The result follows directly from the CDF. For any $x$ in the support, strict monotonicity means that $F_X^{-1}(U)\le x$ exactly when $U\le F_X(x)$. Therefore,

P(F_X^{-1}(U)\le x) =P(U\le F_X(x)) =F_X(x),

because a uniform random variable on $(0,1)$ has probability $F_X(x)$ of being at most $F_X(x)$. This argument establishes why inverse-CDF sampling works; the numerical method only approximates the inverse closely enough for the intended precision.

Evaluating the CDF and Finding Its Root

Write $z=x/\theta$. It is convenient to evaluate the unit-scale gamma CDF at $z$, then recover the desired variate using $x=\theta z$. For a noninteger shape, a numerical library typically evaluates the regularized incomplete gamma function. One series representation, useful for explaining the calculation when $z$ is moderate, is

P(\alpha,z) =\frac{1}{\Gamma(\alpha)}\int_0^z t^{\alpha-1}e^{-t}\,dt =\frac{e^{-z}z^\alpha}{\Gamma(\alpha+1)} \left[ 1+\frac{z}{\alpha+1} +\frac{z^2}{(\alpha+1)(\alpha+2)} +\cdots \right].

In this notation, $P(\alpha,z)$ denotes the regularized lower incomplete gamma function, so $F_X(x)=P(\alpha,x/\theta)$. In practice, use a reliable special-function routine rather than summing a long series by hand, particularly when the shape or argument is large.

A dependable root-finding method is bisection. Begin with a lower bound $L$ and upper bound $H$ satisfying $F_X(L)\le U\le F_X(H)$. Evaluate the CDF at the midpoint, retain the half-interval that still brackets the target, and repeat. The bracket width is halved at every iteration. A safeguarded Newton method can converge faster, using $f_X(x)$ as the derivative of $F_X(x)$, but it must preserve a valid bracket because an unrestricted Newton step can leave the support or overshoot.

Draw the uniform input.
Generate $U$ on $(0,1)$. Use an independent uniform draw for each gamma variate required.

Set the target equation.
For shape $\alpha$ and scale $\theta$, solve $P(\alpha,z)=U$ for $z>0$.

Find and refine a bracket.
Choose bounds whose CDF values lie on opposite sides of $U$, then apply bisection or a safeguarded root solver.

Convert to the requested scale.
Set $X=\theta z$. Check the CDF residual and the final bracket width against the required accuracy.

Worked Numerical Inversions

Worked Example: Shape 2 With a Numerical Root Search

A hypothetical simulation requires $X\sim\operatorname{Gamma}(2,3)$, where the second parameter is the scale. A uniform draw is $U=0.5000$. Numerically invert the CDF and report $X$ to three decimal places.

Method: bracketed inversion. For shape 2 and unit scale, the CDF is $P(2,z)=1-e^{-z}(1+z)$. This gives a convenient way to check the numerical root. Start with $L=1.5$ and $H=1.75$:

P(2,1.5)=1-e^{-1.5}(2.5)\approx0.4422, \qquad P(2,1.75)=1-e^{-1.75}(2.75)\approx0.5221.

Thus $0.4422<0.5000<0.5221$, so the target is bracketed. The first midpoint is $1.625$, where $P(2,1.625)\approx0.4831$. The root is therefore above $1.625$. The next midpoint of the retained interval $[1.625,1.75]$ is $1.6875$, where $P(2,1.6875)\approx0.5029$; the root lies below $1.6875$. Further bisection gives $z\approx1.678347$.

Convert from unit scale to scale 3:

X=\theta z=3(1.678347)=5.035041.

Rounded to three decimal places, the simulated variate is 5.035. As a check, substituting $z\approx1.678347$ in $1-e^{-z}(1+z)$ gives approximately $0.5000$, the input uniform. The method avoids assuming that a gamma inverse has a simple closed form; the familiar shape-2 CDF here is a check on the numerical procedure.

Worked Example: Noninteger Shape and Scale Conversion

For a hypothetical model, let $X\sim\operatorname{Gamma}(2.5,4)$ and suppose the generated uniform is $U=0.584119$. Use numerical inversion to report $X$ to three decimal places.

Method: invert the unit-scale CDF, then apply the scale. Solve $P(2.5,z)=0.584119$. To see the CDF evaluation at $z=2.5$, the series has prefactor

\frac{e^{-2.5}(2.5)^{2.5}}{\Gamma(3.5)} \approx 0.244083

and the bracketed series is approximately $2.393115$. Their product is approximately $0.584119$. A numerical root search therefore returns $z\approx2.500$. For example, evaluations near $z=2.49$ and $z=2.51$ are approximately $0.5817$ and $0.5866$, respectively, which bracket the target.

The original variable has scale 4, so

X=4(2.500)\approx10.000.

The answer to three decimal places is 10.000. Substitution gives $F_X(10)=P(2.5,10/4)=P(2.5,2.5)\approx0.584119$, checking the inversion. The example illustrates that the numerical root is in unit-scale coordinates until the scale conversion is made.

Worked Example: Small Shape and a Bisection Accuracy Check

A hypothetical positive-amount model has $X\sim\operatorname{Gamma}(0.5,4)$, and the uniform draw is $U=0.5000$. Use numerical inversion and report the variate to three decimal places.

Method: bisection on the unit-scale CDF. For shape $0.5$, $P(0.5,z)=\operatorname{erf}(\sqrt z)$, where $\operatorname{erf}$ is the error function. At the initial bounds,

P(0.5,0.20)\approx0.4729,\qquad P(0.5,0.25)\approx0.5205.

Since $0.4729<0.5000<0.5205$, the root is between $0.20$ and $0.25$. The midpoint $z=0.225$ gives $P(0.5,0.225)\approx0.4978$, so the root is above $0.225$. Continued bisection yields $z\approx0.227468$. Apply the scale:

$$ X=4(0.227468)=0.909872. $$

Rounded to three decimal places, $X=0.910$. The check $P(0.5,0.227468)\approx0.5000$ confirms that the unit-scale root matches the uniform draw. This case also shows why a bracket is useful: it makes the direction of each refinement explicit even when the CDF evaluation uses a special function.

Accuracy, Assumptions, and Common Traps

Numerical inversion approximates a mathematical inverse, so the solver needs a stated stopping rule. One option is to stop when the bracket width is below a chosen tolerance; another is to stop when the CDF residual $|F_X(x)-U|$ is small. A small residual is useful, but it does not by itself guarantee a small error in $x$ when the density is very low. In the interior of the support, a local approximation explains the relationship:

F_X(x+\Delta x)-F_X(x)\approx f_X(x)\Delta x, \qquad \Delta x\approx\frac{\Delta F}{f_X(x)}.

Thus, when $f_X(x)$ is small, a small CDF error can correspond to a relatively large error in the variate. For ordinary simulation, choose tolerances appropriate to the application and retain enough numerical precision until the final reported value. For extremely small upper-tail probabilities, direct CDF evaluation can lose precision through subtraction from 1; use a numerical routine designed for the relevant tail if such probabilities are required.

The main parameterization trap is confusing scale $\theta$ with rate $\beta=1/\theta$. Another is to solve for $z$ and report it as $X$ without converting back to the specified scale.

Common mistake: In the shape-2 example, the numerical root $z\approx1.678347$ solves the unit-scale equation. Reporting $1.678$ as the requested variate silently omits the scale. The correction is $X=\theta z=3(1.678347)\approx5.035$. If the problem instead supplies a rate $\beta$, convert it to scale $\theta=1/\beta$ before applying the scale transformation.

The inverse-transform argument assumes the input is uniform on $(0,1)$, and the gamma parameters must satisfy $\alpha>0$ and $\theta>0$. A finite-precision generator may produce values extremely close to 0 or 1; an implementation must handle its numerical library's endpoint conventions rather than attempt to evaluate an infinite endpoint quantile. For typical draws away from the endpoints, a positive bracket and a reliable gamma CDF routine provide a transparent and checkable procedure.

Check Your Understanding

For each problem, use the stated shape-scale convention. Answers are rounded only as indicated.

A hypothetical model has $X\sim\operatorname{Gamma}(2,2)$ and $U=0.5000$. Use the shape-2 unit-scale root $z\approx1.678347$ to find $X$, to three decimals. (Answer: $3.357$)
For $X\sim\operatorname{Gamma}(1,5)$ and $U=0.2000$, invert the CDF and report $X$ to four decimals. (Answer: $1.1157$)
For $X\sim\operatorname{Gamma}(3,2)$, the unit-scale CDF is $P(3,z)=1-e^{-z}(1+z+z^2/2)$. If $U=1-5e^{-2}$, find $X$. (Answer: $4$)
A solver for $X\sim\operatorname{Gamma}(2,0.5)$ returns the unit-scale root $z\approx1.678347$ for $U=0.5000$. What variate should be reported to three decimals? (Answer: $0.839$)
A numerical search gives a root bracket $[1.24,1.25]$ for the unit-scale equation $P(\alpha,z)=U$, with $\theta=3$. What bracket does this imply for the requested gamma variate $X$, and what is its width? (Answer: $[3.72,3.75]$, width $0.03$)

What You'll Learn