Non-Uniform Random Variate Generation: Inversion & Transformations

The foundational theorem underpinning non-uniform variate generation is the Probability Integral Transform.

1. Continuous Probability Integral Transform Theorem

Theorem (Inverse CDF Transform): Let $X$ be a continuous random variable with cumulative distribution function (CDF) $F(x) = P(X \le x)$ that is strictly increasing and continuous on its support. If $U \sim \mathcal{U}(0, 1)$, then the transformed random variable:

has cumulative distribution function exactly equal to $F(x)$.

Mathematical Proof: Since $F$ is strictly increasing, the inverse function $F^{-1}: (0, 1) \to \mathbb{R}$ is uniquely defined and strictly monotonic. For any real $x$:

Applying the strictly increasing function $F$ to both sides of the inequality inside the probability:

Because $U \sim \mathcal{U}(0, 1)$, its cumulative distribution function is $P(U \le u) = u$ for all $u \in [0, 1]$. Therefore:

This establishes that $P(X \le x) \equiv F(x)$, proving the theorem.

Conversely, if $X$ is any continuous random variable with CDF $F(x)$, the transformed variable $Y = F(X)$ follows a standard uniform distribution: $Y \sim \mathcal{U}(0, 1)$.

2. Generalized Inverse for Discrete and Non-Monotonic CDFs

When the CDF $F(x)$ has flat regions (zero probability density) or discrete step jumps (discrete random variables), the standard inverse does not exist. Actuaries and mathematical statisticians define the Generalized Inverse (Quantile Function):

Theorem (Generalized Inverse Transform): For any cumulative distribution function $F(x)$ (continuous, discrete, or mixed), if $U \sim \mathcal{U}(0, 1)$, then $X = F^{-}(U)$ has distribution function $F(x)$.

Proof: By definition of the infimum, the event $F^{-}(u) \le x$ occurs if and only if $u \le F(x)$. Hence:

This universal property ensures that any arbitrary one-dimensional random variable can be sampled via uniform quantile inversion.

Whenever $F(x)$ can be inverted in closed algebraic form, the inverse transform method is optimal in terms of exactness and monotonic synchronization:

1. The Exponential Distribution ($\text{Exp}(\lambda)$)

Let $X \sim \text{Exp}(\lambda)$ with rate parameter $\lambda > 0$. The CDF is:

Setting $U = F(X)$:

Because $U \sim \mathcal{U}(0, 1)$, the complementary variable $1 - U$ also follows $\mathcal{U}(0, 1)$ identically. Thus, in computational implementations, the inversion formula is universally simplified to:

saving one floating-point subtraction per generated variate.

2. The Weibull Distribution ($\text{Weibull}(\alpha, \beta)$)

The Weibull CDF with shape $\alpha > 0$ and scale $\beta > 0$ is:

Inverting $U = F(X)$:

3. The Pareto Distribution ($\text{Pareto}(\alpha, x_m)$)

Widely used in extreme value theory and insurance catastrophe modeling:

Inverting $U = F(X)$:

4. The Cauchy Distribution ($\text{Cauchy}(x_0, \gamma)$)

The standard Cauchy density $f(x) = \frac{1}{\pi (1 + x^2)}$ integrates to:

Inverting yields:

For a discrete random variable taking values $\{x_1, x_2, \dots, x_k\}$ with probability mass function $P(X = x_j) = p_j$ (where $\sum_{j=1}^k p_j = 1$), the generalized inverse corresponds to sequential search through cumulative probability intervals:

To sample $X$:

1. Generate $U \sim \mathcal{U}(0, 1)$.
2. Return $x_j$ where $j$ is the unique index such that:

While binary search reduces lookup time to $\mathcal{O}(\log k)$, large-scale discrete sampling (such as word distributions in natural language processing with $k = 100,000$ categories) demands $\mathcal{O}(1)$ constant-time sampling. This is achieved by Walker's Alias Algorithm (Alastair Walker, 1977).

The Alias Method Mechanics

The alias method represents any discrete probability distribution over $k$ outcomes as an equiprobable mixture of $k$ two-point distributions. Each bucket $j \in \{1, 2, \dots, k\}$ contains:

• A primary outcome $j$ with probability weight $q_j \in [0, 1]$.
• An alias outcome $A_j \in \{1, \dots, k\}$ with complementary probability weight $1 - q_j$.

Sampling Algorithm ($\mathcal{O}(1)$ time):

1. Generate $U \sim \mathcal{U}(0, 1)$.
2. Compute index $I = \lfloor k \cdot U \rfloor + 1 \in \{1, 2, \dots, k\}$.
3. Compute fractional remainder $V = k \cdot U - (I - 1) \in [0, 1)$.
4. If $V < q_I$, return outcome $I$; otherwise, return alias outcome $A_I$.

Walker's alias algorithm requires exactly two floating-point operations and one conditional comparison, completely independent of the alphabet size $k$.

The standard normal distribution $\mathcal{N}(0, 1)$ with density $\phi(z) = \frac{1}{\sqrt{2\pi}} e^{-z^2/2}$ has a cumulative distribution function $\Phi(z) = \int_{-\infty}^z \phi(t)\,dt$ that cannot be expressed in terms of elementary functions, rendering direct algebraic inversion impossible without rational Chebyshev approximations.

In 1958, George Box and Mervin Muller devised an exact geometric transformation that converts a pair of independent standard uniform variates into a pair of independent standard Gaussian variates.

Derivation via Bivariate Polar Coordinates

Consider two independent standard normal random variables $Z_1, Z_2 \sim \mathcal{N}(0, 1)$. Their joint probability density is:

Transform to polar coordinates $(R, \Theta)$ in $\mathbb{R}^2$:

where $R = \sqrt{z_1^2 + z_2^2} \ge 0$ and $\Theta = \arctan(z_2 / z_1) \in [0, 2\pi)$.

The Jacobian determinant of the transformation is:

The joint density in polar coordinates is:

Because the joint density factors into a function of $\theta$ alone and a function of $r$ alone, $R$ and $\Theta$ are mutually independent random variables:

1. $\Theta$ is uniformly distributed on $[0, 2\pi)$:

1. The radius $R$ has cumulative distribution function:

Let $U_2 = 1 - e^{-R^2/2} \sim \mathcal{U}(0, 1)$. Inverting:

Substituting back into Cartesian coordinates produces the Box-Muller Equations:

Both $Z_1$ and $Z_2$ are rigorously independent, standard normally distributed variates.

The Marsaglia-Bray Polar Rejection Method

While mathematically elegant, the standard Box-Muller transform requires evaluating trigonometric functions ($\cos$ and $\sin$), which are computationally expensive on CPU hardware. George Marsaglia and Thomas Bray (1964) eliminated the trigonometric evaluations via acceptance-rejection in the unit disk:

1. Generate $V_1, V_2 \sim \mathcal{U}(-1, 1)$ by setting $V_i = 2U_i - 1$.
2. Compute squared radial distance $W = V_1^2 + V_2^2$.
3. If $W \ge 1$ or $W = 0$, reject and return to Step 1.
4. If $W < 1$ (the point lies inside the unit disk):

Substituting into the Box-Muller formula yields:

The acceptance probability is the area of the unit disk divided by the bounding square: $P(\text{accept}) = \frac{\pi}{4} \approx 0.7854$. Because no trigonometric calls are executed, the polar method is significantly faster.

In multidimensional statistical modeling, financial econometrics, and spatial statistics, actuaries and researchers must generate correlated random vectors:

where $\mathbf{\mu} \in \mathbb{R}^p$ is the mean vector and $\mathbf{\Sigma} \in \mathbb{R}^{p \times p}$ is a symmetric, positive semi-definite covariance matrix.

The Cholesky Factorization Algorithm

By the spectral theorem of linear algebra, any symmetric positive-definite matrix $\mathbf{\Sigma}$ can be uniquely decomposed into the product of a lower triangular matrix $\mathbf{L}$ with strictly positive diagonal entries and its transpose:

Let $\mathbf{Z} = (Z_1, \dots, Z_p)^T$ be a vector of $p$ independent standard normal variates ($Z_j \sim \mathcal{N}(0, 1)$), so that:

Define the affine linear transformation:

We verify the mean and covariance of the transformed vector:

1. Expectation:

2. Covariance:

Because any linear combination of Gaussian random variables is itself Gaussian, $\mathbf{X}$ is guaranteed to follow the exact multivariate normal distribution $\mathcal{N}_p(\mathbf{\mu}, \mathbf{\Sigma})$.

When conducting agent-based simulations, Markov chain simulations with large discrete state spaces, or categorical sampling in natural language models, one must repeatedly draw variates from a discrete probability mass function defined on $K$ categories:

A standard cumulative distribution function (CDF) inversion algorithm constructs the running prefix sums $S_j = \sum_{k=1}^j p_k$ and performs a binary search for a uniform draw $U \sim \mathcal{U}(0, 1)$. This requires $O(\log_2 K)$ search operations per generated variate. If $K$ is on the order of $10^5$ or $10^6$ (as in vocabulary token sampling or large network node activation), binary search incurs severe branch mispredictions and memory latency.

In 1974, Alastair J. Walker introduced the Alias Method, which was subsequently refined into an optimal $O(K)$ setup and $O(1)$ sampling algorithm by Michael D. Vose in 1991.

1. Theoretical Foundation of the Alias Decomposition

The mathematical insight of the Alias Method is that any discrete probability distribution with $K$ outcomes can be expressed as an equiprobable mixture of $K$ two-point distributions.

Imagine $K$ rectangular urns (or bins) of equal width $1/K$ and total height $1.0$. The total area across all $K$ bins is $K \times (1/K) = 1.0$. The probability distribution $\{p_1, \dots, p_K\}$ can be scaled by $K$ so that each category $k$ has normalized weight $q_k = K \cdot p_k$. The average normalized weight is:

Some categories have $q_k < 1.0$ (they are "underfull"), while others have $q_k \ge 1.0$ (they are "overfull"). The Alias Theorem guarantees that one can redistribute the excess mass from overfull categories into underfull categories such that every bin $i \in \{1, \dots, K\}$ contains mass from at most two distinct categories:

1. The primary outcome $i$, occupying probability weight $F_i \in (0, 1]$.
2. An "alias" outcome $A_i \in \{1, \dots, K\}$, occupying the remaining weight $1 - F_i$.

2. The Vose Construction Algorithm

The preprocessing phase partitions categories into two worklists:

• Small: indices with $q_k < 1.0$.
• Large: indices with $q_k \ge 1.0$.

While both lists are non-empty:

1. Remove an index $s$ from Small and an index $l$ from Large.
2. Set the threshold probability table entry $F_s = q_s$.
3. Set the alias table entry $A_s = l$.
4. Reduce the remaining weight of category $l$ by the shortfall of $s$:

1. If the updated $q_l < 1.0$, insert $l$ into Small; otherwise, return $l$ to Large.

Because exactly one element is permanently removed from the active pools at each iteration, the construction terminates in at most $K$ steps, yielding exact $O(K)$ computational and memory complexity.

3. Constant Time Sampling Step

To generate a single realization $X \sim \{p_1, \dots, p_K\}$ in real-time simulation:

1. Generate a single uniform variate $U \sim \mathcal{U}(0, 1)$.
2. Compute the continuous coordinate $V = K \cdot U$.
3. Let $i = \lfloor V \rfloor + 1 \in \{1, \dots, K\}$ be the randomly selected bin index.
4. Let $W = V - \lfloor V \rfloor \in [0, 1)$ be the fractional position within that bin.
5. If $W < F_i$, return outcome $i$; otherwise, return outcome $A_i$.

This execution path consists of a single floating-point multiplication, a floor conversion, one array lookup, and a single branch, achieving strict $O(1)$ time complexity regardless of whether $K = 10$ or $K = 10,000,000$.

The Ratio-of-Uniforms (RoU) method, developed by Albert J. Kinderman and John F. Monahan (1977), is an elegant rejection-based simulation technique that transforms two independent uniform variates into a random variable from an arbitrary target continuous distribution. Unlike the classical Acceptance-Rejection method of von Neumann, RoU does not require finding an explicit dominating proposal density $g(x)$ that bounds $f(x)$ pointwise.

1. Mathematical Formulation

Let $X$ be a continuous random variable with probability density function:

where $h(x) \ge 0$ is a computable unnormalized density function and $Z = \int_{\mathbb{R}} h(x)\,dx$ is an intractable normalizing constant.

Consider the planar region $C_h \subset \mathbb{R}^2$ defined by:

Theorem (Kinderman-Monahan): If a pair $(U, V)$ is distributed uniformly over the region $C_h$, then the ratio:

has exact probability density function $f(x) = h(x) / Z$.

Proof: The area of the region $C_h$ is given by the double integral:

Perform the bivariate change of variables: $x = v/u$ and $u = u$. The Jacobian matrix of the transformation $(u, v) \mapsto (u, x)$ where $v = ux$ is:

Substituting this into the area integral:

Because $(U, V)$ is uniform on $C_h$, its joint density is constant:

Applying the transformation $(U, V) \to (U, X)$ with $V = UX$:

Integrating out the nuisance coordinate $u$ to find the marginal density of $X$:

which is identically the target probability density function.

2. Computational Bounding Box Construction

To sample uniformly from $C_h$, one constructs the minimal enclosing bounding box $[0, u_{\max}] \times [v_{\min}, v_{\max}]$, where:

Provided that $h(x)$ and $x^2 h(x)$ are bounded on $\mathbb{R}$, this bounding box is finite. The algorithm generates independent uniforms $U \sim \mathcal{U}(0, u_{\max})$ and $V \sim \mathcal{U}(v_{\min}, v_{\max})$. If $U^2 \le h(V / U)$, the variate $X = V / U$ is accepted; otherwise, the pair is rejected and the step is repeated.

For standard normal variates $h(x) = e^{-x^2/2}$, we obtain $u_{\max} = 1$, $v_{\min} = -\sqrt{2/e} \approx -0.8578$, and $v_{\max} = \sqrt{2/e} \approx 0.8578$. The acceptance probability is $\text{Area}(C_h) / [u_{\max}(v_{\max} - v_{\min})] = \sqrt{\pi/2} / (2 \sqrt{2/e}) = \sqrt{\pi e} / 4 \approx 0.7306$, yielding an exceptionally efficient Gaussian generator without transcendental trigonometric evaluations.

In spatial statistics, geostatistics, and functional data analysis, simulations require generating continuous sample paths from a Gaussian Process (GP):

where $m(t) = \mathbb{E}[X(t)]$ is the mean function and $K(t, s) = \mathrm{Cov}(X(t), X(s))$ is a symmetric positive semi-definite covariance kernel (such as the Matérn, Squared Exponential, or Ornstein-Uhlenbeck kernel).

1. The Standard Cholesky Factorization Approach

When evaluating the Gaussian process at $N$ discrete locations $\mathbf{t} = (t_1, \dots, t_N)^T$, the joint distribution is multivariate normal:

The classical method computes the Cholesky decomposition of the covariance matrix:

where $\mathbf{L}$ is a lower-triangular matrix. Given an independent standard normal vector $\mathbf{Z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_N)$, the simulated Gaussian process realization is:

While mathematically exact and applicable to irregularly spaced coordinates, Cholesky factorization requires $O(N^3)$ arithmetic operations to compute the decomposition and $O(N^2)$ memory to store the matrix. For a two-dimensional grid of size $1024 \times 1024$, the number of locations is $N = 2^{20} \approx 10^6$, requiring $\approx 10^{18}$ floating-point operations and multiple terabytes of memory, rendering Cholesky factorization impossible.

2. Fast Circulant Embedding and FFT Simulation

When sample locations form a regular equispaced lattice on a stationary domain (where $K(t, s) = \kappa(t - s)$ depends only on the spatial displacement vector), the covariance matrix $\mathbf{\Sigma}$ exhibits a symmetric Toeplitz structure (in 1D) or Block Toeplitz with Toeplitz Blocks (BTTB) structure (in higher dimensions).

Dietrich and Newsam (1993) and Wood and Chan (1994) introduced the Circulant Embedding technique to reduce the computational complexity from $O(N^3)$ to $O(N \log N)$ using the Fast Fourier Transform (FFT).

1. Embedding into a Circulant Matrix: The $N \times N$ Toeplitz covariance matrix $\mathbf{\Sigma}$ is embedded into a larger $M \times M$ circulant matrix $\mathbf{C}$ (typically with $M = 2N$):

where $\mathbf{B}$ is chosen such that the rows of $\mathbf{C}$ are periodic circular shifts of each other.

2. Spectral Decomposition via Discrete Fourier Transform: Every circulant matrix is diagonalized by the unitary Discrete Fourier Transform matrix $\mathbf{F}$:

where $\mathbf{\Lambda} = \mathrm{diag}(\lambda_0, \dots, \lambda_{M-1})$ are the eigenvalues of $\mathbf{C}$, which are computed in $O(M \log M)$ operations by taking the FFT of the first column of $\mathbf{C}$:

3. Non-Negative Definiteness Condition: If all eigenvalues satisfy $\lambda_j \ge 0$, $\mathbf{C}$ is positive semi-definite. A complex normal vector $\mathbf{W} \sim \mathcal{CN}(\mathbf{0}, \mathbf{I}_M)$ is generated with independent real and imaginary standard Gaussian components:

Extracting the first $N$ elements of the real part of $\mathbf{Y}$ yields an exact realization of the Gaussian process $\mathbf{X}$ with covariance $\mathbf{\Sigma}$. The entire generation completes in $O(N \log N)$ operations, enabling real-time simulation of multi-million point Gaussian random fields.

In survival analysis, econometric modeling of censored earnings, and truncated regression (Tobit models), statisticians must draw variates from standard distributions conditioned on belonging to a restricted interval:

When the truncation interval $[a, b]$ falls in the central body of the distribution, naive acceptance-rejection—drawing $X \sim f$ and discarding draws outside $[a, b]$—works well. However, when simulating rare tail events where $a \gg 0$ (e.g., simulating standard normal variates conditioned on $Z \ge 6$), the probability of landing in the interval is $P(Z \ge 6) \approx 9.86 \times 10^{-10}$. Naive rejection sampling requires an average of over one billion simulation draws per accepted variate, resulting in complete computational stalling.

1. Inversion of Truncated Normal CDFs

Using the probability integral transform on the truncated CDF:

Solving for $x$ yields:

While mathematically exact, this formula suffers from catastrophic numerical floating-point underflow when $a$ is large. In standard 64-bit IEEE double precision, $\Phi(6) = 1.0 - \epsilon$ evaluates to $1.0$, causing $\Phi(b) - \Phi(a)$ to round to zero and producing NaN or infinite values.

2. Christian Robert's Optimal Translated Exponential Proposal

Christian P. Robert (1995) designed an optimal rejection algorithm for simulating truncated normal variates $Z \sim \mathcal{N}(0, 1) \mathbb{I}(Z \ge a)$ that achieves bounded, near-unity acceptance rates for arbitrarily large truncation thresholds $a$.

Robert proposed utilizing a translated exponential distribution:

The acceptance ratio is:

To maximize the acceptance probability, we minimize the supremum of the ratio over $z \ge a$. The quadratic exponent achieves its maximum at $z = \lambda$ if $\lambda \ge a$, or at $z = a$ if $\lambda < a$. Robert proved that the optimal scale parameter $\lambda^*$ is:

Under this optimal parameter, the acceptance ratio simplifies remarkably:

The simulation algorithm proceeds as follows:

1. Generate $E_1, E_2 \sim \text{Exp}(1)$ independently.
2. Set proposed value $Z = a + E_1 / \lambda^*$.
3. If $E_2 \ge \frac{(Z - \lambda^*)^2}{2}$, accept $Z$; otherwise reject and repeat.

Because the acceptance probability satisfies $P(\text{accept}) \ge 0.76$ for all $a \ge 0$, and approaches $1.0$ as $a \to \infty$, this algorithm produces exact extreme tail draws in milliseconds with zero numerical underflow.

Point processes model the stochastic arrival of discrete events in continuous time $t \in [0, T]$. In seismology (aftershocks), high-frequency financial microstructure (order book transactions), and epidemiology (contagion propagation), event arrivals exhibit temporal clustering: the occurrence of an event increases the instantaneous probability of subsequent events. Alan G. Hawkes (1971) formalized this behavior in the Hawkes Self-Exciting Process.

The conditional intensity function $\lambda(t \mid \mathcal{H}_t)$ given event history $\mathcal{H}_t = \{t_i : t_i < t\}$ is:

where $\mu > 0$ is the baseline background rate, $\alpha > 0$ is the excitation magnitude, and $\beta > \alpha$ is the exponential decay rate. The ratio $n = \alpha / \beta < 1$ represents the branching ratio, ensuring stability.

The Lewis-Shedler Thinning Algorithm for Hawkes Processes

Because $\lambda(t)$ jumps discontinuously at each arrival time and decays continuously between arrivals, direct inversion of the cumulative intensity is intractable. P. A. W. Lewis and G. S. Shedler (1979) introduced the Thinning Algorithm (a continuous-time Poisson analog of rejection sampling):

1. Set current time $t = 0$ and event counter $k = 0$.
2. Compute the maximum instantaneous intensity immediately after the most recent event:

1. Generate a candidate inter-arrival time from a homogeneous Poisson process with rate $\lambda_{\max}$:

If $t_{\text{cand}} > T$, terminate simulation.

1. Evaluate the exact decaying intensity at the candidate time:

1. Draw $U \sim \mathcal{U}(0, 1)$. If $U \le \frac{\lambda(t_{\text{cand}})}{\lambda_{\max}}$, the candidate arrival is accepted:

Otherwise, the candidate arrival is "thinned" (rejected), and we set $t \leftarrow t_{\text{cand}}$ without recording an arrival.

1. Repeat from Step 2 until the time horizon $T$ is reached.

While the Box-Muller transform and the polar method generate standard Gaussian variates effectively, both algorithms require evaluating expensive transcendental functions: Box-Muller computes square roots, natural logarithms, and trigonometric sines and cosines, consuming dozens of clock cycles per generated pair.

George Marsaglia and Wai Wan Tsang (2000) introduced the Ziggurat Algorithm, which generates continuous random variates (particularly standard Normal and Exponential variables) using almost exclusively table lookups and integer comparisons, executing in fewer than 10 CPU clock cycles on modern hardware.

1. Geometric Construction of the Ziggurat

Let $f(x)$ be a monotonically decreasing probability density function on $[0, \infty)$ (such as the positive half of the standard normal density $f(x) = \sqrt{2/\pi} e^{-x^2/2}$).

The area under $f(x)$ is partitioned into $M$ horizontal rectangular slices (typically $M = 256$, fitting into an 8-bit index). Slices are indexed $i = 0, 1, \dots, M-1$ from bottom to top:

• Each rectangular layer $i$ extends horizontally from $0$ to $x_i$ and vertically from $y_{i-1}$ to $y_i = f(x_i)$.
• Crucially, the endpoints $x_i$ are precomputed such that every layer has identical area $V$:

The bottom slice $i = 0$ is special: it consists of a central rectangle of width $x_0$ and height $y_0 = f(x_0)$, plus an infinite right-hand tail $T = \int_{x_0}^\infty f(x)\,dx$ such that the combined area equals $V$.

Because the layers have equal area $V$, selecting a layer $i$ uniformly at random from $\{0, 1, \dots, M-1\}$ gives an exact equiprobable mixture.

2. The Rapid Sampling Loop

To generate a standard normal variate:

1. Generate an independent 32-bit random integer $J$.
2. Extract the lower 8 bits of $J$ to select layer index $i \in \{0, \dots, 255\}$.
3. Use the remaining 24 bits of $J$ (converted to a uniform fraction $U \in [0, 1)$) to propose a horizontal coordinate:

4. Fast Interior Acceptance Test: If $X < x_{i+1}$, the proposed point lies strictly inside the core interior of the ziggurat layer (underneath the density curve). Accept $X$ immediately!

5. Slow Boundary Test: If $X \ge x_{i+1}$, the proposed point falls in the outer margin near the density boundary. Only then does the algorithm evaluate:

If $Y \le f(X)$, accept $X$; otherwise, reject and repeat.

6. Tail Draw: If layer $i = 0$ was selected and $X \ge x_0$, sample directly from the tail using Robert's transformed exponential tail algorithm.

Because $M = 256$ slices provide an exceptionally tight stepped approximation to the smooth Gaussian curve, the fast interior acceptance test in Step 4 succeeds over $98.5\%$ of the time. The ziggurat algorithm avoids transcendental evaluations on $98.5\%$ of all generated variates, making it the fastest known algorithm for standard continuous distributions.

In spatial navigation, crystallography, structural bioinformatics (protein backbone torsion angles), and paleomagnetism, random variables are not unconstrained Euclidean vectors, but unit directions residing on non-Euclidean manifolds such as the unit sphere $\mathbb{S}^{d-1} = \{ \mathbf{x} \in \mathbb{R}^d : \|\mathbf{x}\|_2 = 1 \}$ or the rotation group $\mathrm{SO}(3)$.

1. The Von Mises-Fisher Distribution on the Sphere

The Von Mises-Fisher (vMF) distribution models directional clustering on $\mathbb{S}^{d-1}$ around a mean direction $\boldsymbol{\mu} \in \mathbb{S}^{d-1}$ with concentration parameter $\kappa \ge 0$:

where $C_d(\kappa) = \frac{\kappa^{d/2 - 1}}{(2\pi)^{d/2} I_{d/2 - 1}(\kappa)}$ is the normalizing constant involving the modified Bessel function of the first kind.

To simulate a vMF variate $\mathbf{X}$ efficiently without multidimensional rejection (Wood 1994):

1. Decompose the target vector into the parallel component $W = \boldsymbol{\mu}^T \mathbf{X} \in [-1, 1]$ and an orthogonal tangent vector $\mathbf{V} \in \mathbb{S}^{d-2}$:

where $\mathbf{Q}$ is an orthonormal basis spanning the orthogonal complement of $\boldsymbol{\mu}$.

1. The scalar coordinate $W$ has density $f_W(w) \propto (1 - w^2)^{(d-3)/2} e^{\kappa w}$, which is sampled using an acceptance-rejection scheme with an envelope based on the Beta distribution.
2. The orthogonal unit vector $\mathbf{V}$ is sampled uniformly over the equatorial sphere $\mathbb{S}^{d-2}$ by drawing standard Gaussian vectors and normalizing to unit length: $\mathbf{V} = \mathbf{Z} / \|\mathbf{Z}\|_2$ with $\mathbf{Z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_{d-1})$.

In computational linguistics and generative modeling, sampling next-token predictions involves drawing from categorical distributions with vocabulary sizes $K \in [32,000, 256,000]$.

Comparing the execution profile across $10^7$ independent categorical draws:

• Cumulative Binary Search: Requires $O(\log_2 K) \approx 15$ to 18 memory lookups and branch evaluations per draw, triggering extensive CPU branch mispredictions and taking 14.2 seconds.
• Vose's Alias Method: Incurs an initial one-time table setup cost of $0.08$ seconds ($O(K)$), but executes each subsequent sample in strict $O(1)$ time with a single multiplication and array lookup, completing all $10^7$ draws in 0.38 seconds—a 37-fold speedup.

For static or infrequently updated probability vectors, the Alias Method represents the optimal industry standard.

Unit 2 examined the transformation of uniform random variables into arbitrary non-uniform continuous and discrete distributions. We derived the Probability Integral Transform and established its exact inversion mechanics for closed-form cumulative distribution functions. We formulated the Acceptance-Rejection method of John von Neumann, proving that the expected number of proposal evaluations is exactly the bounding constant $c = \sup [f(x)/g(x)]$. We detailed Box-Muller and Marsaglia-Bray polar transformations for Gaussian generation, alongside Marsaglia and Tsang's Ziggurat algorithm, which achieves over $98.5\%$ direct interior acceptance through precomputed horizontal planar layers. For large discrete categorical spaces, we analyzed Vose's formulation of the Walker Alias Method, showing that pre-allocating an alias table and probability threshold table in $O(K)$ time permits subsequent draws in strict $O(1)$ time. Finally, we explored the Kinderman-Monahan Ratio-of-Uniforms method and fast circulant embedding for generating high-dimensional Gaussian process paths via Fast Fourier Transforms.

In constrained econometric models and spatial Gaussian processes with inequality barriers, one must sample from a high-dimensional truncated Gaussian distribution $\mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{\Sigma})$ subject to linear constraints $\mathbf{a} \le \mathbf{B} \mathbf{X} \le \mathbf{b}$. Naive rejection sampling fails completely because the probability of an unconstrained Gaussian draw satisfying dozens of joint linear inequalities is virtually zero ($P \le 10^{-15}$).

Zdravko Botev (2017) developed an exact minimax tilting algorithm that combines exponential twisting with smooth unconstrained optimization to achieve acceptance rates bounded away from zero regardless of dimension $d$. Alternatively, componentwise Gibbs sampling updates each scalar coordinate $X_j$ conditional on all remaining coordinates $\mathbf{X}_{-j}$ from a univariate truncated Gaussian distribution using Robert's optimal translated exponential proposal. This allows simulation of thousands of correlated variables under strict convex polyhedral bounds in milliseconds.