How to Sample a Gaussian Distribution: Methods, Algorithms, and Performance

Scene Description

Gaussian distribution, also known as the normal distribution, is a fundamental probability distribution in mathematics, physics, and engineering. In practice we often need to draw samples from it. Understanding the underlying algorithms deepens knowledge of probability and statistics, and comparing different sampling methods gives a concrete impression of their performance.

Problem Description

If you were asked to implement a Gaussian random number generator, how would you do it?

Answer and Analysis

Assume a random variable z follows the standard normal distribution N (0,1). By setting x = σ·z + μ , the variable x follows a Gaussian distribution N (μ, σ²). Thus any Gaussian can be obtained from a standard normal via scaling and shifting, so we only need to sample the standard normal. All sampling methods start from a uniform random number generator, typically implemented with a linear congruential generator.

Common sampling methods include Inverse Transform, Rejection Sampling, Importance Sampling (and its resampling variant), and Markov Chain Monte Carlo. For the Gaussian case we consider the following approaches:

Inverse Transform Method

The inverse of the error function erf (x) has no closed‑form expression, making direct inverse transform cumbersome.

Box‑Muller Algorithm

To avoid the non‑elementary inverse, Box‑Muller generates two independent standard normal variables x and y from two uniform variables. By converting to polar coordinates (r, θ) and sampling r using the inverse CDF of its distribution and θ uniformly on [0, 2π), we obtain (x, y) that follow the standard normal. The procedure is illustrated below:

Box‑Muller requires trigonometric functions, which are relatively costly.

Marsaglia Polar Method

This method avoids trigonometric calls, making it faster. The sampling steps are shown:

Rejection Sampling

We can choose a reference distribution with an easy inverse CDF, such as the exponential distribution (λ=1), to cover the target normal distribution. Because the exponential is defined for x ≥ 0 while the normal spans (‑∞, +∞), we exploit symmetry to map samples between half‑ and full‑axes. The acceptance probability should be maximized by keeping the scaling factor M small.

Ziggurat Algorithm

This is an efficient rejection‑sampling technique that approximates the target distribution with multiple stacked rectangles. Although the implementation looks more involved, it is highly efficient. See the diagram for the Ziggurat layout:

Summary and Extension

Many Gaussian sampling methods exist; we listed the most common ones. In interviews, candidates need only know one or two methods, allowing interviewers to probe deeper into theory, advantages, disadvantages, and performance. Extensions include sampling truncated Gaussian distributions and handling high‑dimensional variables, where rejection sampling becomes inefficient and alternative strategies are required.