Week 2: Random Variable Generation

STAT 536 — Statistical Computing

Soumojit Das

Fall 2026

The Fundamental Challenge

Generating from an Arbitrary Distribution

The Problem. How do we generate samples from arbitrary distributions?

Standing assumption

We can generate \(U \sim \Uniform(0,1)\) using pseudo-random number generators.

Our Toolkit

  1. Inverse Transform: Use \(F^{-1}\) when available
  2. Transformation: Exploit relationships between distributions
  3. Accept-Reject: General method for complex distributions

The first two are this week. Accept-reject is next week.

The Inverse Transform

Inverse Transform Method: The Foundation

Theorem (Probability Integral Transform)

Let \(F\) be any CDF and let \(F^{-}(u) = \inf\{x : F(x) \geq u\}\). Then:

  1. If \(U \sim \Uniform(0,1)\), then \(F^{-}(U) \sim F\)
  2. If \(X \sim F\) and \(F\) is continuous, then \(F(X) \sim \Uniform(0,1)\)

Essentially, this says that the CDF is a two-way bridge between uniform draws and target draws: part (1) turns a uniform into a draw from any distribution by way of its generalized inverse, and part (2) reverses the process when \(F\) is continuous.

The Generalized Inverse

The theorem needs no inverse function in the usual sense. Define

\[F^{-}(u) = \inf\{x : F(x) \geq u\}, \qquad 0 < u < 1 .\]

Right-continuity of \(F\) gives the equivalence that drives everything below:

\[F^{-}(u) \leq x \iff u \leq F(x) .\]

We can think of this as saying: asking whether the generated value is at most \(x\) is the same as asking whether the uniform draw is at most \(F(x)\), which is why the proof of part (1) is one line.

When \(F\) is continuous and strictly increasing, \(F^{-} = F^{-1}\). Otherwise \(F^{-}\) is what the algorithm actually computes, and the discrete case is the example.

Inverse Transform: Proof of (1)

For any \(x \in \mathbb{R}\), the equivalence of the previous slide gives

\[P[F^{-}(U) \leq x] = P[U \leq F(x)] \quad (\text{right-continuity of } F)\]

\[P[U \leq F(x)] = F(x) \quad (\text{since } U \sim \Uniform(0,1))\]

Therefore \(F^{-}(U)\) has CDF \(F\).

Inverse Transform: Exponential Example

Generate \(X \sim \Exp(\lambda)\).

  1. CDF: \(F(x) = 1 - e^{-\lambda x}\) for \(x \geq 0\)

  2. Set \(u = F(x) = 1 - e^{-\lambda x}\)

  3. Solve for \(x\):

    \[ \begin{aligned} u &= 1 - e^{-\lambda x}\\ 1 - u &= e^{-\lambda x}\\ \log(1-u) &= -\lambda x\\ x &= -\frac{1}{\lambda}\log(1-u) \end{aligned} \]

Simplifying with \(1 - U\)

  1. Since \(1-U \sim \Uniform(0,1)\) when \(U \sim \Uniform(0,1)\):

    \[X = -\frac{1}{\lambda}\log(U) \sim \Exp(\lambda)\]

Reading a Draw Off the CDF

Discrete Distributions via Inverse Transform

Given: \(P[X = x_i] = p_i\) for \(x_1 < x_2 < \cdots < x_k\), with \(\sum p_i = 1\).

Method 1: Linear Search

  1. Generate \(U \sim \Uniform(0,1)\)
  2. Set \(F_0 = 0\)
  3. For \(i = 1, 2, \ldots, k\):
    • \(F_i = F_{i-1} + p_i\)
    • If \(U \leq F_i\), return \(X = x_i\)

Complexity: \(O(k)\) per sample

Method 2: Binary Search

  1. Precompute \(F_i = \sum_{j=1}^i p_j\)
  2. Generate \(U \sim \Uniform(0,1)\)
  3. Binary search for \(i\) such that \(F_{i-1} < U \leq F_i\)
  4. Return \(X = x_i\)

Complexity: \(O(\log k)\) per sample

The indexing rule \(F_{i-1} < U \leq F_i\) is exactly \(X = F^{-}(U)\).

Six Distributions You Can Invert in Closed Form

When Inverse Transform Fails

Limitations.

  • No closed form for \(F^{-1}\) (e.g., Normal, Gamma)
  • \(F^{-1}\) exists but is expensive to compute
  • Numerical inversion introduces errors

The Normal Has No Closed-Form Inverse

\[F(x) = \Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-t^2/2}\, dt\]

  • No closed form for \(\Phi\) or \(\Phi^{-1}\)
  • Numerical methods exist but are complex
  • Better methods available (Box–Muller)

Solution. Use transformation methods or accept-reject algorithms.

Three Routes to the Same \(\Exp(1)\) Sample

The third panel uses accept-reject, which is next week.

Transformation Methods

Transformation Methods: Known Relationships

Key idea. Use known relationships between distributions to generate complex from simple.

Theorem (General Transformation). If \(X\) has PDF \(f_X\) and \(Y = g(X)\) where \(g\) is strictly monotone and continuously differentiable with \(g' \neq 0\) on the support of \(X\), then

\[f_Y(y) = f_X(g^{-1}(y)) \left|\frac{d}{dy}g^{-1}(y)\right| .\]

Essentially, this says that the density of a transformed variable is the original density evaluated at the pre-image, scaled by how much \(g^{-1}\) stretches or compresses locally.

Simple Transformations

  • If \(X \sim \Normal(0,1)\), then \(\mu + \sigma X \sim \Normal(\mu, \sigma^2)\)
  • If \(X \sim \Uniform(0,1)\), then \(a + (b-a)X \sim \Uniform(a,b)\)
  • If \(X \sim \Exp(1)\), then \(X/\lambda \sim \Exp(\lambda)\)

Special Relationships: Sums and Products

Theorem (Sum of Exponentials \(\to\) Gamma). If \(X_1, \ldots, X_\alpha \iid \Exp(\beta)\), then

\[Y = \sum_{i=1}^\alpha X_i \sim \Gam(\alpha, \beta)\]

when \(\alpha\) is a positive integer (rate parameterization throughout: \(\E[Y] = \alpha/\beta\)).

We can think of this as saying that a \(\Gam(\alpha, \beta)\) draw, when \(\alpha\) is a positive integer, is the sum of \(\alpha\) consecutive independent exponential waiting times: the time to the \(\alpha\)-th arrival of a rate-\(\beta\) Poisson process.

Algorithm for \(\Gam(\alpha, \beta)\), \(\alpha \in \mathbb{N}\).

  1. Generate \(U_1, \ldots, U_\alpha \iid \Uniform(0,1)\)
  2. Return \(Y = -\frac{1}{\beta}\sum_{i=1}^\alpha \log(U_i)\)

Special cases. \(\chi^2_k = \Gam(k/2, 1/2)\), and this algorithm applies when \(k\) is even; \(\Exp(\lambda) = \Gam(1, \lambda)\).

Gamma from a Sum of Exponentials

Beta from Gamma Random Variables

Theorem (Gamma Ratio \(\to\) Beta). If \(X \sim \Gam(a, \lambda)\) and \(Y \sim \Gam(b, \lambda)\) are independent, then

\[\frac{X}{X + Y} \sim \Betad(a, b) .\]

Essentially, this says that the share of the total contributed by one of two independent gamma variables, when they share the same rate, follows a beta distribution.

Proof sketch. Use the transformation \((U,V) = (X/(X+Y), X+Y)\) with Jacobian \(J = V\).

Algorithm for \(\Betad(a, b)\), \(a, b \in \mathbb{N}\).

  1. Generate \(U_1, \ldots, U_a \iid \Uniform(0,1)\)
  2. Generate \(V_1, \ldots, V_b \iid \Uniform(0,1)\)
  3. \(X = -\sum_{i=1}^a \log(U_i)\), \(Y = -\sum_{i=1}^b \log(V_i)\)
  4. Return \(Z = X/(X+Y)\)

Box–Muller Transformation for Normal

Theorem (Box–Muller)

If \(U_1, U_2 \iid \Uniform(0,1)\), then

\[ \begin{aligned} X_1 &= \sqrt{-2\log U_1}\, \cos(2\pi U_2) \\ X_2 &= \sqrt{-2\log U_1}\, \sin(2\pi U_2) \end{aligned} \]

are independent \(\Normal(0,1)\) random variables.

We can think of this as saying that two uniform draws, passed through the right polar-coordinate recipe, produce two independent standard normals.

Box–Muller: Geometric Intuition

  • \((X_1, X_2)\) has spherical symmetry in \(\mathbb{R}^2\)
  • \(R^2 = X_1^2 + X_2^2 \sim \chi^2_2 = \Exp(1/2)\)
  • \(\Theta = \operatorname{atan2}(X_2, X_1) \sim \Uniform(0, 2\pi)\)
  • Transform from Cartesian to polar coordinates

Recall the \(\Exp(\lambda)\) inverse transform from the first act of this lecture? At \(\lambda = 1/2\), it gives \(R^2 = -2\log U_1 \sim \chi^2_2\), so the squared radius is an inverse-transform exponential draw.

Box–Muller: The Change of Variables

Define the transformation:

\[ \begin{aligned} R &= \sqrt{X_1^2 + X_2^2} = \sqrt{-2\log U_1}\\ \Theta &= \operatorname{atan2}(X_2, X_1) = 2\pi U_2 \end{aligned} \]

Inverting: \(u_1 = e^{-(x_1^2 + x_2^2)/2}\) and \(u_2 = \frac{1}{2\pi}\operatorname{atan2}(x_2, x_1)\).

Compute the Jacobian:

\[J = \det \begin{pmatrix} \frac{\partial u_1}{\partial x_1} & \frac{\partial u_1}{\partial x_2}\\[2pt] \frac{\partial u_2}{\partial x_1} & \frac{\partial u_2}{\partial x_2} \end{pmatrix} = -\frac{1}{2\pi} e^{-(x_1^2 + x_2^2)/2}\]

Box–Muller: The Joint Density Factors

Joint density of \((X_1, X_2)\):

\[f_{X_1,X_2}(x_1,x_2) = |J| \cdot 1 = \frac{1}{2\pi} e^{-(x_1^2 + x_2^2)/2}\]

This factors as two independent \(\Normal(0,1)\) densities.

From Uniform Pairs to Normal Pairs

Marsaglia’s Polar Method

Motivation. Avoid expensive trigonometric functions in Box–Muller.

Algorithm.

  1. Repeat:

    • Generate \(V_1, V_2 \iid \Uniform(-1,1)\)
    • Set \(S = V_1^2 + V_2^2\)

    until \(0 < S \leq 1\).

  2. Set

    \[ \begin{aligned} X_1 &= V_1\sqrt{\frac{-2\log S}{S}}\\ X_2 &= V_2\sqrt{\frac{-2\log S}{S}} \end{aligned} \]

Marsaglia’s Polar Method: Efficiency

  • Expected attempts: \(4/\pi \approx 1.27\)
  • Avoids \(\sin\), \(\cos\) computation
  • Often faster despite rejection step

The number comes from area: the disc occupies \(\pi/4 \approx 0.785\) of the square, so that is the acceptance probability.

Accepted and Rejected Draws in the Unit Square

What You Actually Do in Practice

To generate from a nonstandard target distribution \(F\):

  1. Write down \(F\). If \(F^{-}\) is closed-form, return \(F^{-}(U)\). The PIT licenses this step.
  2. For a discrete target, compute cumulative probabilities \(F_i = \sum_{j \leq i} p_j\), search for the interval \(F_{i-1} < U \leq F_i\), and return \(x_i\). Same theorem.
  3. Otherwise, look for a relationship to a distribution you can already draw from. The General Transformation theorem does the density bookkeeping.

R’s built-in generators already cover the standard families, each with a tuned algorithm of its own: rnorm() inverts \(\Phi\) numerically, and rgamma() uses accept-reject, next week’s method. The recipe matters when the target is nonstandard: a truncated or mixture distribution, a posterior you derived yourself, or any density without a built-in generator.

The Multivariate Normal

Multivariate Normal via a Linear Transformation

Theorem. If \(\mathbf{z} = (z_1, \ldots, z_p)' \sim \Normal_p(\mathbf{0}, \mathbf{I}_p)\) and \(\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}'\) is positive definite (Cholesky, \(\mathbf{L}\) lower triangular), then

\[\mathbf{x} = \boldsymbol{\mu} + \mathbf{L}\mathbf{z} \sim \Normal_p(\boldsymbol{\mu}, \boldsymbol{\Sigma}) .\]

Essentially, this says that multiplying independent standard normals by the Cholesky factor of \(\boldsymbol{\Sigma}\) and shifting by \(\boldsymbol{\mu}\) produces a draw with the desired mean and covariance structure.

Mean and Covariance of the Transform

A linear transformation of a normal vector is normal, so only the mean and covariance need checking.

\[\E[\mathbf{x}] = \E[\boldsymbol{\mu} + \mathbf{L}\mathbf{z}] = \boldsymbol{\mu} + \mathbf{L}\,\E[\mathbf{z}] = \boldsymbol{\mu}\]

\[\Cov(\mathbf{x}) = \Cov(\mathbf{L}\mathbf{z}) = \mathbf{L}\,\Cov(\mathbf{z})\,\mathbf{L}' = \mathbf{L}\mathbf{I}_p\mathbf{L}' = \boldsymbol{\Sigma}\]

Generating a Multivariate Normal Draw

Algorithm.

  1. Compute the Cholesky factor \(\mathbf{L}\) with \(\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}'\)
  2. Generate \(z_i \iid \Normal(0,1)\) for \(i = 1, \ldots, p\)
  3. Return \(\mathbf{x} = \boldsymbol{\mu} + \mathbf{L}\mathbf{z}\)

In R

chol() returns the upper triangular factor \(\mathbf{R}\) with \(\boldsymbol{\Sigma} = \mathbf{R}'\mathbf{R}\), so take L <- t(chol(Sigma)).

Cholesky Reshapes a Standard Normal Cloud

Three Generation Strategies, Three Roles

Result Statement Role
PIT (part 1) \(F^{-}(U) \sim F\) Uniform to any distribution with a tractable \(F^{-}\)
General Transformation \(f_Y(y) = f_X(g^{-1}(y))\,\lvert(g^{-1})'(y)\rvert\) A known draw to a related distribution
MVN Theorem \(\boldsymbol{\mu} + \mathbf{L}\mathbf{z} \sim \Normal_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) Independent normals to a correlated normal vector

Where This Leaves Us

  • \(F^{-}\) carries a uniform draw to any distribution whose CDF you can invert.
  • A change of variables carries a pair of uniforms to a pair of normals.
  • A Cholesky factor carries independent normals to a correlated normal vector.

Every method this week traces back to a uniform draw, and all three need structure you can exploit: a tractable inverse, or a relationship to a distribution you can already sample.

Next: when neither \(F^{-1}\) nor a transformation is available, accept-reject.