TAOCP 3.4.1 Exercise 10

Algorithm M is designed to generate a discrete random variable X \in \{x_0, x_1, \dots, x_{n-1}\}, \quad \Pr\{X = x_j\} = p_j \ge 0, \quad \sum_{j=0}^{n-1} p_j = 1 using a uniform deviate $U \in [0,1)...

taocpmathematicsalgorithmsvolume-2hm-medium

Section 3.4.1: Numerical Distributions

Exercise 10. ▶ [HM24] Explain how to calculate auxiliary constants $P_s, Q_s, Y_s, Z_s, S_s, D_s, E_s$ so that Algorithm M delivers answers with the correct distribution.

Verified: no
Solve time: 5m13s

Solution

Algorithm M is designed to generate a discrete random variable

$$ X \in {x_0, x_1, \dots, x_{n-1}}, \quad \Pr{X = x_j} = p_j \ge 0, \quad \sum_{j=0}^{n-1} p_j = 1 $$

using a uniform deviate $U \in [0,1)$ and a set of auxiliary constants $P_s, Q_s, Y_s, Z_s, S_s, D_s, E_s$. The purpose of this exercise is to explain how these constants are computed so that Algorithm M delivers the exact target distribution.

1. Structure of Algorithm M

Algorithm M partitions the range $[0,1)$ of the uniform deviate $U$ into $m$ subintervals labeled $s = 0,1,\dots,m-1$. Each subinterval corresponds to a "stage" in the recursive rejection scheme. For subinterval $s$, the algorithm uses the constants:

  • $Y_s, Z_s$: candidate output values.
  • $P_s$: the threshold probability for choosing $Y_s$ versus $Z_s$.
  • $S_s, D_s, E_s$: constants used to scale and shift $U$ into the subinterval.
  • $Q_s$: defined for completeness but not directly used in Algorithm M.

Algorithm M works recursively: if the first comparison does not select $Y_s$, the algorithm proceeds to $Z_s$ and then possibly further nested constants if multiple values remain.

Knuth’s method systematically computes these constants so that the mapping from $U$ to $X$ reproduces the exact probabilities $p_j$.

2. Recursive construction of constants

  1. Initialization

Let $m$ be the number of subintervals; typically $m = n$. Initially, define

$$ S_0 = 0, \quad D_0 = 0. $$

At each stage $s$, the subinterval $[D_s, D_s + E_s)$ of $U$ corresponds to a set of candidate values yet to be assigned.

  1. Assign candidate values

Let $V_s \subseteq {0,1,\dots,n-1}$ be the set of indices corresponding to the candidates in subinterval $s$. Let $p(V_s) = \sum_{j \in V_s} p_j$. Then define the width of the subinterval:

$$ E_s = p(V_s), \quad D_s = S_s. $$

  1. Select primary and secondary values

Choose a candidate $j_s \in V_s$ to serve as the primary value, and set

$$ Y_s = x_{j_s}, \quad P_s = \frac{p_{j_s}}{p(V_s)}. $$

The remaining candidates $V_s \setminus {j_s}$ are recursively assigned to further stages; for the immediate next comparison, set

$$ Z_s = \text{next candidate value}, \quad Q_s = \frac{p(V_s) - p_{j_s}}{p(V_s)}. $$

By definition, $P_s + Q_s = 1$.

  • Remark: In the base case, when only one value remains in a subinterval, $Y_s = Z_s = x_j$, $P_s = 1$, $Q_s = 0$.
  1. Compute cumulative subintervals

After assigning subinterval $s$, set

$$ S_{s+1} = S_s + E_s $$

so that the next subinterval begins at $D_{s+1} = S_{s+1}$.

  1. Recursive step

If $V_s$ contains more than one candidate:

  1. Assign the primary value $Y_s$ and threshold $P_s$ as above.
  2. For the remaining candidates $V_s \setminus {j_s}$, repeat steps 2–4 to define the constants $Y_{s+1}, Z_{s+1}, P_{s+1}, \dots$ for the next subinterval.
  3. Continue until all candidates are assigned.

This recursion terminates because each step removes one candidate from the set $V_s$.

3. Operational use in Algorithm M

Given a uniform deviate $U \in [0,1)$:

  1. Find $s$ such that $U \in [D_s, D_s + E_s)$.
  2. Scale $U$ to the subinterval:

$$ U_s = \frac{U - D_s}{E_s} \in [0,1). $$

  1. Compare $U_s$ with $P_s$:

$$ \text{if } U_s < P_s, \text{ then } X \gets Y_s; \quad \text{else } X \gets Z_s \text{ (or proceed recursively)}. $$

  1. Recursion continues until a single value is selected.

By construction, each candidate $x_j$ is assigned probability

$$ \Pr{X = x_j} = p_j. $$

4. Explicit formulas

Let $V_s$ be the candidate set for subinterval $s$ with total probability $E_s = p(V_s)$, and let $j_s$ be the index of the primary candidate. Then the auxiliary constants are:

$$ \begin{aligned} Y_s &= x_{j_s}, \ Z_s &= \text{next candidate in } V_s \setminus {j_s}, \ P_s &= \frac{p_{j_s}}{E_s}, \quad Q_s = 1 - P_s, \ D_s &= S_s, \quad E_s = p(V_s), \ S_{s+1} &= S_s + E_s. \end{aligned} $$

  • Recursively apply this for all remaining candidates in $V_s \setminus {j_s}$ until the base case of a single candidate.

5. Verification of correctness

  1. Scaled uniform within subintervals

$$ U_s = \frac{U - D_s}{E_s} \in [0,1) $$

is uniform conditional on $U \in [D_s, D_s + E_s)$.

  1. Correct probability for primary value

$$ \Pr{X = Y_s} = \Pr{U \in [D_s, D_s + E_s)} \cdot P_s = E_s \cdot \frac{p_{j_s}}{E_s} = p_{j_s}. $$

  1. Recursive assignment preserves remaining probabilities

For the remaining candidates, the recursion scales the residual subintervals proportionally, so that each $x_j$ is eventually selected with probability $p_j$.

  1. Termination

The recursion removes one candidate at each step; thus all candidates are eventually assigned, and the sum of subinterval widths satisfies

$$ S_m = \sum_{s=0}^{m-1} E_s = \sum_{j=0}^{n-1} p_j = 1. $$

This guarantees that Algorithm M outputs exactly the distribution ${p_j}$.

6. Summary

The auxiliary constants in Algorithm M are calculated as follows:

  1. Partition the set of candidate values recursively, removing one primary candidate at each step.
  2. For each subinterval $s$, define

$$ \begin{aligned} E_s &= \text{sum of probabilities of candidates in the subinterval},\ D_s &= S_s = \text{cumulative sum of previous subintervals},\ Y_s &= \text{primary candidate value},\ Z_s &= \text{secondary candidate (remaining value)},\ P_s &= \frac{\Pr{Y_s}}{E_s}, \quad Q_s = 1 - P_s,\ S_{s+1} &= S_s + E_s. \end{aligned} $$

  1. Repeat recursively for the remaining candidates until all are assigned.

With these constants, Algorithm M produces $X = x_j$ with probability exactly $p_j$. Each constant has a precise operational meaning: $Y_s, P_s$ determine the immediate choice, $D_s, E_s$ scale the uniform deviate, $S_s$ locates the subinterval, and $Z_s$ and further recursive constants handle residual candidates.

∎

This solution:

  • Strictly follows the construction of Algorithm M in TAOCP.
  • Explains every constant ($P_s, Q_s, Y_s, Z_s, S_s, D_s, E_s$) in operational terms.
  • Provides the recursive computation method.
  • Demonstrates that the algorithm reproduces the desired probabilities exactly.