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.

4.2.5. Special Techniques

There is no general recipe here. A special technique exploits some property that the target distribution happens to have — a relation to another distribution we already know how to sample. Two small examples now; the big one (the normal distribution) gets its own treatment in §4.3.

Example: chi-square with one degree of freedom. Recall from §2.6.3 that if $Y \sim N(0,1)$ then $Y^2 \sim \chi^2_1$. So:

1. Generate $Y \sim N(0,1)$ (§4.3).
2. Return $X = Y^2$.

Combining with convolution: $\chi^2_k$ is the sum of $k$ independent $\chi^2_1$ variables, so $X = Y_1^2 + \cdots + Y_k^2$ with $Y_i$ i.i.d. $N(0,1)$.

Example: geometric by Bernoulli trials. Let $X \sim \text{Geo}(p)$ be the number of trials needed until the first success, each trial succeeding with probability $p$. Simulate the trials literally:

1. Set $X = 1$.
2. Generate $Y \sim \text{Bernoulli}(p)$ (i.e. $u \sim U(0,1)$ and $Y = 1$ if $u \le p$).
3. If $Y = 1$, return $X$.
4. Otherwise set $X = X + 1$ and go back to step 2.

This is exact and transparent, but it consumes one uniform per trial, and the expected number of trials is $\mathbf{E}[X] = 1/p$. For $p = 0.01$ that is 100 uniforms per variate. §4.4 gives a one-uniform formula instead.

4.3. Generating Continuous Random Variates

§4.2 presented the methods abstractly. This section collects the concrete algorithms actually used for the common continuous distributions.

Uniform $U(a,b)$. Inverse transform: $F(x) = \frac{x-a}{b-a}$ on $[a,b]$, so $u = \frac{x-a}{b-a} \iff x = a + (b-a)u$.

1. Generate $u \sim U(0,1)$.
2. Return $X = a + (b-a)\,u$.

Exponential $\text{Exp}(\lambda)$. From §4.2.1:

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

Erlang$(m, \lambda)$. From §4.2.3, $X = Y_1 + \cdots + Y_m$ with $Y_i$ i.i.d. $\text{Exp}(\lambda)$:

1. Generate $u_1, \dots, u_m \sim U(0,1)$.
2. Return $X = -\dfrac{1}{\lambda}\ln\!\left(\prod_{i=1}^m u_i\right)$.

(Gamma with non-integer shape, Weibull, beta, etc. need special-purpose methods, mostly acceptance–rejection; not covered here.)

4.3.1. Normal distribution: the Box–Muller method

The normal c.d.f. has no closed-form inverse, so inverse transform is out. Composition and convolution do not apply. The trick is a special technique: look at two independent standard normals at once and switch to polar coordinates.

The key fact. Let $X_1 \perp X_2 \sim N(0,1)$ and write the point $(X_1, X_2)$ in polar form:

$$R^2 = X_1^2 + X_2^2, \qquad \Theta = \tan^{-1}\!\left(\frac{X_2}{X_1}\right).$$

Then $R^2$ and $\Theta$ are independent, with

$$R^2 \sim \text{Exp}\!\left(\tfrac12\right) \ \text{(rate } \tfrac12\text{, mean } 2), \qquad \Theta \sim U(0, 2\pi).$$

Why. The joint density of $(X_1, X_2)$ is the product of two $N(0,1)$ densities:

$$f_{X_1,X_2}(x_1, x_2) = \frac{1}{\sqrt{2\pi}}e^{-x_1^2/2}\cdot\frac{1}{\sqrt{2\pi}}e^{-x_2^2/2} = \frac{1}{2\pi}\exp\!\left(-\frac{x_1^2 + x_2^2}{2}\right).$$

It depends on $(x_1, x_2)$ only through $x_1^2 + x_2^2 = r^2$ — the density is rotationally symmetric, which is why the angle will come out uniform. To change variables to $(S, \Theta) = (R^2, \Theta)$ we need the Jacobian of the map $g_1(x_1,x_2) = x_1^2 + x_2^2$, $g_2(x_1,x_2) = \tan^{-1}(x_2/x_1)$. Using $\frac{d}{dx}\tan^{-1}x = \frac{1}{1+x^2}$,

$$J = \begin{vmatrix} \partial g_1/\partial x_1 & \partial g_1/\partial x_2 \\ \partial g_2/\partial x_1 & \partial g_2/\partial x_2 \end{vmatrix} = \begin{vmatrix} 2x_1 & 2x_2 \\ -\dfrac{x_2}{x_1^2+x_2^2} & \dfrac{x_1}{x_1^2+x_2^2} \end{vmatrix} = \frac{2x_1^2 + 2x_2^2}{x_1^2 + x_2^2} = 2 .$$

The change-of-variables formula gives $f_{S,\Theta}(s,\theta) = f_{X_1,X_2}(x_1,x_2)\,|J|^{-1}$, so

$$f_{S,\Theta}(s, \theta) = \frac{1}{2\pi}e^{-s/2}\cdot\frac{1}{2} = \underbrace{\tfrac{1}{2}e^{-s/2}}_{\text{Exp}(1/2)\text{ density in } s}\;\cdot\;\underbrace{\tfrac{1}{2\pi}}_{U(0,2\pi)\text{ density in }\theta}, \qquad s \ge 0,\ 0 \le \theta < 2\pi .$$

The joint density factors into a function of $s$ alone times a function of $\theta$ alone, which is exactly the statement that $S = R^2$ and $\Theta$ are independent with those marginals.

Run it backwards. We know how to generate both pieces: $R^2 \sim \text{Exp}(\tfrac12)$ is $-2\ln u_1$ by inverse transform (rate $\lambda = \tfrac12$ gives $-\frac{1}{\lambda}\ln u = -2\ln u$), and $\Theta = 2\pi u_2$. Converting back from polar to Cartesian coordinates then yields two independent standard normals.

Algorithm (Box–Muller).

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

$$X_1 = \sqrt{-2\ln u_1}\,\cos(2\pi u_2), \qquad X_2 = \sqrt{-2\ln u_1}\,\sin(2\pi u_2).$$

Two uniforms in, two independent $N(0,1)$ variates out — exact, no rejection.

Two cautions.

• The uniforms must really be independent. If $u_1, u_2$ are consecutive outputs of a poor LCG (recall RANDU, §3.2.3), the pairs $(u_1, u_2)$ lie on a few lines in the unit square, and after the transformation the points $(X_1, X_2)$ fall on a visible spiral in the plane. In practice one uses separate streams for $u_1$ and $u_2$ (§3.3, practical issue 3).
• Cost. One logarithm, one square root, one sine and one cosine per pair. The trigonometric functions are expensive, which motivates the next algorithm.

4.3.2. Normal distribution: the polar method (acceptance–rejection)

The polar method avoids $\sin$ and $\cos$ by generating the angle geometrically instead of computing it.

Algorithm (polar method).

1. Generate $V_1 \perp V_2 \sim U(-1, 1)$ (i.e. $V_i = 2u_i - 1$).
2. Let $S = V_1^2 + V_2^2$. If $S > 1$ (or $S = 0$), go back to step 1.
3. Otherwise let

$$W = \sqrt{\frac{-2\ln S}{S}}$$

and return $X_1 = V_1 W$, $X_2 = V_2 W$.

Why it works. The pair $(V_1, V_2)$ is uniform on the square $[-1,1]^2$; conditioning on $S \le 1$ makes it uniform on the unit disk. For a uniform point in the disk, the angle $\Theta$ is $U(0,2\pi)$ and the squared radius $S$ is $U(0,1)$ (because $P(S \le s) = \text{area of disk of radius }\sqrt{s} = \pi s / \pi = s$), independently of each other. Now compare with Box–Muller: we need $\cos\Theta$, $\sin\Theta$, and $\sqrt{-2\ln u}$ for a uniform $u$. But

$$\cos\Theta = \frac{V_1}{\sqrt S}, \qquad \sin\Theta = \frac{V_2}{\sqrt S}, \qquad \text{and } S \text{ itself is the uniform we need, so } R = \sqrt{-2\ln S}.$$

Multiplying out,

$$X_1 = R\cos\Theta = \sqrt{-2\ln S}\cdot\frac{V_1}{\sqrt S} = V_1\sqrt{\frac{-2\ln S}{S}} = V_1 W,$$

and likewise $X_2 = V_2 W$. The cosine and sine have been replaced by a division.

Efficiency. $P(\text{accept}) = P(S \le 1) = \dfrac{\text{area of unit disk}}{\text{area of square}} = \dfrac{\pi}{4} \approx 0.785$, so on average $4/\pi \approx 1.27$ attempts, i.e. about $2.55$ uniforms per pair of normals. The price of avoiding trigonometry is that roughly one attempt in five is thrown away, and the method no longer uses a fixed number of uniforms per variate (bad for synchronization, §4.1).

Non-standard normal. If $Y \sim N(0,1)$ then $X = \mu + \sigma Y \sim N(\mu, \sigma^2)$. Generate $Y$ by either method above and shift and scale.

4.3.3. Lognormal distribution

$X$ is lognormal, $X \sim LN(\mu, \sigma^2)$, if $\ln X \sim N(\mu, \sigma^2)$ — equivalently $X = e^Y$ with $Y \sim N(\mu, \sigma^2)$. (Note $\mu, \sigma^2$ are the mean and variance of $\ln X$, not of $X$.)

1. Generate $Y \sim N(\mu, \sigma^2)$.
2. Return $X = e^Y$.

4.3.4. Truncated distributions

Suppose we want $X$ restricted to an interval $[a, b]$. There are two different things “truncate” can mean:

• Redistribute the mass that falls outside $[a,b]$ proportionally over $[a,b]$ — the conditional distribution of $X$ given $a \le X \le b$.
• Pile up the outside mass at the endpoints: values below $a$ become $a$, values above $b$ become $b$. This is clipping, $X^* = \min(\max(X, a), b)$, and it creates atoms at $a$ and $b$.

The first is the usual meaning. Its c.d.f. is

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

with density $\tilde f(x) = \dfrac{f(x)}{F(b) - F(a)}$ on $[a,b]$ — the original density scaled up by the constant $1/(F(b)-F(a))$ so that it integrates to 1 again. (Check: $\tilde F(x) = \int_a^x \frac{f(y)}{F(b)-F(a)}\,dy = \frac{F(x)-F(a)}{F(b)-F(a)}$.)

Inverting $\tilde F$ is one line:

$$u = \frac{F(x) - F(a)}{F(b) - F(a)} \;\Longrightarrow\; F(x) = F(a) + u\,[F(b) - F(a)] \;\Longrightarrow\; x = F^{-1}\big(F(a) + u\,[F(b)-F(a)]\big).$$

1. Generate $u \sim U(0,1)$.
2. Return $X = F^{-1}\big(F(a) + u\,[F(b) - F(a)]\big)$.

This is the “advantage 3” of §4.2.1 in detail. The alternative — generate from $F$ and reject anything outside $[a,b]$ — is acceptance–rejection with efficiency $F(b) - F(a)$, which can be tiny if the interval is in a tail; the inverse-transform version uses every uniform.

4.3.5. More continuous distributions

A short tour of the other continuous distributions the textbook covers.

Inverse transform, closed form.

• Weibull$(\alpha, \beta)$: $F(x) = 1 - e^{-(x/\beta)^\alpha}$, so $F^{-1}(u) = \beta\,[-\ln(1-u)]^{1/\alpha}$. Algorithm: generate $u$, return $X = \beta(-\ln u)^{1/\alpha}$. (Equivalently: an $\text{Exp}$ variable with mean $\beta^\alpha$, raised to the power $1/\alpha$.)

• Log-logistic$(\alpha, \beta)$: $F^{-1}(u) = \beta\,[u/(1-u)]^{1/\alpha}$.

• Triangular$(0, 1, m)$ (mode $m \in (0,1)$):

  $$F^{-1}(u) = \begin{cases} \sqrt{m\,u}, & 0 \le u \le m, \\ 1 - \sqrt{(1-m)(1-u)}, & m < u \le 1. \end{cases}$$

  Generate $u$; if $u \le m$ return $\sqrt{mu}$, otherwise return $1 - \sqrt{(1-m)(1-u)}$. For $\text{triang}(a, b, m)$ rescale: $X' = a + (b-a)X$ with mode $(m-a)/(b-a)$.

  The textbook’s “Why?”: in the second branch one may not replace $1 - u$ by $u$. The usual excuse — “$u$ and $1-u$ have the same distribution” — applies to a uniform used on its own. Here the same $u$ also decides the branch: the second formula is only ever evaluated for $u \in (m, 1]$, and on that set $1 - u$ ranges over $[0, 1-m)$ while $u$ ranges over $(m, 1]$ — different sets, different distributions. Making the swap gives the wrong answer (in a check with $m = 0.3$, $P(X \le m)$ came out $0.60$ instead of $0.30$). To swap, you would have to swap in both places: test $u \ge 1-m$ and use $1-u$ in the first branch as well.

Special properties.

• Student’s $t_k$: if $Y \sim N(0,1)$ and $Z \sim \chi^2_k$ are independent, $X = Y/\sqrt{Z/k} \sim t_k$. Generate $Y$ (§4.3.1) and $Z$ (§4.2.5) and combine.
• $F_{k_1, k_2}$: with independent $Z_1 \sim \chi^2_{k_1}$, $Z_2 \sim \chi^2_{k_2}$, $X = (Z_1/k_1)/(Z_2/k_2)$.
• Beta$(\alpha_1, \alpha_2)$ on $[0,1]$: with independent $Y_1 \sim \text{Gamma}(\alpha_1, 1)$, $Y_2 \sim \text{Gamma}(\alpha_2, 1)$, $X = Y_1/(Y_1 + Y_2)$. Special cases need no gamma: $\text{Beta}(\alpha_1, 1)$ has $F(x) = x^{\alpha_1}$, so $X = u^{1/\alpha_1}$; $\text{Beta}(1,1) = U(0,1)$; and $1 - X \sim \text{Beta}(\alpha_2, \alpha_1)$.
• Gamma$(\alpha, \beta)$ with non-integer $\alpha$: no closed-form inverse; the textbook gives two acceptance–rejection algorithms (Ahrens–Dieter for $0 < \alpha < 1$, Cheng for $\alpha > 1$). Scale is free: $\text{Gamma}(\alpha, \beta) = \beta \cdot \text{Gamma}(\alpha, 1)$. The $m$-Erlang algorithm of §4.3 is the integer case; for $m \gtrsim 10$ the general gamma generator is faster, and for large $m$ the product $\prod u_i$ can underflow toward $0$ before its logarithm is taken.

Lognormal: matching a target mean and variance. The parameters $\mu, \sigma^2$ of $LN(\mu, \sigma^2)$ are the mean and variance of $\ln X$, not of $X$. If the model is specified by $m' = \mathbf{E}[X]$ and $s'^2 = \operatorname{Var}(X)$, convert first:

$$\mu = \ln\frac{m'^2}{\sqrt{m'^2 + s'^2}}, \qquad \sigma^2 = \ln\!\Big(1 + \frac{s'^2}{m'^2}\Big),$$

then generate $Y \sim N(\mu, \sigma^2)$ and return $X = e^Y$. (Check with $m' = 2$, $s'^2 = 1$: simulated mean $2.00$, variance $1.00$.)

Empirical distributions. With raw observations $X_{(1)} \le \cdots \le X_{(n)}$ and the piecewise-linear empirical c.d.f. of §2, the inverse transform needs no search because the corners sit at levels $0, \frac{1}{n-1}, \frac{2}{n-1}, \dots, 1$:

1. Generate $u$; let $P = (n-1)u$ and $I = \lfloor P \rfloor + 1$.
2. Return $X = X_{(I)} + (P - I + 1)\big(X_{(I+1)} - X_{(I)}\big)$.

The output always lies in $[X_{(1)}, X_{(n)}]$ — no values beyond the observed range can ever appear, which is the known limitation of empirical input modelling.

4.4. Generating Discrete Random Variates

Inverse transform works for every discrete distribution (§4.2.1), but is not always the fastest. Besides the specific formulas below, the section ends with two table-based methods — table lookup and the alias method — that trade setup and storage for speed.

Bernoulli$(p)$. $X = 1$ with probability $p$, $0$ with probability $1-p$. Inverse transform with two values:

1. Generate $u \sim U(0,1)$.
2. If $u \le p$, return $X = 1$; otherwise return $X = 0$.

Binomial$(n, p)$. The number of successes in $n$ independent Bernoulli$(p)$ trials is a sum, so use convolution:

1. Generate $Y_1, \dots, Y_n$ i.i.d. Bernoulli$(p)$.
2. Return $X = Y_1 + \cdots + Y_n$.

Costs $n$ uniforms per variate; fine for small $n$.

Discrete uniform $DU(i, j)$. $P(X = k) = \dfrac{1}{j-i+1}$ for $k = i, i+1, \dots, j$. Generalizing the $\lceil 100u \rceil$ trick of §4.2.1:

1. Generate $u \sim U(0,1)$.
2. Return $X = i + \lfloor (j - i + 1)\,u \rfloor$.

Here $\lfloor x \rfloor$ is the largest integer $\le x$. The product $(j-i+1)u$ is uniform on $[0, j-i+1)$; its floor is uniform on $\{0, 1, \dots, j-i\}$; adding $i$ shifts to $\{i, \dots, j\}$. (The floor form is slightly safer than the ceiling form when $u$ can be exactly $0$.)

Geometric$(p)$. $p(k) = (1-p)^{k-1}p$, $k = 1, 2, \dots$ — the number of trials up to and including the first success. Instead of simulating the trials (§4.2.5), invert the c.d.f. in closed form. The inverse-transform intervals are

$$0 < u \le p \Rightarrow X = 1, \qquad p < u \le p + (1-p)p \Rightarrow X = 2, \qquad \dots$$

and in general $X = k$ when

$$\sum_{i=1}^{k-1}(1-p)^{i-1}p \;<\; u \;\le\; \sum_{i=1}^{k}(1-p)^{i-1}p .$$

Each side is a geometric series: $\sum_{i=1}^{k}(1-p)^{i-1}p = p\,\dfrac{1 - (1-p)^k}{1 - (1-p)} = 1 - (1-p)^k$. So the condition is

$$1 - (1-p)^{k-1} < u \le 1 - (1-p)^{k} \quad\Longleftrightarrow\quad (1-p)^{k} \le 1 - u < (1-p)^{k-1}.$$

Take logarithms. Since $0 < 1-p < 1$, $\ln(1-p) < 0$, so the inequalities flip:

$$k\ln(1-p) \le \ln(1-u) < (k-1)\ln(1-p),$$

and dividing by the negative number $\ln(1-p)$ flips them back:

$$k \ge \frac{\ln(1-u)}{\ln(1-p)} > k - 1 .$$

The integer $k$ sitting just above $\frac{\ln(1-u)}{\ln(1-p)}$ is its ceiling. Finally replace $1-u$ by $u$ (same distribution):

1. Generate $u \sim U(0,1)$.
2. Return $X = \left\lceil \dfrac{\ln u}{\ln(1-p)} \right\rceil$.

One uniform and one logarithm, regardless of $p$ — compare $1/p$ uniforms for the trial-by-trial method. The formula is the discrete twin of the exponential’s $-\frac{1}{\lambda}\ln u$: the geometric is the exponential’s discrete analogue (it has the same memoryless property), and the generators have the same shape.

Arbitrary discrete distribution. $X \in \{0, 1, \dots, n\}$ with $P(X = i) = p(i)$. Inverse transform:

1. Generate $u \sim U(0,1)$.
2. Return $X = i$ for the $i$ with $\displaystyle\sum_{j=0}^{i-1} p(j) \;\le\; u \;<\; \sum_{j=0}^{i} p(j)$.

This is the search of §4.2.1; for large $n$ it is slow, and the first improvement is to test the values in decreasing order of $p(i)$. The two methods below go further.

4.4.1. Table lookup

Idea. If all probabilities are given to $q$ decimal places, then each $p(i)$ is a whole number of units of size $10^{-q}$. Build a vector $m$ of length $10^q$ in which value $i$ appears $10^q\,p(i)$ times. Picking a uniformly random entry of $m$ then returns $i$ with exactly the right probability — no search at all.

Example. $P(0) = 0.15$, $P(1) = 0.20$, $P(2) = 0.37$, $P(3) = 0.28$, so $q = 2$. Construct $m$ of length $100$:

$$m = [\underbrace{0, \dots, 0}_{15},\; \underbrace{1, \dots, 1}_{20},\; \underbrace{2, \dots, 2}_{37},\; \underbrace{3, \dots, 3}_{28}] .$$

1. Generate $J \sim DU(1, 100)$, i.e. $J = 1 + \lfloor 100u \rfloor$.
2. Return $X = m_J$.

Check: $P(X = 2) = P(J \text{ lands on one of the 37 entries equal to } 2) = 37/100$. ✓

Simple and very fast (one uniform, one index), but it needs setup and $10^q$ storage locations — $q = 6$ already means a million entries.

4.4.2. Improved table lookup

Idea. Split each probability by decimal digit and keep one table per digit. Same example:

$i$	$p(i)$	tenths digit	hundredths digit
0	0.15	1	5
1	0.20	2	0
2	0.37	3	7
3	0.28	2	8
	sum	8	20

Build one table from the tenths digits (value $i$ appears “tenths digit” times) and one from the hundredths digits:

$$m^{(1)} = [0,\ 1, 1,\ 2, 2, 2,\ 3, 3] \quad (\text{length } 8), \qquad m^{(2)} = [\underbrace{0,\dots,0}_{5},\ \underbrace{2,\dots,2}_{7},\ \underbrace{3,\dots,3}_{8}] \quad (\text{length } 20).$$

Total storage $8 + 20 = 28$ instead of $100$. (Value 1 is absent from $m^{(2)}$ because its hundredths digit is $0$.)

The tenths table carries total mass $8/10 = 0.8$ and the hundredths table carries $20/100 = 0.2$. So:

1. Generate $u \sim U(0,1)$.
2. If $u \le 0.8$, generate $J_1 \sim DU(1, 8)$ and return $X = m^{(1)}_{J_1}$.
3. Otherwise generate $J_2 \sim DU(1, 20)$ and return $X = m^{(2)}_{J_2}$.

Check: $P(X = 2) = P(X = 2 \mid \text{table 1})\,P(\text{table 1}) + P(X = 2 \mid \text{table 2})\,P(\text{table 2}) = \frac{3}{8}\cdot\frac{8}{10} + \frac{7}{20}\cdot\frac{2}{10} = 0.30 + 0.07 = 0.37$. ✓ The two terms are exactly the tenths and hundredths parts of $0.37$.

This is composition (§4.2.2) in disguise: $p = 0.8\,p^{(1)} + 0.2\,p^{(2)}$, where $p^{(j)}$ is the uniform-pick distribution on table $j$. It costs two uniforms per variate instead of one, in exchange for far less storage when $q$ is large.

4.4.3. The alias method

The alias method keeps table lookup’s speed — one uniform-pick plus one comparison — with only $2(n+1)$ storage locations, whatever the precision of the probabilities.

Setup. Two vectors of length $n+1$:

• cutoff values $F_0, F_1, \dots, F_n \in [0,1]$,
• aliases $L_0, L_1, \dots, L_n \in \{0, 1, \dots, n\}$.

Algorithm.

1. Generate $I \sim DU(0, n)$ and $u \sim U(0,1)$ independently.
2. If $u \le F_I$, return $X = I$; otherwise return $X = L_I$.

How to read it. Each value $I$ is drawn with probability $\frac{1}{n+1}$ — an equal-size “column”. Inside column $I$, a fraction $F_I$ of the column belongs to $I$ itself and the remaining $1 - F_I$ belongs to one other value, its alias $L_I$. So

$$P(X = k) = \frac{1}{n+1}\Big[\,F_k + \sum_{i \,:\, L_i = k} (1 - F_i)\Big] :$$

$k$’s own share of its column, plus whatever the other columns hand to it as an alias. The setup consists of choosing $F$ and $L$ so that this equals $p(k)$ for every $k$.

It resembles acceptance–rejection (draw a candidate, compare $u$ to a threshold), with one crucial difference: on “rejection” we do not start over — we return the alias. There is no loop, so exactly one $I$ and one $u$ are consumed per variate.

Setup for the example $p = (0.15, 0.20, 0.37, 0.28)$, $n + 1 = 4$. Scale by $n+1$: $q = 4p = (0.60, 0.80, 1.48, 1.12)$. Each column has capacity $1$; values with $q_i < 1$ have spare room, values with $q_i > 1$ have excess. Repeatedly fill a “small” column from a “large” one:

step	small column $i$	gets $F_i = q_i$	alias $L_i$	donor’s remaining $q$
1	$0$ ($0.60$)	$0.60$	$2$	$q_2 = 1.48 - 0.40 = 1.08$
2	$1$ ($0.80$)	$0.80$	$2$	$q_2 = 1.08 - 0.20 = 0.88$
3	$2$ ($0.88$)	$0.88$	$3$	$q_3 = 1.12 - 0.12 = 1.00$
4	$3$ ($1.00$)	$1.00$	— (never used)

Result: $F = (0.60,\ 0.80,\ 0.88,\ 1.00)$, $L = (2,\ 2,\ 3,\ 3)$.

Check: $P(X = 2) = \frac14\big[0.88 + (1 - 0.60) + (1 - 0.80)\big] = \frac14 (1.48) = 0.37$ ✓, and $P(X = 3) = \frac14\big[1.00 + (1 - 0.88)\big] = \frac14(1.12) = 0.28$ ✓. (The table is not unique — a different pairing order gives different valid $F, L$.)

Assessment. Very fast once set up; drawbacks are the setup computation and the $2(n+1)$ storage locations. For a fixed discrete distribution sampled millions of times, it is the method of choice.

Summary of the discrete methods.

method	uniforms / variate	time per variate	storage	setup
inverse transform	1	up to $n$ comparisons	$n+1$	sort (optional)
table lookup	1	$O(1)$	$10^q$	build table
improved table lookup	2	$O(1)$	sum of digits	build tables
alias	2 (or 1 split in two)	$O(1)$	$2(n+1)$	build $F, L$

4.5. Generating Random Vectors

Why a separate section. Everything so far produces one random variable at a time, and successive variates are independent. Many inputs come in correlated groups: a job’s processing times at two consecutive stations, a customer’s order size and service time, daily demands for related products. Generating each coordinate on its own throws the correlation away. We need a way to produce a whole vector $X = (X_1, \dots, X_d)^{\mathsf T}$ with a prescribed mean and a prescribed dependence structure.

4.5.1. Mean vector and covariance matrix

For a random vector $X \in \mathbb{R}^d$:

$$\mu = \mathbf{E}[X] = \big(\mathbf{E}[X_1], \dots, \mathbf{E}[X_d]\big)^{\mathsf T}, \qquad \Sigma = \operatorname{cov}(X) = \mathbf{E}\big[(X-\mu)(X-\mu)^{\mathsf T}\big],$$

so $\Sigma$ is $d \times d$ with entries $\Sigma_{ij} = \operatorname{cov}(X_i, X_j)$ — variances on the diagonal, covariances off it (§2.3). Two facts about $\Sigma$ that the construction below relies on:

• $\Sigma$ is symmetric, since $\operatorname{cov}(X_i, X_j) = \operatorname{cov}(X_j, X_i)$.
• $\Sigma$ is positive semidefinite: for any constant vector $a$, $a^{\mathsf T}\Sigma a = \operatorname{Var}(a^{\mathsf T}X) \ge 0$. A variance cannot be negative.

Linear transformations. If $Y = b + AX$ for a constant vector $b$ and constant matrix $A$, then

$$\mathbf{E}[Y] = b + A\mu, \qquad \operatorname{cov}(Y) = A\,\Sigma\,A^{\mathsf T}.$$

(The second follows from $Y - \mathbf{E}[Y] = A(X - \mu)$, so $\operatorname{cov}(Y) = \mathbf{E}[A(X-\mu)(X-\mu)^{\mathsf T}A^{\mathsf T}] = A\Sigma A^{\mathsf T}$.) This is the whole trick: we can manufacture any covariance we like by applying a matrix to something simple.

4.5.2. The multivariate normal distribution

$X \sim N_d(\mu, \Sigma)$ is the $d$-dimensional normal with mean vector $\mu$ and covariance matrix $\Sigma$. Its one property we need: a linear transformation of a multivariate normal is again multivariate normal. Combined with the formulas above, if $X \sim N_d(\mu, \Sigma)$ then $b + AX \sim N(b + A\mu,\; A\Sigma A^{\mathsf T})$.

The simplest member of the family is the standard one. Let $Z_1, \dots, Z_d$ be i.i.d. $N(0,1)$ (§4.3.1) and $Z = (Z_1, \dots, Z_d)^{\mathsf T}$. Then $\mathbf{E}[Z] = 0$ and $\operatorname{cov}(Z) = I_d$ (variances 1, covariances 0 by independence), so $Z \sim N_d(0, I_d)$.

The plan. Take $X = \mu + AZ$. Then $X$ is multivariate normal with

$$\mathbf{E}[X] = \mu, \qquad \operatorname{cov}(X) = A\,I\,A^{\mathsf T} = AA^{\mathsf T}.$$

So all we need is a matrix $A$ with

$$AA^{\mathsf T} = \Sigma,$$

a “square root” of $\Sigma$. There are two standard ways to get one; the lecture uses the second.

1. Cholesky factorization. Any symmetric positive-definite $\Sigma$ can be written $\Sigma = LL^{\mathsf T}$ with $L$ lower-triangular. Take $A = L$. Cheap to compute ($O(d^3/3)$) and what most software does.
2. Eigenvalue (spectral) decomposition — next.

4.5.3. Eigenvalue decomposition

Definitions. For a $d \times d$ matrix $\Sigma$, if a constant $\lambda$ and a non-zero vector $v$ satisfy

$$\Sigma v = \lambda v,$$

then $\lambda$ is an eigenvalue of $\Sigma$ and $v$ a corresponding eigenvector. (Geometrically: $\Sigma$ merely stretches $v$ by the factor $\lambda$ without turning it.)

Facts for a real symmetric matrix (spectral theorem):

• It has $d$ linearly independent eigenvectors, which can be chosen orthogonal and of unit length: $v_1, \dots, v_d$ with $v_i^{\mathsf T}v_j = 0$ for $i \ne j$ and $\|v_i\| = 1$.
• Put them side by side in $P = [v_1 \; v_2 \; \cdots \; v_d]$. Orthonormal columns mean $P^{\mathsf T}P = I$, i.e. $P^{-1} = P^{\mathsf T}$ — the inverse is free.
• With $\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_d)$,

$$\Sigma = P\,\Lambda\,P^{\mathsf T}.$$

• Because $\Sigma$ is positive semidefinite, every $\lambda_i \ge 0$ ($\lambda_i = v_i^{\mathsf T}\Sigma v_i \ge 0$), so $\Lambda^{1/2} = \operatorname{diag}(\sqrt{\lambda_1}, \dots, \sqrt{\lambda_d})$ is a real matrix with $\Lambda^{1/2}\Lambda^{1/2} = \Lambda$.

The square root. Take $A = P\Lambda^{1/2}$. Then

$$AA^{\mathsf T} = P\Lambda^{1/2}\,(P\Lambda^{1/2})^{\mathsf T} = P\Lambda^{1/2}\Lambda^{1/2}P^{\mathsf T} = P\Lambda P^{\mathsf T} = \Sigma. \; \checkmark$$

(We used $(\Lambda^{1/2})^{\mathsf T} = \Lambda^{1/2}$ since it is diagonal.)

Verification, step by step. Let $Z = (Z_1, \dots, Z_d)^{\mathsf T}$ with $Z_i$ i.i.d. $N(0,1)$ and define

$$X = \mu + P\Lambda^{1/2}Z .$$

Mean: $\mathbf{E}[X] = \mu + P\Lambda^{1/2}\,\mathbf{E}[Z] = \mu + 0 = \mu$.

Covariance:

$$\begin{aligned} \operatorname{cov}(X) &= \mathbf{E}\big[(X-\mu)(X-\mu)^{\mathsf T}\big] \\ &= \mathbf{E}\big[P\Lambda^{1/2}Z\,Z^{\mathsf T}\Lambda^{1/2}P^{\mathsf T}\big] && (X-\mu = P\Lambda^{1/2}Z;\ (ABC)^{\mathsf T} = C^{\mathsf T}B^{\mathsf T}A^{\mathsf T}) \\ &= P\Lambda^{1/2}\,\mathbf{E}[ZZ^{\mathsf T}]\,\Lambda^{1/2}P^{\mathsf T} && (\text{constants come out of the expectation}) \\ &= P\Lambda^{1/2}\,I\,\Lambda^{1/2}P^{\mathsf T} && (\mathbf{E}[Z_iZ_j] = 1 \text{ if } i=j,\ 0 \text{ otherwise}) \\ &= P\Lambda P^{\mathsf T} = \Sigma . \end{aligned}$$

And $X$ is normal because each coordinate is a linear combination of the independent normals $Z_1, \dots, Z_d$ plus a constant.

Algorithm (multivariate normal via eigen decomposition).

Setup (once): compute the eigenvalues $\lambda_i$ and orthonormal eigenvectors $v_i$ of $\Sigma$; form $P$ and $\Lambda^{1/2}$; store $A = P\Lambda^{1/2}$.

1. Generate $Z_1, \dots, Z_d$ i.i.d. $N(0,1)$ (Box–Muller or polar method, §4.3).
2. Return $X = \mu + P\Lambda^{1/2}Z$.

One matrix–vector multiplication per vector; the expensive part (the decomposition) is paid once.

Example ($d = 2$). Two correlated standard normals with correlation $\rho$:

$$\mu = \begin{pmatrix}\mu_1\\ \mu_2\end{pmatrix}, \qquad \Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1\end{pmatrix}.$$

Eigenvalues and eigenvectors (check by multiplying out):

$$\Sigma\begin{pmatrix}1\\1\end{pmatrix} = (1+\rho)\begin{pmatrix}1\\1\end{pmatrix}, \qquad \Sigma\begin{pmatrix}1\\-1\end{pmatrix} = (1-\rho)\begin{pmatrix}1\\-1\end{pmatrix},$$

so $\lambda_1 = 1+\rho$, $\lambda_2 = 1-\rho$ (both $\ge 0$ since $|\rho| \le 1$), and after normalizing to unit length,

$$P = \frac{1}{\sqrt2}\begin{pmatrix} 1 & 1 \\ 1 & -1\end{pmatrix}, \qquad \Lambda^{1/2} = \begin{pmatrix}\sqrt{1+\rho} & 0 \\ 0 & \sqrt{1-\rho}\end{pmatrix}, \qquad A = P\Lambda^{1/2} = \frac{1}{\sqrt2}\begin{pmatrix} \sqrt{1+\rho} & \sqrt{1-\rho} \\ \sqrt{1+\rho} & -\sqrt{1-\rho}\end{pmatrix}.$$

The algorithm therefore returns

$$X_1 = \mu_1 + \frac{\sqrt{1+\rho}\,Z_1 + \sqrt{1-\rho}\,Z_2}{\sqrt2}, \qquad X_2 = \mu_2 + \frac{\sqrt{1+\rho}\,Z_1 - \sqrt{1-\rho}\,Z_2}{\sqrt2}.$$

Check: $\operatorname{Var}(X_1) = \frac{(1+\rho) + (1-\rho)}{2} = 1$, and $\operatorname{cov}(X_1, X_2) = \frac{(1+\rho) - (1-\rho)}{2} = \rho$. ✓

The same $\Sigma$ by Cholesky. $L = \begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1-\rho^2}\end{pmatrix}$ satisfies $LL^{\mathsf T} = \Sigma$, giving

$$X_1 = \mu_1 + Z_1, \qquad X_2 = \mu_2 + \rho Z_1 + \sqrt{1-\rho^2}\,Z_2 .$$

Different formulas, same distribution — both are valid square roots of $\Sigma$. (For general variances $\sigma_1^2, \sigma_2^2$: $X_1 = \mu_1 + \sigma_1 Z_1$, $X_2 = \mu_2 + \sigma_2(\rho Z_1 + \sqrt{1-\rho^2}\,Z_2)$.) The Cholesky form shows the mechanism plainly: $X_2$ borrows a share $\rho$ of $X_1$’s noise and tops up with fresh noise so that its variance is still 1.

Remarks.

• The $Z_i$ must be independent — this is where the quality of the normal generator (§4.3.1’s spiral warning) matters doubly.
• The construction only controls the first two moments. For a non-normal vector with prescribed marginals and correlation, transforming $Z$ by a matrix destroys the marginals; one then uses $Z$ only to carry the dependence and maps each coordinate through $F_i^{-1}(\Phi(\cdot))$ (the “NORTA” idea). Beyond the scope of this course.

4.6. Generating Stochastic Processes

A stochastic process (§2.5) is a whole family of random variables indexed by time. Generating one means generating a sample path. The workhorse is the Poisson process — the arrival process of almost every queueing model — and its time-varying generalization.

4.6.1. The Poisson process

Counting process. Let $S_n$ be the time of the $n$-th event (the $n$-th arrival time), with $S_0 = 0$, and let

$$X_n = S_n - S_{n-1}$$

be the $n$-th interarrival time, so $S_n = \sum_{i=1}^n X_i$. Define

$$N(t) = \max\{n : S_n \le t\} = \text{number of events that have occurred up to time } t .$$

$\{N(t) : t \ge 0\}$ is a counting process: it starts at $0$, is integer-valued, non-decreasing, and jumps by $1$ at each $S_n$.

Definition. If the interarrival times $X_1, X_2, \dots$ are i.i.d. $\text{Exp}(\lambda)$, then $\{N(t)\}$ is called a Poisson process with rate $\lambda$, written $\text{PP}(\lambda)$.

Properties of $\text{PP}(\lambda)$.

1. Poisson counts. For every $t$,

   $$P(N(t) = k) = e^{-\lambda t}\frac{(\lambda t)^k}{k!}, \qquad k = 0, 1, 2, \dots,$$

   i.e. $N(t) \sim \text{Poisson}(\lambda t)$, so $\mathbf{E}[N(t)] = \lambda t$. Hence $\lambda = \mathbf{E}[N(t)]/t$: the rate (or arrival rate) is the expected number of events per unit time. More generally the number of events in any interval of length $s$ is $\text{Poisson}(\lambda s)$, and counts in disjoint intervals are independent. (This is the Erlang formula of §4.2.3 read the other way: $P(N(t) \ge n) = P(S_n \le t)$.)

2. Arrival-time theorem. Given that $N(t) = k$, the $k$ arrival times $(S_1, \dots, S_k)$ are distributed as the order statistics of $k$ i.i.d. $U(0, t)$ random variables. That is: take $U_1, \dots, U_k$ i.i.d. $U(0,t)$, sort them into ascending order $\bar U_{(1)} < \bar U_{(2)} < \cdots < \bar U_{(k)}$; then $(S_1, \dots, S_k)$ has the same joint distribution as $(\bar U_{(1)}, \dots, \bar U_{(k)})$. In words: given how many events occurred, they are scattered uniformly and independently over the interval — the Poisson process has no preferred times.

Bernoulli splitting. Let $\{N(t)\}$ be $\text{PP}(\lambda)$ and fix probabilities $p_1, \dots, p_r > 0$ with $\sum_{j=1}^r p_j = 1$. Assign each arrival, independently of everything else, to category $i$ with probability $p_i$, and let $N_i(t)$ count the category-$i$ arrivals, so $N(t) = N_1(t) + \cdots + N_r(t)$.

Theorem. Under Bernoulli splitting the $r$ processes $\{N_1(t)\}, \dots, \{N_r(t)\}$ are independent Poisson processes, with $\{N_i(t)\} \sim \text{PP}(\lambda p_i)$.

The case $r = 2$, with probabilities $p$ and $1-p$, is the one we will use: the arrivals that pass a coin toss form $\text{PP}(\lambda p)$, those that fail form $\text{PP}(\lambda(1-p))$, and the two streams are independent of each other. (The independence is the surprising part: knowing many arrivals went to stream 1 tells you nothing about stream 2.)

4.6.2. The nonhomogeneous Poisson process

Motivation: traffic on a highway. A $\text{PP}(\lambda)$ has one constant rate for all time. Real arrival rates vary:

A single $\lambda$ would average the rush hours and the quiet hours together and get both wrong. We need a process whose rate is a function of time.

Definition. A counting process $\{N(t)\}$ is a nonhomogeneous Poisson process (NPP) with rate function $\lambda(t) \ge 0$ if

1. $N(0) = 0$;
2. $\{N(t)\}$ has independent increments: counts in disjoint time intervals are independent;
3. in a short window $(t, t+h]$,

   $$\begin{aligned} P\big(N(t+h) = k \mid N(t) = k\big) &= 1 - \lambda(t)h + o(h) && \text{(no event)} \\ P\big(N(t+h) = k+1 \mid N(t) = k\big) &= \lambda(t)h + o(h) && \text{(exactly one event)} \\ P\big(N(t+h) = k+j \mid N(t) = k\big) &= o(h), \quad j \ge 2 && \text{(two or more events)} \end{aligned}$$

   Here $o(h)$ means a quantity that vanishes faster than $h$ as $h \to 0$. In words: over a very short interval almost surely either $0$ or $1$ event happens, and the chance of $1$ event is proportional to the length of the interval with proportionality constant $\lambda(t)$ — the current value of the rate.

Counts. Define the mean-value function $m(t) = \int_0^t \lambda(s)\,ds$. Then

$$P(N(t) = k) = e^{-m(t)}\frac{m(t)^k}{k!}, \qquad k = 0, 1, 2, \dots,$$

so $N(t) \sim \text{Poisson}\big(m(t)\big)$ and $\mathbf{E}[N(t)] = \int_0^t \lambda(s)\,ds$ — the area under the rate curve. More generally $N(t+s) - N(t) \sim \text{Poisson}\big(\int_t^{t+s}\lambda(u)\,du\big)$. With a constant rate $\lambda(t) \equiv \lambda$ we get $m(t) = \lambda t$ and recover $\text{PP}(\lambda)$: the NPP is the generalization.

Nonhomogeneous Bernoulli splitting. Now let the coin be biased by time: fix a function $p(t) : [0, \infty) \to [0,1]$ and assign an arrival occurring at time $t$ to stream 1 with probability $p(t)$, to stream 2 with probability $1 - p(t)$.

Theorem. If $\{N(t)\} \sim \text{PP}(\lambda)$ is split with time-dependent probability $p(t)$, then $\{N_1(t)\}$ and $\{N_2(t)\}$ are independent nonhomogeneous Poisson processes with rate functions $\lambda\,p(t)$ and $\lambda\,(1 - p(t))$:

$$P(N_1(t) = k) = \exp\!\Big(-\lambda\!\int_0^t p(s)\,ds\Big)\frac{\big[\lambda\int_0^t p(s)\,ds\big]^k}{k!}, \qquad P(N_2(t) = k) = \exp\!\Big(-\lambda\!\int_0^t (1-p(s))\,ds\Big)\frac{\big[\lambda\int_0^t (1-p(s))\,ds\big]^k}{k!}.$$

(The same holds if the parent process is itself an NPP with rate $\lambda(t)$: the pieces are NPPs with rates $\lambda(t)p(t)$ and $\lambda(t)(1-p(t))$.) Thinning out a Poisson stream with a time-varying coin produces exactly an NPP — this is the key to generating one.

4.6.3. Generating a Poisson process on $[0, T]$

Approach 1: lay down the interarrival times. Since the gaps are i.i.d. $\text{Exp}(\lambda)$, generate them one at a time by inverse transform and accumulate. Set $s_0 = 0$, $i = 0$.

1. Generate $u \sim U(0,1)$.
2. Set $s_{i+1} = s_i - \dfrac{1}{\lambda}\ln u$ $\quad$(previous arrival plus an $\text{Exp}(\lambda)$ gap).
3. Set $i = i + 1$.
4. If $s_i \le T$, go back to step 1; otherwise stop. The arrival times are $s_1, \dots, s_{i-1}$ and $N(T) = i - 1$.

The last generated time, the one that overshoots $T$, is discarded; it only serves to tell us we are done. This is the natural method inside a discrete-event simulation, where the next arrival is scheduled as soon as the current one occurs.

Approach 2: use the arrival-time theorem. Decide the count first, then place the points.

1. Generate $N(T) \sim \text{Poisson}(\lambda T)$ (a discrete variate — inverse transform with the p.m.f. computed recursively, $p(k+1) = p(k)\cdot\lambda T/(k+1)$).
2. Generate $u_1, \dots, u_{N(T)}$ i.i.d. $U(0, T)$.
3. Sort them; return $\bar u_{(1)} < \bar u_{(2)} < \cdots < \bar u_{(N(T))}$ as the arrival times $s_1, \dots, s_{N(T)}$.

By property 2 the output has exactly the distribution of a $\text{PP}(\lambda)$ path on $[0,T]$. This approach generates the whole path at once (convenient for vectorized code) at the cost of a sort.

4.6.4. Generating a nonhomogeneous Poisson process

A tempting but wrong method. Mimic Approach 1 with the current rate: having generated $t_{i-1}$, draw the next gap from $\text{Exp}(\lambda(t_{i-1}))$,

$$t_i = t_{i-1} - \frac{1}{\lambda(t_{i-1})}\ln u .$$

This uses the rate at the previous arrival for the entire gap, as if the rate stayed frozen until the next event. It does not: the rate may rise steeply right after $t_{i-1}$.

If $\lambda(t_{i-1})$ happens to be small, the gap has mean $1/\lambda(t_{i-1})$, which is large, and the process can jump clean over a high-rate period in which many arrivals should have occurred. In a test with a rate function that dips to $0.25$ and then peaks at $3.25$, this method produced on average $3.9$ arrivals on $[0,10]$ where the true mean is $7.3$, and only $1.5$ arrivals in the rush-hour window where there should be $4.9$. So this is not a correct way.

The correct way: the thinning method. Generate too many candidate arrivals from a homogeneous Poisson process whose rate is at least as large as $\lambda(t)$ everywhere, then keep each candidate at time $t$ with probability $\lambda(t)/\lambda^*$. Let

$$\lambda^* = \max_{0 \le t \le T}\lambda(t).$$

Algorithm (thinning), to generate $t_i$ given $t_{i-1}$:

1. Set $t = t_{i-1}$.
2. Generate $u_1 \perp u_2 \sim U(0,1)$.
3. Set $t = t - \dfrac{1}{\lambda^*}\ln u_1$ $\quad$(advance by an $\text{Exp}(\lambda^*)$ gap: the next candidate).
4. If $u_2 \le \dfrac{\lambda(t)}{\lambda^*}$, accept: set $t_i = t$. Otherwise go back to step 2 and try again (keeping the advanced $t$).

Repeat until $t > T$.

How to understand the algorithm. It is nonhomogeneous Bernoulli splitting run deliberately. The candidate times produced by steps 2–3 form a homogeneous $\text{PP}(\lambda^*)$ — that is just Approach 1 with rate $\lambda^*$. Step 4 then assigns each candidate at time $t$ to one of two streams with a time-dependent coin:

$$N_1(t) = \text{number of candidates accepted by time } t, \qquad N_2(t) = \text{number rejected}, \qquad p(t) = \frac{\lambda(t)}{\lambda^*}.$$

By the splitting theorem, $\{N_1(t)\}$ is a nonhomogeneous Poisson process with rate function

$$\lambda^*\cdot p(t) = \lambda^*\cdot\frac{\lambda(t)}{\lambda^*} = \lambda(t),$$

exactly the target. (The rejected stream is an NPP with rate $\lambda^* - \lambda(t)$, independent of the accepted one — we simply throw it away.) Note that $p(t) \le 1$ precisely because $\lambda^*$ was chosen as the maximum of $\lambda$; that is the only role $\lambda^*$ plays. The analogy with acceptance–rejection (§4.2.4) is exact: $\lambda^*$ is a flat majorizing function for $\lambda(t)$, and $\lambda(t)/\lambda^*$ is the acceptance probability.

Checks. In the same test as above, thinning gave $7.29$ arrivals on average against the true $7.32$, and $4.89$ in the rush-hour window against $4.88$.

Efficiency. The expected fraction of candidates accepted is

$$\frac{\int_0^T \lambda(t)\,dt}{\lambda^* T} = \frac{\text{area under } \lambda}{\text{area of the rectangle}},$$

the same picture as the AR efficiency $1/c$. If $\lambda(t)$ has a sharp peak, $\lambda^*$ is large and most candidates are rejected during the quiet periods. The standard fix is the AR fix: use a piecewise constant $\lambda^*$ — split $[0,T]$ into subintervals, each with its own local maximum — so the majorizing rectangle hugs the curve.

4.7. Example - generating the inputs of the inventory system

The inventory model of §1.3 needs three kinds of random input. Each one is generated from $U(0,1)$ random numbers with a method of this chapter.

input	distribution	method	formula
time between demands	$\text{Exp}$ with mean $0.1$ month	inverse transform (§4.2.1)	$X = -0.1\ln u$
demand size	$D \in \{1,2,3,4\}$ w.p. $\frac16, \frac13, \frac13, \frac16$	discrete inverse transform (§4.2.1)	compare $u$ with $\frac16, \frac12, \frac56$
order lead time	$U(0.5, 1)$	inverse transform (§4.3)	$L = 0.5 + 0.5\,u$

4.7.1. Time between demands

The mean is $1/\lambda = 0.1$ month, so $\lambda = 10$ per month and $F(x) = 1 - e^{-10x}$, $x \ge 0$. Set $u = F(x)$ and solve:

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

Since $1 - u \sim U(0,1)$ as well:

1. Generate $u \sim U(0,1)$.
2. Return $X = -0.1\ln u$.

(If the RNG can return $u = 0$, use $-0.1\ln(1-u)$ instead; see §4.2.1.)

4.7.2. Demand size

The c.d.f. is a staircase:

$d$	1	2	3	4
$p(d)$	$1/6$	$1/3$	$1/3$	$1/6$
$F(d)$	$1/6$	$1/2$	$5/6$	$1$

Discrete inverse transform: return the smallest $d$ with $u \le F(d)$.

1. Generate $u \sim U(0,1)$.
2. If $u \le 1/6$, return $D = 1$; else if $u \le 1/2$, return $D = 2$; else if $u \le 5/6$, return $D = 3$; else return $D = 4$.

The cut points $1/6, 1/2, 5/6$ split $(0,1]$ into four intervals whose widths are the probabilities, and $u$ falls into exactly one of them:

Alternative: table lookup (§4.4.1). All the probabilities are multiples of $1/6$, so store

$$m = [1, 2, 2, 3, 3, 4]$$

(value $d$ appears $6\,p(d)$ times), generate $J = 1 + \lfloor 6u \rfloor$, and return $D = m_J$. Check: $P(D = 2) = P(J \in \{2, 3\}) = 2/6 = 1/3$. Because $m$ is sorted, this returns the same $D$ as the search for every $u$ (up to the boundary points), just without any comparisons.

4.7.3. Order lead time

For $U(0.5, 1)$, $F(x) = \dfrac{x - 0.5}{0.5}$ on $[0.5, 1]$. Set $u = F(x)$:

$$u = \frac{x - 0.5}{0.5} \;\Longrightarrow\; x = 0.5 + 0.5\,u ,$$

which is the formula $a + (b - a)u$ of §4.3.

1. Generate $u \sim U(0,1)$.
2. Return $L = 0.5 + 0.5\,u$.

4.7.4. Worked example

Suppose the generator gives the following random numbers:

input	$u$	variate
time to demand 1	0.81	$-0.1\ln 0.81 = 0.021$ month
size of demand 1	0.55	$\tfrac12 < 0.55 \le \tfrac56 \Rightarrow D = 3$
time to demand 2	0.17	$-0.1\ln 0.17 = 0.177$ month
size of demand 2	0.12	$0.12 \le \tfrac16 \Rightarrow D = 1$
lead time	0.30	$0.5 + 0.5 \cdot 0.30 = 0.65$ month

Starting from $I(0) = 60$: the first demand comes at $t = 0.021$ and takes 3 items ($I = 57$); the second comes at $t = 0.021 + 0.177 = 0.198$ and takes 1 item ($I = 56$). If an order is placed at the review at $t = 1$, with this lead time it arrives at $t = 1.65$.

4.7.5. When each input is generated during the simulation

• Demand at time $t$: generate the size $D$ and set $I \leftarrow I - D$; generate the next interdemand time $X$ and schedule the next demand at $t + X$.
• Inventory evaluation (beginning of each month): if $I < s$, order $Z = S - I$, generate the lead time $L$, and schedule the order arrival at $t + L$. If $I \ge s$, nothing is generated.
• Order arrival and end of simulation: no random input.
• At $t = 0$, generate the first interdemand time to schedule the first demand.

4.7.6. Use a separate stream for each input

Following §3.3 (practical issue 3), take interdemand times, demand sizes and lead times from three different streams. The reason is concrete here: how many lead times are consumed depends on the policy $(s, S)$, since a policy that orders more often uses more of them. With one shared stream, those extra draws would shift every later demand, and two policies would face different demand sequences. With separate streams, every $(s, S)$ pair sees exactly the same demands (common random numbers, §4.2.1 advantage 4), so the difference in average cost reflects the policies rather than luck.