4. Generating Random Variates

How to turn i.i.d. U(0,1) random numbers into observations from the distributions a model actually needs — the inverse transform and composition methods, and what makes a variate-generation algorithm good.

Goal. Study algorithms that produce observations (realizations) from some desired input distribution, given a stream of i.i.d. $U(0,1)$ random numbers.

Chapter 3 produced the raw material: uniform random numbers. This chapter turns them into random variates from whatever distribution the model calls for — e.g. the exponential interarrival and service times of an $M/M/1$ queue.

4.1. Structure

• It is essential to have a statistically reliable underlying RNG. Every algorithm below assumes its input really is i.i.d. $U(0,1)$; a defective generator corrupts every variate built from it.

Factors in choosing an algorithm.

1. Exactness. An algorithm is exact if the variates it produces have exactly the desired distribution, not merely approximately.

Example of an inexact method. To generate $Z \sim N(0,1)$:

1. Generate $u_1, \dots, u_{12}$ i.i.d. $U(0,1)$.
2. Return $Z = u_1 + u_2 + \cdots + u_{12} - 6$.

The moments come out right:

$$\mathbf{E}[Z] = 12 \cdot \tfrac{1}{2} - 6 = 0, \qquad \operatorname{Var}(Z) = 12 \cdot \tfrac{1}{12} = 1,$$

and by the CLT the sum is approximately normal. But it is not exact. Since each $u_i \in [0,1]$, the output always satisfies

$$Z \in [-6, 6],$$

whereas a true $N(0,1)$ variable has range $(-\infty, \infty)$. The algorithm can never produce a value beyond $\pm 6$, so it truncates both tails — exactly the region that matters most in reliability, risk, and rare-event studies. It is also wasteful: 12 random numbers per variate. Avoid it; use an exact method (§4.2.5) instead.

2. Efficiency.

• memory usage;
• execution time, including setup time (some algorithms precompute tables; that cost is paid once but must be counted).

3. Simplicity: easy to understand, implement, and debug. There is a trade-off between efficiency and simplicity — the fastest algorithm is rarely the clearest.

4. If possible, use exactly one random number per random variate. Reasons:

• Simplicity.
• Synchronization. When comparing two systems with common random numbers (§3.3, practical issue 3), we want the $i$-th variate in system A to be built from the same $u_i$ as in system B. If one system consumes a variable number of uniforms per variate, the two streams drift out of step and the comparison loses its sharpness.
• Monotonicity. With one uniform per variate and $X = F^{-1}(u)$, the variate is a monotone function of $u$. This is what makes variance-reduction techniques such as antithetic variates work.

4.2. RVG Algorithms

Direct methods — build the variate straight from the definition of the distribution:

1. Inverse transform — universally applicable in principle, but may be hard to implement.
2. Composition ⎫
3. Convolution ⎭ — apply only to input distributions of a specific form (mixtures and sums, respectively).

When these fail — that is, when $F^{-1}$ has no usable closed form and the distribution is neither a convenient mixture nor a convenient sum (the normal, gamma, and beta distributions are the standard examples) — fall back on:

4. Acceptance–rejection method.
5. Special techniques, e.g. for normal random variables.

4.2.1. Inverse Transform Method

Continuous case. Let $X$ be a continuous r.v. with c.d.f. $F$.

Algorithm.

1. Generate $u \sim U(0,1)$.
2. Return $X = F^{-1}(u)$.

The catch is step 2: $F^{-1}$ must be computed, and for many distributions (normal, gamma, beta) it has no closed form.

Proof that the output has c.d.f. $F$.

$$\begin{aligned} P(X \le x) &= P\left(F^{-1}(u) \le x\right) \\ &= P\left(F(F^{-1}(u)) \le F(x)\right) && \text{applying } F \text{ to both sides; valid since } F \text{ is non-decreasing} \\ &= P\left(u \le F(x)\right) \\ &= F(x) && \text{since } u \sim U(0,1). \end{aligned}$$

The last step is the definition of the uniform distribution: $P(u \le a) = a$ for any $a \in [0,1]$, and here $a = F(x)$ is indeed in $[0,1]$ because it is a probability. So the c.d.f. of the generated $X$ is $F$, exactly — the method is exact.

Example: $X \sim \text{Exp}(\lambda)$. Here $F(x) = 1 - e^{-\lambda x}$ for $x > 0$. Set $u = F(x)$ and solve:

$$u = 1 - e^{-\lambda x} \;\Longrightarrow\; e^{-\lambda x} = 1 - u \;\Longrightarrow\; x = -\frac{1}{\lambda}\ln(1 - u).$$

Algorithm.

1. Generate $u \sim U(0,1)$.
2. Return $X = -\dfrac{1}{\lambda}\ln(1-u)$.

Variant. Replace step 2 by “2′. Return $X = -\dfrac{1}{\lambda}\ln u$.”

This is valid because $1 - u \sim U(0,1)$ whenever $u \sim U(0,1)$, so the two versions produce the same distribution. It saves one subtraction per variate — a genuine but small gain in a tight loop.

Two cautions, though:

• Monotonicity is lost. Version 2 is increasing in $u$; version 2′ is decreasing. Under 2, a large $u$ gives a long interarrival time in both of two systems being compared; under 2′ the correspondence is reversed. Where synchronization or antithetic variates matter, prefer 2.
• The value $u = 0$. An LCG can output $u = 0$ (when $z_i = 0$), and $\ln 0$ is undefined. Version 2 never has this problem, since $1 - u \in (0, 1]$. If 2′ is used, $u = 0$ must be screened out.

Efficiency problem. The number of comparisons the chain performs equals $I$, the index it returns. Since the algorithm returns $x_i$ exactly when $I = i$, and $P(I = i) = P(X = x_i) = p(x_i)$, the number of comparisons $N$ is a random variable with the same distribution as the index, so by the definition of expectation,

$$\mathbf{E}[N] = \sum_{i=1}^{n} i \cdot P(N = i) = \sum_{i=1}^{n} i\, p(x_i).$$

If the probable values sit late in the list, this is close to $n$: tiny jumps at $x_1, x_2, \dots$ and the big jumps only at the end.

Two fixes.

1. Binary search. Since $F(x_1) \le F(x_2) \le \cdots \le F(x_n)$ is a sorted array, locate $u$ by bisection instead of scanning: $O(\log_2 n)$ comparisons instead of $O(n)$, whatever the probabilities.

2. Reorder the comparisons. Test the values in decreasing order of probability. Relabel so that $p_{(1)} \ge p_{(2)} \ge \cdots \ge p_{(n)}$ and build the cumulative sums for that new order; then $\mathbf{E}[N] = \sum_i i\, p_{(i)}$ is as small as possible. The reordering is done once in setup, so it costs nothing at run time.

Example. $n = 6$ with $p = (0.05,\, 0.05,\, 0.10,\, 0.10,\, 0.20,\, 0.50)$:

order of testing	$\mathbf{E}[\text{comparisons}]$
as given: $x_1, x_2, \dots, x_6$	$\sum_i i\,p_i = 4.85$
decreasing probability: $x_6, x_5, x_3, x_4, x_1, x_2$	$\sum_i i\,p_{(i)} = 2.15$

More than twice as fast, for no cost beyond sorting once.

Does reordering change the simulation? No. The output distribution is exactly the same — each $x_i$ is still returned with probability $p(x_i)$, because the intervals of $(0,1]$ assigned to the values have the same widths $p(x_i)$; only their positions along $(0,1]$ are permuted. The only thing that changes is the expected number of comparisons, i.e. the run time.

The one side effect: reordering destroys the monotone relationship between $u$ and $X$, so the output is no longer $F^{-1}(u)$ as a function of $u$. If synchronization or antithetic variates matter (see “Advantages” below), keep the natural order or use binary search, which preserves it.

Example: discrete uniform. $X \sim DU(1,100)$, i.e. $X \in \{1, 2, \dots, 100\}$ with $p(i) = \frac{1}{100}$ for every $i$, so $F(i) = \frac{i}{100}$.

The general algorithm becomes:

1. Generate $u \sim U(0,1)$.
2. if u ≤ 0.01 return X = 1; else if u ≤ 0.02 return X = 2; … ; else return X = 100

But here the chain can be replaced by a formula. The first $i$ with $u \le i/100$ is characterized by

$$\frac{i-1}{100} < u \le \frac{i}{100} \quad\Longleftrightarrow\quad i - 1 < 100u \le i \quad\Longleftrightarrow\quad i = \lceil 100u \rceil,$$

so step 2 collapses to

$$\boxed{\;X = \lceil 100 u \rceil\;}$$

— one multiplication and one ceiling, instead of up to 100 comparisons. (Equal probabilities are what make this possible; unequal ones need the search.)

Generalized inverse transform. To cover continuous, discrete, and mixed distributions in one statement:

1. Generate $u \sim U(0,1)$.
2. Return $X = \min\{x : F(x) \ge u\}$.

This is the same rule as before, written so it always makes sense. An ordinary inverse $F^{-1}$ requires $F$ to be strictly increasing and continuous; a discrete $F$ is flat in places (so no unique inverse) and jumps in places (so some values of $u$ have no $x$ with $F(x) = u$ at all). Taking the smallest $x$ whose $F(x)$ has reached $u$ handles both defects, and reduces to the usual $F^{-1}(u)$ when $F$ is continuous and strictly increasing.

Advantages and disadvantages of the inverse transform method.

Disadvantages.

1. Must invert the c.d.f., which may have no closed form (normal, gamma, beta).
2. May not be the most efficient approach even when the inverse exists.

Advantages.

3. Easy to generate from truncated distributions. Truncating means restricting $X$ to an interval $[a,b]$ and renormalizing:

   $$F^*(x) = \frac{F(x) - F(a)}{F(b) - F(a)}, \qquad a \le x \le b,$$

   which is the original density chopped off outside $[a,b]$ and scaled up so the area is again 1.

   Inverting $F^*$ costs nothing extra: generate $u \sim U(0,1)$ and return

   $$X = F^{-1}\Big(F(a) + u\,\big[F(b) - F(a)\big]\Big).$$

   In words: instead of drawing $u$ uniformly over all of $(0,1)$, draw it uniformly over the slice $\big(F(a), F(b)\big)$ of the vertical axis. Every draw is usable — no rejection, one uniform per variate.

4. Facilitates variance reduction when comparing alternative systems. Let $X_1 \sim F_1$ be the output of system 1 and $X_2 \sim F_2$ that of system 2, and suppose we want the difference in mean cost,

   $$\mathbf{E}[X_1] - \mathbf{E}[X_2] = \mathbf{E}[X_1 - X_2],$$

   so that $X_1 - X_2$ is an unbiased estimator of the gap. Generate both by inverse transform, $X_1 = F_1^{-1}(u_1)$ and $X_2 = F_2^{-1}(u_2)$. Since

   $$\operatorname{Var}(X_1 - X_2) = \operatorname{Var}(X_1) + \operatorname{Var}(X_2) - 2\operatorname{cov}(X_1, X_2),$$

   the precision of the comparison is controlled entirely by how we couple $u_1$ and $u_2$:

   coupling	$\operatorname{cov}(X_1, X_2)$	$\operatorname{Var}(X_1 - X_2)$
   1. $u_1 \perp u_2$ (independent)	$0$	$\operatorname{Var}(X_1) + \operatorname{Var}(X_2)$
   2. $u_1 = u_2$ (common random numbers)	$> 0$	smaller
   3. $u_1 = 1 - u_2$ (antithetic)	$< 0$	larger

   Case 2 is the useful one. Both $F_1^{-1}$ and $F_2^{-1}$ are non-decreasing, so a large $u$ pushes both $X_1$ and $X_2$ up and a small $u$ pushes both down: the two outputs move together, $\operatorname{cov}(X_1,X_2) > 0$, and the variance of the difference drops. Note $X_1 \ne X_2$ still, since $F_1 \ne F_2$. Monotonicity is exactly what the inverse transform provides and what a rejection-based method does not — hence “facilitates”.

   (Case 3 is not useless in general; antithetic variates reduce variance when estimating a single mean by averaging $X_1$ and $X_2$, where the negative correlation helps. It is only harmful here, where we take a difference.)

4.2.2. Composition

The name. The target distribution is composed of simpler ones — the way a mixture’s chemical composition lists its ingredients. The method is also called the mixture method: decompose $F$ into pieces you know how to sample from, pick a piece at random, then sample from it.

Setup. Use this when $F$ can be written as a convex combination of other c.d.f.s $F_1, F_2, \dots$ (finitely or infinitely many):

$$F(x) = \sum_{j \ge 1} p_j F_j(x), \qquad \sum_{j \ge 1} p_j = 1, \quad p_j \ge 0 .$$

Equivalently, in terms of densities or mass functions,

$$f(x) = \sum_{j \ge 1} p_j f_j(x) \qquad\text{or}\qquad p(x) = \sum_{j \ge 1} p_j\, p^{(j)}(x).$$

Algorithm.

1. Generate a positive random integer $J$ with $P(J = j) = p_j$ (a discrete variate — use the inverse transform of §4.2.1).
2. Return $X$ generated from the c.d.f. $F_J$.

Proof. Conditioning on $J$,

$$P(X \le x) = \sum_{j \ge 1} P(X \le x \mid J = j)\, P(J = j) = \sum_{j \ge 1} F_j(x)\, p_j = F(x). \qquad \blacksquare$$

Example: right trapezoidal distribution.

$$f(x) = \begin{cases} a + 2(1-a)x, & 0 \le x \le 1 \quad (0 < a < 1), \\ 0, & \text{otherwise.}\end{cases}$$

Split the formula into its two terms and read off the weights:

$$f(x) = \underbrace{a}_{p_1} \cdot \underbrace{1}_{f_1(x)} \;+\; \underbrace{(1-a)}_{p_2} \cdot \underbrace{2x}_{f_2(x)} ,$$

where $f_1$ is the $U(0,1)$ density and $f_2(x) = 2x$ on $[0,1]$ is a triangular density with $F_2(x) = x^2$, so $F_2^{-1}(u) = \sqrt{u}$.

Algorithm.

1. Generate $u_1, u_2$ independent $U(0,1)$.
2. If $u_1 \le a$, return $X = u_2$; else return $X = \sqrt{u_2}$.

Here $u_1$ chooses the component ($J = 1$ with probability $a$) and $u_2$ generates from it.

Example: symmetric triangular distribution on $[-1,1]$.

$$f(x) = \begin{cases} 1 + x, & -1 \le x \le 0, \\ 1 - x, & 0 \le x \le 1, \\ 0, & \text{otherwise.}\end{cases}$$

Each half of the triangle has area $\tfrac12$, so take $p_1 = p_2 = \tfrac12$ with

$$f_1(x) = 2(1+x) \ \text{ on } [-1,0], \qquad f_2(x) = 2(1-x) \ \text{ on } [0,1]$$

(each doubled so that it integrates to 1 on its own half). For the right piece, $F_2(x) = 2x - x^2$, and solving $u = 2x - x^2$ gives $x = 1 - \sqrt{1-u}$, i.e. $1 - \sqrt{u}$ after replacing $1-u$ by $u$. The left piece is its mirror image.

Algorithm (composition).

1. Generate $u_1, u_2$ independent $U(0,1)$.
2. If $u_1 \le \tfrac12$, return $X = \sqrt{u_2} - 1$; else return $X = 1 - \sqrt{1 - u_2}$.

The same distribution by inverse transform. The c.d.f. is

$$F(x) = \begin{cases} \dfrac{(1+x)^2}{2}, & -1 \le x \le 0, \\[6pt] 1 - \dfrac{(1-x)^2}{2}, & 0 \le x \le 1,\end{cases}$$

and $F(0) = \tfrac12$, so $u < \tfrac12$ corresponds to the left half. Solving $u = \frac{(1+x)^2}{2}$ gives $x = -1 + \sqrt{2u}$; solving $u = 1 - \frac{(1-x)^2}{2}$ gives $x = 1 - \sqrt{2(1-u)}$.

1. Generate $u \sim U(0,1)$.
2. If $u < \tfrac12$, return $X = -1 + \sqrt{2u}$; else return $X = 1 - \sqrt{2(1-u)}$.

Efficiency comparison (expected operations per variate):

	$u$’s	comparisons	additions	multiplications	square roots
inverse	1	1	1.5	1	1
composition	2	1	1.5	0	1

How the counts are obtained. Each row averages the two branches, which are taken with probability $\tfrac12$ each; a subtraction counts as an addition. For the inverse method, the left branch $-1 + \sqrt{2u}$ uses one multiplication ($2u$), one square root and one addition, while the right branch $1 - \sqrt{2(1-u)}$ uses one multiplication ($2 \times$), one square root and two additions ($1-u$ and $1 - \sqrt{\cdot}$) — so the average is $\tfrac12(1) + \tfrac12(2) = 1.5$ additions, with 1 multiplication and 1 square root either way, plus the single comparison $u < \tfrac12$ and the single uniform. For composition, $\sqrt{u_2} - 1$ has one addition and $1 - \sqrt{1-u_2}$ has two, giving the same $1.5$ average and the same one square root and one comparison, but no multiplication — the factor 2 never appears, because doubling the density was absorbed into choosing a half. The cost is that $u_1$ is spent on the coin toss and $u_2$ on the value, so two uniforms are consumed instead of one.

Conclusion. The two methods differ in exactly two places: composition saves one multiplication but spends one extra random number. Generating a random number is far more expensive than a multiplication, so the inverse transform is more efficient here. (It also keeps $X$ monotone in $u$, which composition does not — see §4.2.1, advantage 4.)

4.2.3. Convolution

Setup. Let $Y_1, Y_2, \dots, Y_m$ be i.i.d. with common c.d.f. $G$, and define

$$X = Y_1 + Y_2 + \cdots + Y_m .$$

The c.d.f. $F$ of $X$ is called the $m$-fold convolution of $G$ with itself. Use this method whenever the target distribution is known to arise as such a sum.

Algorithm (to generate from $F$):

1. Generate $Y_1, \dots, Y_m$ independently from $G$.
2. Return $X = Y_1 + Y_2 + \cdots + Y_m$.

Why the proof is trivial. Unlike the inverse transform or composition, there is nothing to verify: $F$ was defined as the distribution of the sum, and the algorithm literally forms that sum from independent draws. Correctness is by construction. The work in this method lies elsewhere — in recognizing that the desired $F$ is a convolution of something samplable.

Example: sum of two exponentials. Let $X_1 \perp X_2 \sim \text{Exp}(\lambda)$ and $Z = X_1 + X_2$. Then

$$\begin{aligned} F_Z(z) = P(Z \le z) &= P(X_1 + X_2 \le z) \\ &= \int_0^\infty P\!\left(X_1 + X_2 \le z \mid X_2 = x\right) f_{X_2}(x)\,dx && \text{(1)} \\ &= \int_0^\infty P\!\left(X_1 \le z - x\right) \lambda e^{-\lambda x}\,dx && \text{(2)} \\ &= \int_0^z \left(1 - e^{-\lambda(z-x)}\right) \lambda e^{-\lambda x}\,dx && \text{(3)} \\ &= \int_0^z \lambda e^{-\lambda x}\,dx \;-\; \int_0^z \lambda e^{-\lambda z}\,dx && \text{(4)} \\ &= \left(1 - e^{-\lambda z}\right) - \lambda z\, e^{-\lambda z} && \text{(5)} \\ &= 1 - \sum_{i=0}^{1} e^{-\lambda z}\frac{(\lambda z)^i}{i!} . && \text{(6)} \end{aligned}$$

Explanation of each step.

1. Condition on $X_2$ (law of total probability, continuous version). To find $P(A)$, split according to the value $x$ that $X_2$ takes, weight each case by the density $f_{X_2}(x)$, and integrate over all $x$ — the continuous analogue of $P(A) = \sum_j P(A \mid B_j)\,P(B_j)$. The range is $[0, \infty)$ because $X_2 \ge 0$. The purpose is to turn a two-variable event into a one-variable one.
2. Use the condition, then drop it. Given $X_2 = x$, the event $X_1 + X_2 \le z$ is $X_1 + x \le z$, i.e. $X_1 \le z - x$. Then, because $X_1 \perp X_2$, knowing $X_2 = x$ says nothing about $X_1$, so $P(X_1 \le z - x \mid X_2 = x) = P(X_1 \le z - x)$ — the condition disappears. Also substitute $f_{X_2}(x) = \lambda e^{-\lambda x}$.
3. Substitute the c.d.f. of $X_1$. $P(X_1 \le z - x) = F_{X_1}(z - x) = 1 - e^{-\lambda(z-x)}$, valid when $z - x \ge 0$. When $x > z$ the argument is negative and $F_{X_1} = 0$ ($X_1$ is never negative), so the integrand vanishes there; hence the upper limit drops from $\infty$ to $z$.
4. Expand and simplify. Multiply out $(1 - e^{-\lambda(z-x)})\,\lambda e^{-\lambda x}$ and split into two integrals. In the second, $e^{-\lambda(z-x)} \cdot e^{-\lambda x} = e^{-\lambda z + \lambda x - \lambda x} = e^{-\lambda z}$: the $x$ cancels, leaving a constant.
5. Evaluate. First integral: $\int_0^z \lambda e^{-\lambda x}\,dx = \left[-e^{-\lambda x}\right]_0^z = 1 - e^{-\lambda z}$. Second: the integrand is the constant $\lambda e^{-\lambda z}$, so the integral is that constant times the interval length $z$.
6. Rewrite as a sum. Factor out $e^{-\lambda z}$: $1 - e^{-\lambda z}(1 + \lambda z)$. The two terms in the bracket are $\frac{(\lambda z)^0}{0!} = 1$ and $\frac{(\lambda z)^1}{1!} = \lambda z$, i.e. $i = 0, 1$ of $\sum_i \frac{(\lambda z)^i}{i!}$. Writing it this way exposes the pattern: for a sum of $n$ exponentials the same calculation, repeated, gives the terms $i = 0, \dots, n-1$.

Generalization: the Erlang distribution. Let $X_1, \dots, X_n$ be i.i.d. $\text{Exp}(\lambda)$ and $Z = X_1 + \cdots + X_n$. Repeating the argument gives

$$F_Z(z) = 1 - \sum_{i=0}^{n-1} e^{-\lambda z}\frac{(\lambda z)^i}{i!}, \qquad z \ge 0,$$

and $Z$ is called an Erlang$(n, \lambda)$ variable — the special case of the Gamma$(n, \lambda)$ distribution with integer shape parameter.

Reading the formula. The sum is exactly $P(N(z) \le n-1)$ for a Poisson variable $N(z)$ with mean $\lambda z$. This is the Poisson-process statement of §2.5: the $n$-th arrival occurs by time $z$ if and only if at least $n$ arrivals have occurred by then. So “waiting for the $n$-th event” and “counting events” are two views of the same process.

Algorithm for Erlang$(n, \lambda)$.

1. Generate $Y_1, \dots, Y_n$ i.i.d. $\text{Exp}(\lambda)$ (each by inverse transform, §4.2.1).
2. Return $X = Y_1 + Y_2 + \cdots + Y_n$.

This is exact and trivially correct, but it costs $n$ uniforms and $n$ logarithms per variate. One improvement is free:

$$X = \sum_{i=1}^{n} \left(-\frac{1}{\lambda}\ln u_i\right) = -\frac{1}{\lambda}\ln\!\left(\prod_{i=1}^{n} u_i\right),$$

so multiply the uniforms first and take a single logarithm. The $n$ uniforms are still needed, which is why convolution becomes unattractive for large $n$ — and why the Gamma distribution with non-integer shape needs a different method (acceptance–rejection, §4.2.4).

Composition vs. convolution. Both write the target in terms of simpler distributions, but in opposite ways:

	composition	convolution
decomposition	$F = \sum_j p_j F_j$ (weighted mixture)	$X = Y_1 + \cdots + Y_m$ (sum)
what you do	pick one component at random, sample it	sample all components, add them
uniforms used	1 for the choice + those for one component	those for all $m$ components

Example: symmetric triangular distribution, third method. Recall (§4.2.2)

$$f(x) = \begin{cases} x + 1, & -1 \le x < 0, \\ -x + 1, & 0 \le x \le 1, \\ 0, & \text{otherwise.} \end{cases}$$

Let $u_1 \perp u_2 \sim U(0,1)$ and $Z = u_1 + u_2$. The pair $(u_1, u_2)$ is uniform on the unit square (joint density $= 1$), so a probability is an area:

$$F_Z(z) = P(Z \le z) = P(u_1 + u_2 \le z) = \iint\limits_{\substack{0 \le u_1 \le 1,\; 0 \le u_2 \le 1 \\ u_1 + u_2 \le z}} 1 \; du_1\, du_2 = \text{area of the part of the square below the line } u_2 = z - u_1 .$$

• Case $0 \le z \le 1$. The line cuts off a right triangle at the origin with legs $z$ and $z$:

  $$F_Z(z) = \tfrac{1}{2} z^2 .$$

• Case $1 \le z \le 2$. Now the line cuts off the top-right corner, a right triangle with legs $2 - z$. Shaded area $=$ whole square minus that corner:

  $$F_Z(z) = 1 - \tfrac{1}{2}(2 - z)^2 = -\tfrac{1}{2}z^2 + 2z - 1 .$$

Differentiating gives the density of $Z$:

$$f_Z(z) = F_Z'(z) = \begin{cases} z, & 0 \le z \le 1, \\ -z + 2, & 1 \le z \le 2, \end{cases}$$

a triangle on $[0, 2]$ with peak at $z = 1$. Subtracting 1 slides it to $[-1, 1]$: $X = Z - 1$ has $f_X(x) = f_Z(x + 1)$, which is $x + 1$ on $[-1, 0]$ and $-(x+1) + 2 = 1 - x$ on $[0, 1]$ — exactly the target.

Algorithm (convolution).

1. Generate $u_1 \perp u_2 \sim U(0,1)$.
2. Return $X = u_1 + u_2 - 1$.

Why this beats both earlier methods. Add the new row to the table from §4.2.2:

	$u$’s	comparisons	additions	multiplications	square roots
inverse	1	1	1.5	1	1
composition	2	1	1.5	0	1
convolution	2	0	2	0	0

Convolution uses the same two uniforms as composition but eliminates the comparison, the multiplication, and above all the square root, at the price of half an extra addition. A square root is by far the most expensive operation in the table — many times the cost of an addition, and more than the cost of drawing a uniform — so removing it outweighs the second uniform, and convolution is the most efficient of the three for this distribution. (Against inverse transform, the trade is one extra uniform plus half an addition in exchange for dropping a comparison, a multiplication and a square root; the square root alone settles it.)

4.2.4. Acceptance–Rejection (AR)

Where it sits. Methods 1–3 (§4.2.1–4.2.3) are direct: they compute $X$ from the uniforms by a formula. Acceptance–rejection is the indirect method from the list in §4.2: it does not compute $X$; it generates candidates from an easier distribution and keeps only some of them. The arrow in the lecture outline (“rejection → indirect”) just records that classification.

When to use it. When the desired p.d.f. $f$ has a form too complex to invert or decompose.

Idea. Use a simple, easy-to-sample p.d.f. — a surrogate — that approximates $f$ from above, then correct for the difference by discarding a fraction of the draws.

Setup. Choose a majorizing function $t$ with

$$t(x) \ge f(x) \quad \text{for all } x .$$

Then

$$c = \int_{-\infty}^{\infty} t(x)\,dx \;\ge\; \int_{-\infty}^{\infty} f(x)\,dx = 1,$$

so $t$ is not a density (its area exceeds 1), but its normalization

$$r(x) = \frac{t(x)}{c}$$

is a valid p.d.f. ($r \ge 0$ and $\int r = c/c = 1$). $r$ is the surrogate we actually sample from.

AR algorithm.

1. Generate $Y \sim r$.
2. Generate $u \sim U(0,1)$, independent of $Y$.
3. If $u \le \dfrac{f(Y)}{t(Y)}$, return $X = Y$ (accept).
4. Otherwise discard $Y$ and go back to step 1.

On step 1: $Y$ is generated from the surrogate $r$, not from $f$ — that is the whole point. Since $r$ was chosen to be simple, this is done by one of the direct methods (typically inverse transform). $Y$ is only a candidate for $X$; steps 3–4 decide whether it becomes one.

On step 3: the ratio $f(Y)/t(Y)$ lies in $[0,1]$ because $t \ge f$. Comparing it with a uniform means: accept $Y$ with probability $f(Y)/t(Y)$. Where $t$ hugs $f$ this is near 1 (almost always accept); where $t$ is far above $f$ it is small (usually reject). The rejections thin out exactly the regions where $r$ over-samples relative to $f$.

Proof that the output has c.d.f. $F$. The returned $X$ is a candidate $Y$ given that it was accepted, so

$$P(X \le x) = P(Y \le x \mid Y \text{ accepted}) = \frac{P(Y \le x,\; Y \text{ accepted})}{P(Y \text{ accepted})} .$$

We compute the denominator first (it is also the efficiency, comment 2 below):

$$\begin{aligned} P(Y \text{ accepted}) &= P\!\left(u \le \tfrac{f(Y)}{t(Y)}\right) \\ &= \int_{-\infty}^{\infty} P\!\left(u \le \tfrac{f(Y)}{t(Y)} \;\middle|\; Y = y\right) r(y)\,dy && \text{(condition on } Y\text{)} \\ &= \int_{-\infty}^{\infty} P\!\left(u \le \tfrac{f(y)}{t(y)}\right) r(y)\,dy && (u \perp Y \text{, so the condition drops}) \\ &= \int_{-\infty}^{\infty} \frac{f(y)}{t(y)} \cdot \frac{t(y)}{c}\,dy && (P(u \le a) = a;\; r = t/c) \\ &= \frac{1}{c}\int_{-\infty}^{\infty} f(y)\,dy = \frac{1}{c} . \end{aligned}$$

The numerator is the same calculation with the extra event $Y \le x$, which after conditioning on $Y = y$ just restricts the integral to $y \le x$:

$$P(Y \le x,\; Y \text{ accepted}) = \int_{-\infty}^{x} \frac{f(y)}{t(y)} \cdot \frac{t(y)}{c}\,dy = \frac{1}{c}\int_{-\infty}^{x} f(y)\,dy = \frac{F(x)}{c} .$$

Dividing,

$$P(X \le x) = \frac{F(x)/c}{1/c} = F(x). \qquad\blacksquare$$

Note the mechanism: $t(y)$ cancels between the acceptance probability $f/t$ and the surrogate density $t/c$, and $c$ cancels between numerator and denominator. Whatever shape $t$ has, the output is exactly $f$ — the choice of $t$ affects only speed, never correctness.

Comments.

1. Generating from $r$ must be much easier than generating from $f$; otherwise there is no point in using a surrogate. In practice, choose simple majorizing functions — piecewise constant or piecewise linear.

2. The algorithm loops until a pair $(u, Y)$ with $u \le f(Y)/t(Y)$ turns up. From the proof,

   $$P(Y \text{ accepted}) = \frac{1}{c} = \frac{1}{\text{area under } t}, \qquad P(\text{reject}) = 1 - \frac{1}{c}.$$

   Why the average number of trials is $c$: each pass through steps 1–3 is an independent trial with success probability $p = 1/c$ (fresh $Y$ and $u$ each time). The number of trials until the first success is therefore Geometric$(p)$, whose mean is $1/p = c$. So $c$ uniforms-pairs are consumed per variate on average — a large $c$ means a slow generator.

3. $P(Y \text{ accepted}) = 1/c$ is called the efficiency of the AR method. We want it as close to 1 as possible, i.e. $t$ as tight around $f$ as possible — but a tighter $t$ is usually a more complicated one, harder to sample from. This is the same simplicity vs. efficiency trade-off as in §4.1.

Example. $f(x) = \frac{3}{2}x^2$ on $(-1, 1)$ (and $0$ elsewhere). Check: $\int_{-1}^{1} \frac{3}{2}x^2\,dx = 1$. The maximum of $f$ is $\frac32$, at $x = \pm 1$.

First attempt: a flat majorizing function. Take $t(x) = \frac{3}{2}$ on $(-1,1)$ — the smallest constant that stays above $f$. Then $c = \frac{3}{2} \cdot 2 = 3$ and $r(x) = t(x)/3 = \frac{1}{2}$ on $(-1,1)$, i.e. $r$ is the $U(-1,1)$ density.

1. Generate $Y \sim U(-1, 1)$ (e.g. $Y = 2u' - 1$ with $u' \sim U(0,1)$).
2. Generate $u \sim U(0,1)$, independent of $Y$.
3. If $u \le Y^2$, return $X = Y$. $\;\left(\dfrac{f(Y)}{t(Y)} = \dfrac{\frac32 Y^2}{\frac32} = Y^2\right)$
4. Otherwise go back to step 1.

Efficiency $= 1/c = 1/3$: on average three candidates (six uniforms) per accepted variate. The waste is visible in the figure — the rectangle under $t$ has area 3 while the bowl under $f$ has area 1, and everything in between is rejected.

Improvement: a tighter majorizing function. We need $t \ge f$ but with less area. On $[-1,1]$ we have $|x| \ge x^2$, so

$$t(x) = \tfrac{3}{2}\,|x| \;\ge\; \tfrac{3}{2}\,x^2 = f(x),$$

with equality at $x = 0, \pm 1$. Its area is $c = \int_{-1}^{1} \frac32|x|\,dx = \frac32$, so the efficiency rises to $1/c = 2/3$: on average $1.5$ candidates per variate instead of $3$. The surrogate is $r(x) = t(x)/c = |x|$ on $(-1,1)$, a symmetric “V” density, easy to sample by inverse transform: its right half has c.d.f. $x^2$, so $|Y| = \sqrt{u'}$, and a fair coin picks the sign.

1. Generate $u_1, u_2 \sim U(0,1)$; set $Y = \sqrt{u_2}$ if $u_1 \le \frac12$, else $Y = -\sqrt{u_2}$. Then $Y \sim r$.
2. Generate $u \sim U(0,1)$, independent of $Y$.
3. If $u \le |Y|$, return $X = Y$. $\;\left(\dfrac{f(Y)}{t(Y)} = \dfrac{\frac32 Y^2}{\frac32 |Y|} = |Y|\right)$
4. Otherwise go back to step 1.

Each trial now costs three uniforms and a square root instead of two uniforms, but only $1.5$ trials are needed on average instead of $3$ — the trade-off of comment 3 in miniature.