Normal distribution

The normal distribution, also called the Gaussian distribution, is a continuous probability distribution defined by two parameters, a mean $\mu$ and a variance $\sigma^2$, whose density $f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$ is the symmetric bell curve. It is the most widely used distribution in statistics and machine learning: it is the limiting distribution of sums of many independent finite-variance random variables (the central limit theorem), the maximum-entropy distribution given a fixed mean and variance ^[4], and the noise, prior, and weight-initialization model assumed across modern deep learning. Its density is the bell curve that almost every introductory statistics course pins to its first chapter, and the one that almost every machine learning paper assumes somewhere in the fine print.

In shorthand, a random variable $X$ is normally distributed when $X \sim N(\mu, \sigma^2)$ and its probability distribution has density:

$$f(x; \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$$

The distribution is symmetric around $\mu$, has thin (sub-exponential) tails, and is closed under linear combinations. That last property is what makes it cooperate so well with everything else in statistics. Sums of independent normals stay normal. Linear transformations of normals stay normal. Marginals and conditionals of multivariate normals stay normal. There is no other distribution that is this well-behaved while still being non-trivial, which is part of why it has become the default model for noise, errors, prior beliefs, latent variables, and weight initializations across modern machine learning.

The other reason for its dominance is the central limit theorem, which says that the sum of many independent finite-variance contributions tends to a normal distribution regardless of the shape of the individual pieces. So even when individual sources of randomness are not normal, their aggregates often are. That is why measurement errors, polling results, residuals from a well-specified linear regression, and many natural phenomena all end up looking Gaussian.

Whether they actually are Gaussian is a separate question, and one this article also tries to take seriously. Plenty of real-world data is not normal, and assuming it is can blow up in your face.

What are the parameters of a normal distribution?

Density and parameters

A random variable $X$ follows a normal distribution with mean $\mu$ and variance $\sigma^2 > 0$ if its probability density function is:

$$f(x; \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)$$

for $x \in \mathbb{R}$. The two parameters have the usual interpretation. The mean $\mu$ is a real number that locates the center of the distribution. The variance $\sigma^2$ is positive and controls the spread; the standard deviation $\sigma = \sqrt{\sigma^2}$ is the more common scale measure because it shares units with $x$.

The density integrates to 1 over the real line, which can be verified using the Gaussian integral $\int \exp(-x^2) \, dx = \sqrt{\pi}$ over $(-\infty, \infty)$. The normalizing constant $\frac{1}{\sigma\sqrt{2\pi}}$ exists precisely to make this integral equal 1.

What is the standard normal distribution?

The standard normal distribution is the special case $\mu = 0$, $\sigma^2 = 1$. It is written $N(0, 1)$ and its density is usually denoted $\varphi(z)$:

$$\varphi(z) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{z^2}{2}\right)$$

Any normal X ~ N(μ, σ²) can be standardized by subtracting the mean and dividing by the standard deviation. The standardized variable $Z = (X - \mu) / \sigma$ follows $N(0, 1)$, and the value $z$ is called the z-score of $x$. Statistical tables, and the entire ecosystem of hypothesis tests built around the normal, lean on this trick because there is only one standard normal distribution to tabulate. The same standardize-to-$N(0, 1)$ operation is what machine learning practitioners call feature standardization or z-scoring, a near-universal preprocessing step before fitting linear models, support vector machines, or neural networks.

Cumulative distribution function

The CDF is:

$$\Phi(z) = \int_{-\infty}^{z} \varphi(t) \, dt$$

Unlike the density, the CDF has no closed-form expression in elementary functions. It is usually written in terms of the error function erf:

$$\Phi(z) = \frac{1}{2}\left(1 + \mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)\right)$$

This is mostly a notation issue. Numerical libraries compute Φ via rational approximations, asymptotic expansions for the tails, or piecewise polynomial fits. NumPy uses scipy.special.ndtr, which calls into a C implementation derived from Cephes; PyTorch's torch.special.ndtr does the analogous thing on the GPU.

Mean, median, mode, variance

For the normal distribution, the mean, median, and mode all coincide at $\mu$. This is one of the few distributions where this is true exactly, and it is part of why "average" is so often treated as a synonym for "typical" in everyday speech. The variance is $\sigma^2$, so the standard deviation is $\sigma$. Skewness is 0 because the density is symmetric, and the kurtosis is 3 (excess kurtosis 0). Kurtosis 3 is the reference value other distributions are compared against, which is why excess kurtosis is reported relative to it.

What is the 68-95-99.7 rule?

A useful rule of thumb for the normal distribution, known as the empirical rule or three-sigma rule, is that about 68% of the mass falls within one standard deviation of the mean, about 95% within two, and about 99.7% within three ^[5]. The full quantile picture is a bit more nuanced, and the exact values are 0.6827, 0.9545, and 0.9973 ^[5].

Scroll sideways for more →

The exact figures to ten digits are 0.6826894921 for one sigma, 0.9544997361 for two, and 0.9973002039 for three ^[5]. The last few rows are why physicists insist on "five sigma" for discoveries (about 1 in 3.5 million under the null) and why finance practitioners get nervous when they see "six sigma" events happening in a single decade. If your data is actually Gaussian, six-sigma days should not happen on a human timescale. If they do, your data is not Gaussian.

The round numbers people quote (68, 95, 99.7) are slightly off from the actual values. The two-sigma probability is closer to 95.45%, not 95% ^[5]. The exact two-sided 95% interval uses $1.95996\sigma$, and the exact 99% interval uses $2.5758\sigma$. Most reference tables stop at four decimal places of z because that is all anyone ever needs in practice.

When was the normal distribution discovered?

De Moivre and the binomial limit

The earliest derivation of the normal density appears in Abraham de Moivre's 1733 paper Approximatio ad Summam Terminorum Binomii (a + b)^n in Seriem expansi, a seven-page Latin pamphlet circulated privately on 12 November 1733 and later folded into the 1738 second edition of his The Doctrine of Chances ^[1][7]. De Moivre was studying the binomial distribution, the probability of getting $k$ heads in $n$ flips of a fair or biased coin, and he wanted to evaluate the central terms when $n$ is large.

For a symmetric binomial with $p = 1/2$, the central term $C(n, n/2) \cdot 2^{-n}$ is hard to compute directly when $n$ is large, so de Moivre used Stirling's approximation (which Stirling and de Moivre had both worked on) and showed that the binomial probability mass near the center decays approximately as $\exp(-2k^2 / n)$ where $k$ measures distance from the mean in standard-deviation units. This is the normal density, although de Moivre did not call it that, did not parameterize it in modern (μ, σ²) notation, and saw it as a tool for binomial calculations rather than a distribution in its own right. His pamphlet nonetheless gave the first derivation of the curve, the first method for computing normal tail probabilities, and the first recognition that the standard deviation is the natural unit for measuring deviations from the mean, more than seventy years before Gauss ^[7].

De Moivre's result is the first instance of what would later be called the central limit theorem for the special case of i.i.d. Bernoulli sums.

Laplace and the central limit theorem

The more general statement came from Pierre-Simon Laplace (Laplace). In his 1774 Mémoire sur la probabilité des causes par les événements, and culminating in the 1810 central limit theorem and the 1812 Théorie analytique des probabilités, Laplace proved that sums of independent random variables, suitably normalized, converge to the normal distribution under fairly broad conditions ^[2][8]. He used this result to argue that astronomical observation errors, which are sums of many small independent perturbations, should be expected to follow what we now call a Gaussian distribution.

Laplace's contribution is often underrated next to Gauss's, but the central limit theorem is the deeper reason the normal distribution shows up everywhere. De Moivre's coin-flipping result is a special case; Laplace's theorem is the principle.

Gauss and astronomy

The distribution gets its most common name from Carl Friedrich Gauss (Gauss). In his 1809 Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium, Gauss developed the method of least squares as a tool for fitting orbits to noisy astronomical observations ^[3]. He showed that if observation errors are normally distributed, the maximum-likelihood estimate of the orbital parameters coincides with the least-squares estimate, and that the resulting estimator is in a precise sense optimal.

Gauss arrived at the normal density essentially by reverse-engineering the assumption: if the arithmetic mean of repeated measurements is the most probable value of the true quantity (a postulate that was already widely used by astronomers), then the error distribution must take the form $\exp(-h^2 \varepsilon^2)$ up to a constant ^[3][8]. The argument is short, surprisingly clever, and not quite a proof of the centrality of the normal distribution, but it gave astronomers a principled reason to use least squares and the bell curve became permanently associated with his name.

The Gauss-Laplace synthesis was already complete by the 1810s. By the time the 19th century was over, the normal distribution had become the default model for measurement error in essentially every quantitative science.

Quetelet and the social sciences

Adolphe Quetelet in the 1830s and 1840s did something that, depending on how charitable you are, was either inspired or deeply confused. He took the normal distribution out of astronomy and applied it to human characteristics: heights of conscripts, chest sizes of Scottish soldiers, crime rates, suicide rates. He observed that many of these distributions looked roughly bell-shaped, and proposed his concept of l'homme moyen, the "average man," as a kind of statistical ideal type whose deviations represented social or biological error ^[4].

This was the moment when the normal distribution became a generic model for empirical phenomena rather than a specifically astronomical tool. Quetelet's interpretation of variation as "error around the average man" has not aged well, but his methodological move (that the normal distribution applies to social data as well as physical measurements) was enormously influential and helped launch the field of mathematical statistics.

Why is it called the "normal" distribution?

The word "normal" in this technical sense was coined relatively late. Karl Pearson is usually credited with popularizing it in 1893, although Charles Sanders Peirce, Francis Galton, and Wilhelm Lexis had all used variants of the term in the 1870s and 1880s ^[9]. Pearson himself later expressed regret about the choice, writing in 1920: "Many years ago I called the Laplace-Gaussian curve the normal curve, which name, while it avoids an international question of priority, has the disadvantage of leading people to believe that all other distributions of frequency are in one sense or another 'abnormal.'" ^[10]

He was right. The connotation has caused real damage in fields that uncritically apply the Gaussian as a default. Calling something "normal" suggests that deviations from it are pathological, which is exactly the wrong intuition for distributions like the Cauchy or the t-distribution where the "abnormality" is simply having a heavier tail.

Properties

Closed under linear combinations

If $X \sim N(\mu_X, \sigma_X^2)$ and $Y \sim N(\mu_Y, \sigma_Y^2)$ are independent, then for any constants $a$ and $b$:

$$aX + bY \sim N(a\mu_X + b\mu_Y, a^2\sigma_X^2 + b^2\sigma_Y^2)$$

This property generalizes to any finite number of independent normals, and it is the workhorse identity behind almost every analytic calculation involving Gaussians. It is also why the sample mean of a normal sample is exactly normal, not merely approximately normal: a sum of $n$ independent $N(\mu, \sigma^2)$ variables divided by $n$ is $N(\mu, \sigma^2/n)$ on the nose.

Maximum entropy

Among all probability distributions on the real line with a given mean and variance, the normal distribution has the maximum differential entropy ^[4]. This is a precise sense in which the Gaussian is the "least informative" or "most uncertain" distribution consistent with first and second moments. From the perspective of information theory, if all you know about a quantity is its mean and variance, the most honest distribution to assume is the normal one.

This is one of the deepest justifications for the Gaussian as a default model. It is not just convenient, it is the unique solution to a maximum-entropy problem under second-moment constraints.

Stability and the central limit theorem

The normal distribution is one of only three stable distributions with finite variance (the others being trivial cases). Stability here means that a sum of i.i.d. normals is itself normal up to a location and scale shift. This is the algebraic version of the central limit theorem: if you keep adding finite-variance contributions, you converge to the unique stable law with finite variance, which is the normal.

Conjugacy in Bayesian inference

The normal distribution is its own conjugate prior under Gaussian likelihood with known variance. If you assume $X \mid \mu \sim N(\mu, \sigma^2)$ with $\sigma^2$ known, and you place a prior $\mu \sim N(\mu_0, \tau_0^2)$, then after observing $n$ samples with sample mean $\bar{x}$ the posterior is also normal:

$$\mu \mid \text{data} \sim N(\mu_n, \tau_n^2)$$

where the posterior mean $\mu_n$ is a precision-weighted average of the prior mean and the data mean, and the posterior precision $1/\tau_n^2$ is the sum of the prior precision $1/\tau_0^2$ and the data precision $n/\sigma^2$ ^[11][12]. This is the cleanest possible update rule in Bayesian inference and is the reason Gaussian assumptions show up so often in Bayesian models.

Memoryless tails

The normal distribution has tails that decay faster than any exponential. For large $z$, $P(Z > z) \approx \varphi(z) / z$, which means the density falls off as $\exp(-z^2/2)$. This is much faster decay than the exponential or Laplace distributions, and it is one of the reasons the Gaussian is sometimes a bad model for real data: extreme events that look impossible under a Gaussian assumption may actually be common in heavier-tailed distributions.

The multivariate normal

The extension of the normal distribution to d dimensions is the multivariate normal. A random vector $X \in \mathbb{R}^d$ is multivariate normal with mean vector $\mu \in \mathbb{R}^d$ and covariance matrix $\Sigma \in \mathbb{R}^{d \times d}$ (symmetric positive definite) if its joint density is:

$$f(x; \mu, \Sigma) = \frac{1}{(2\pi)^{d/2} \, \lvert \Sigma \rvert^{1/2}} \exp\left(-\frac{1}{2}(x - \mu)^\top \Sigma^{-1} (x - \mu)\right)$$

The quantity $(x - \mu)^\top \Sigma^{-1} (x - \mu)$ is the squared Mahalanobis distance from $x$ to $\mu$. Level sets of the density are ellipsoids whose principal axes are the eigenvectors of $\Sigma$ and whose lengths are proportional to $\sqrt{\lambda_i}$, where $\lambda_i$ are the eigenvalues.

Marginal and conditional distributions of a multivariate normal are themselves normal, with means and covariances given by closed-form expressions in terms of the original $\mu$ and $\Sigma$ ^[11]. This is the property that makes Gaussian processes work, since a GP can be defined as a collection of random variables of which any finite subset has a joint multivariate normal distribution, and conditioning on observations is just a linear-algebra operation.

The special case $\Sigma = \sigma^2 I$ gives an isotropic Gaussian, which is what people mean when they say "sample from $N(0, I)$." In high dimensions, isotropic Gaussians have a lot of geometry that is not at all obvious from the 1D picture: most of the mass concentrates in a thin spherical shell at radius $\sqrt{d}$, and any two random samples are nearly orthogonal. This concentration of measure is one of the reasons high-dimensional Gaussians behave so differently from intuition built on 2D bell curves.

What is the central limit theorem?

The central limit theorem is the main reason the Gaussian is everywhere. The classical Lindeberg-Lévy version says: if $X_1, X_2, \ldots, X_n$ are i.i.d. random variables with mean $\mu$ and finite variance $\sigma^2$, then the standardized sample mean

$$\sqrt{n} \, (\bar{X}_n - \mu) / \sigma$$

converges in distribution to $N(0, 1)$ as $n \to \infty$ ^[11][12].

The theorem has dozens of generalizations. Lyapunov and Lindeberg-Feller relax the i.i.d. assumption to allow non-identically-distributed variables under a uniform integrability condition. Martingale CLTs handle dependent sequences. Stable distribution theory tells you what happens when the finite-variance assumption fails: the limit is no longer Gaussian but a stable law with infinite variance.

What this means in practice for machine learning: any quantity that is a sum or average of many small independent contributions (gradients across a mini-batch, neuron activations under random weights, errors in an ensemble of predictors) will tend toward a Gaussian as the number of contributions grows. This is why Gaussian noise models often work even when the underlying data-generating process is not itself Gaussian; you are typically modeling residuals or aggregates rather than raw observations.

The theorem can also mislead. The CLT only kicks in once $n$ is large, and the rate of convergence depends on higher moments of the underlying distribution. For heavy-tailed inputs the convergence is painfully slow, and using a Gaussian approximation for $n = 30$ (the textbook rule of thumb) on a power-law-tailed dataset will give wildly wrong tail probabilities.

How do you sample from a normal distribution?

Most computational work with normals starts from a uniform random number generator and converts uniforms to normals. The classical methods are:

Box-Muller

The Box-Muller transform, published by George E. P. Box and Mervin E. Muller in 1958, takes two independent uniform samples $U_1, U_2 \in (0, 1)$ and produces two independent $N(0, 1)$ samples ^[13]:

$$Z_1 = \sqrt{-2 \ln U_1} \cos(2\pi U_2)$$ $$Z_2 = \sqrt{-2 \ln U_1} \sin(2\pi U_2)$$

The derivation is a polar-coordinate change of variables on the joint density of two independent standard normals. It is exact (modulo floating-point), but it requires calls to log, sqrt, sin, and cos, which used to be expensive on older hardware.

Marsaglia polar method

George Marsaglia's polar method is a rejection-sampling variant of Box-Muller that avoids the trigonometric calls. You sample $U_1, U_2$ uniformly in $(-1, 1)$, reject if $U_1^2 + U_2^2 > 1$, and otherwise compute the Box-Muller transform using $U_1 / \sqrt{s}$ and $U_2 / \sqrt{s}$ as the cosine and sine, where $s = U_1^2 + U_2^2$. It is faster than Box-Muller in practice on most hardware but has a rejection rate of about 21% (the probability that a uniform point in the square falls outside the inscribed unit disk is $1 - \pi/4 \approx 0.2146$).

Ziggurat

The Ziggurat algorithm, published by Marsaglia and Wai Wan Tsang in 2000, is the method most modern numerical libraries use under the hood ^[14]. It approximates the normal density with a stack of horizontal rectangles and a tail region, samples uniformly from the union, and accepts most samples without any expensive function evaluation. NumPy's default_rng().standard_normal() and PyTorch's torch.randn both use Ziggurat-style methods, sometimes with vendor-specific tweaks.

For multivariate normals N(μ, Σ), the standard procedure is to compute the Cholesky factor $\Sigma = L L^\top$, draw $Z \sim N(0, I)$ using one of the methods above, and return $X = \mu + L Z$. For very high-dimensional $\Sigma$ this is the bottleneck, and various structured-covariance tricks (low-rank plus diagonal, Toeplitz, banded) are used in practice to avoid the full $O(d^3)$ Cholesky.

Inference

Maximum likelihood

Given $n$ i.i.d. samples $x_1, \ldots, x_n$ from $N(\mu, \sigma^2)$, the maximum likelihood estimators are:

$$\hat{\mu} = \frac{1}{n} \sum x_i = \bar{x}$$ $$\hat{\sigma}^2 = \frac{1}{n} \sum (x_i - \bar{x})^2$$

The MLE for the variance is biased: $\mathbb{E}[\hat{\sigma}^2] = \frac{n - 1}{n} \sigma^2$. The unbiased estimator divides by $n - 1$ instead of $n$:

$$s^2 = \frac{1}{n - 1} \sum (x_i - \bar{x})^2$$

This is Bessel's correction, named after Friedrich Bessel. The factor of $(n - 1)$ compensates for the fact that $\bar{x}$ is itself estimated from the data, which uses up one degree of freedom. The biased and unbiased estimators differ noticeably for small $n$ and become indistinguishable as $n \to \infty$.

For linear regression with Gaussian errors, the MLE for the regression coefficients coincides with ordinary least squares. This is the connection Gauss exploited in 1809: maximizing a Gaussian likelihood is mathematically the same as minimizing the mean squared error ^[3]. Squared error is not just a convenient loss function, it is the negative log likelihood under a Gaussian noise model.

Z-test and t-test

If you know σ in advance, the z-test uses the fact that $(\bar{X}_n - \mu_0) / (\sigma / \sqrt{n}) \sim N(0, 1)$ under the null hypothesis $\mu = \mu_0$. You compare the observed z to the standard normal CDF and read off a p-value.

In practice σ is almost never known, so you estimate it from the sample using $s$, and the resulting test statistic $(\bar{X}_n - \mu_0) / (s / \sqrt{n})$ follows Student's t-distribution with $n - 1$ degrees of freedom. This is the t-test, derived by William Sealy Gosset in 1908 under the pen name "Student." For $n$ much larger than 30 the t-distribution is essentially indistinguishable from the standard normal, but for small samples the t has heavier tails to account for the additional uncertainty in s.

The F-test for variance ratios and the chi-squared test for variance both also assume Gaussianity. All three of these tests are standard fare in statistical practice and they all rest on the same set of distributional assumptions about the underlying data.

Confidence intervals

A two-sided 95% confidence interval for $\mu$ when $\sigma$ is known is:

$$\bar{x} \pm 1.96 \, \sigma / \sqrt{n}$$

When $\sigma$ is unknown, replace 1.96 with the appropriate t-quantile (for $n = 30$ this is about 2.04; for $n = 10$ it is 2.26). The interval is constructed so that, over many repeated samples, 95% of the constructed intervals will contain the true $\mu$.

How do you test whether data is normally distributed?

A persistent problem with normality assumptions is that they are usually not checked. The following table summarizes the standard tests practitioners use to decide whether a normal model is reasonable.

Scroll sideways for more →

The two most common tests in machine learning practice are Shapiro-Wilk (built into scipy.stats.shapiro) and the visual QQ plot. The QQ plot in particular is often more useful than a hypothesis test, because it tells you not just whether the data deviates from normal but where (in the tails, in the center, only in one direction) and by how much.

There is a subtle problem with normality tests in modern data: with $n = 10000$ or more, virtually any real dataset will fail a Shapiro-Wilk test even if the deviation from normal is small enough to be irrelevant for downstream inference. The hypothesis test is not the right question; the right question is whether the deviation matters for whatever you are doing.

How is the normal distribution used in machine learning?

The Gaussian shows up in machine learning in more places than would fit in any one section. The following list covers the major ones.

Scroll sideways for more →

How does Gaussian weight initialization work?

When you initialize the weights of a deep neural network, you almost always sample them from some kind of Gaussian (or symmetric uniform with matched variance). Xavier Glorot and Yoshua Bengio's 2010 paper Understanding the difficulty of training deep feedforward neural networks showed that for tanh networks, sampling weights with variance $2 / (n_{\text{in}} + n_{\text{out}})$ keeps the variance of activations and gradients roughly constant across layers, which prevents the vanishing or exploding gradient problem ^[15]. Kaiming He et al.'s 2015 follow-up Delving Deep into Rectifiers derived the analogous variance for ReLU networks, $2 / n_{\text{in}}$, accounting for the fact that ReLU zeros out about half its inputs, so the variance has to be doubled to compensate ^[16].

The practical effect is dramatic. Without principled initialization, training networks deeper than about ten layers was essentially impossible in 2010. With Glorot or He initialization, plus residual connections and normalization layers, modern networks can be hundreds of layers deep and still trainable. The mathematical content of the initialization is just a Gaussian with carefully chosen variance, but the engineering consequences were enormous. See weight initialization for the full treatment.

Batch normalization

Batch normalization, introduced by Sergey Ioffe and Christian Szegedy in 2015, normalizes the pre-activations of each layer to have approximately zero mean and unit variance over each mini-batch, then applies a learned affine transformation ^[17]. The Gaussian connection is partly motivational (the original paper appealed to "internal covariate shift" and the idea that activations should have a stable distribution) and partly mechanical (the normalization step rescales to $N(0, 1)$-like statistics).

The theoretical justification has shifted over the years. The original "covariate shift" story turned out to be incomplete, and follow-up work (Santurkar et al. 2018, How Does Batch Normalization Help Optimization?) argued that the real benefit comes from smoothing the loss landscape. The empirical fact that batch norm works extremely well in practice has not been seriously contested.

Variational autoencoders

VAEs (Kingma and Welling, 2013, Auto-Encoding Variational Bayes) parameterize the encoder as $q(z \mid x) = N(\mu_\varphi(x), \mathrm{diag}(\sigma_\varphi^2(x)))$ and the prior as $p(z) = N(0, I)$ ^[18]. The KL divergence between two diagonal Gaussians has a closed-form expression, which makes the training objective tractable. The reparameterization trick (sample $\varepsilon \sim N(0, I)$ and set $z = \mu + \sigma \varepsilon$) is what allows the encoder gradients to flow through the sampling step.

Without the Gaussian assumption, none of this works in closed form. You can substitute other distributions (Bernoulli for discrete latents, mixture of Gaussians for multimodal latents), but every alternative pays a tractability cost.

Gaussian processes

A Gaussian process is a stochastic process where any finite collection of values has a joint multivariate normal distribution ^[12]. Specifying a GP requires a mean function $m(x)$ and a kernel function $k(x, x')$, and posterior inference under Gaussian likelihood reduces to a linear system involving the kernel matrix.

GPs are the canonical example of nonparametric Bayesian regression and are widely used in Bayesian optimization, where the GP serves as a surrogate model for an expensive black-box function. The catch is the $O(n^3)$ scaling of the kernel-matrix inversion, which limits exact GPs to a few thousand training points; sparse approximations and inducing-point methods extend this to hundreds of thousands.

How do diffusion models use Gaussian noise?

Diffusion models, and specifically DDPMs (Ho, Jain, Abbeel 2020, Denoising Diffusion Probabilistic Models), define a forward process that gradually corrupts data $x_0$ by adding Gaussian noise according to a variance schedule $\beta_t$ ^[19]:

$$q(x_t \mid x_{t-1}) = N(x_t; \sqrt{1 - \beta_t} \, x_{t-1}, \beta_t I)$$

After $T$ steps the data becomes pure Gaussian noise. The model learns to invert this process by predicting the noise (or equivalently, the score) at each step. Both the forward process and the conditional distributions in the reverse process are Gaussian, which makes the entire framework analytically tractable. Score-based models and continuous-time formulations (Song et al. 2021) extend the Gaussian noise to a stochastic differential equation, but the Gaussian structure is the entire reason the math works out. See diffusion model for a broader overview.

Reparameterization

The reparameterization trick is one of those small ideas with outsized consequences. If you want to compute the gradient of an expectation $\mathbb{E}_{z \sim q_\varphi(z)}[f(z)]$ with respect to $\varphi$, and $q$ is the family $q_\varphi(z) = N(\mu_\varphi, \sigma_\varphi^2)$, you can rewrite the sample as $z = \mu_\varphi + \sigma_\varphi \varepsilon$ where $\varepsilon \sim N(0, I)$ ^[18]. The gradient $\nabla_\varphi$ then passes through the deterministic transformation, and you avoid the high-variance score-function estimator. Almost every modern stochastic latent variable model uses this trick somewhere, and almost every implementation relies on the fact that $N(\mu, \sigma^2)$ decomposes into $\mu + \sigma N(0, 1)$.

Stochastic gradient methods

The stochastic gradient Langevin dynamics algorithm (Welling and Teh, 2011) injects Gaussian noise into SGD updates ^[20]:

$$\theta_{t+1} = \theta_t + \frac{\varepsilon_t}{2} \nabla \log p(\theta_t \mid \text{data}) + N(0, \varepsilon_t I)$$

This turns SGD into an approximate posterior sampler, providing a bridge between optimization and Markov chain Monte Carlo. The Gaussian noise term is what gives the chain its mixing properties, and the magnitude is chosen to balance the deterministic drift against the stochastic exploration.

More generally, SGD itself is sometimes analyzed as a stochastic differential equation with Gaussian-like gradient noise, which produces useful intuitions about generalization and the implicit regularization of large learning rates.

How is the Gaussian used in differential privacy?

The Gaussian mechanism is one of the two foundational tools (alongside the Laplace mechanism) for achieving differential privacy, the formal standard for releasing statistics or training models without leaking information about any single record ^[21]. To make a query $f$ satisfy $(\varepsilon, \delta)$-differential privacy, the mechanism adds independent Gaussian noise to each output coordinate: it returns $f(D) + N(0, \sigma^2)$. Dwork and Roth's standard analysis shows that for $\varepsilon < 1$, choosing

$$\sigma \ge \sqrt{2 \ln(1.25 / \delta)} \, \Delta_2 / \varepsilon$$

is sufficient, where $\Delta_2$ is the $\ell_2$ sensitivity of $f$ (the most any single record can change the output, measured in Euclidean norm) ^[21]. The noise scale grows with sensitivity and shrinks with the privacy budget $\varepsilon$, which makes the privacy-accuracy tradeoff explicit. This is exactly the mechanism behind differentially private stochastic gradient descent (DP-SGD), where clipped per-example gradients are perturbed with Gaussian noise before the optimizer step, and it is why the Gaussian shows up at the center of privacy-preserving machine learning.

Random Fourier features

Ali Rahimi and Ben Recht's 2007 paper Random Features for Large-Scale Kernel Machines showed that for shift-invariant kernels (in particular the Gaussian RBF kernel), the kernel can be approximated by an inner product of random feature maps where the random directions are sampled from a Gaussian distribution ^[22]. This converts kernel methods, which scale $O(n^2)$ in memory and $O(n^3)$ in time, into linear methods that scale $O(nD)$ for $D$ random features.

The Gaussian sampling here is a direct consequence of Bochner's theorem applied to the RBF kernel: the spectral measure of an RBF kernel is itself a Gaussian, so to sample random frequencies that approximate the kernel you sample from that Gaussian.

When should you not assume normality?

The most common failure mode in applied statistics is assuming normality where it does not hold. A few areas where this has caused real problems:

Finance. Asset returns are notoriously not normal. Daily log returns of equities have heavy tails, with extreme moves (the 1987 crash, the 2008 collapse, the 2020 March drawdown) being orders of magnitude more frequent than a Gaussian model predicts. Long-Term Capital Management famously blew up in 1998 in part because their risk models assumed Gaussian returns and the 1998 Russian default produced moves that, under their model, should have happened roughly once in the lifetime of the universe. The standard fix is to use Student's t-distribution, generalized hyperbolic distributions, or stable distributions with infinite variance.

Network traffic and outliers in computing. Inter-arrival times of network packets, file sizes on the web, sizes of scientific computations: these all tend to follow heavy-tailed distributions (Pareto, log-normal, power-law) rather than Gaussians. Modeling them as Gaussian leads to wildly wrong capacity-planning decisions.

High-dimensional ML data. Pixel values in images, embedding norms, gradient magnitudes during training, all of these often have non-Gaussian distributions, and assuming Gaussianity in the wrong place can produce surprising failures. Gradient clipping in deep learning exists partly because gradient magnitudes have heavier tails than a Gaussian model would suggest.

Bayesian inference with heavy-tailed priors. Robust priors like Cauchy or Student-t are often preferable to Gaussian priors when you genuinely do not have strong prior information, because they put more weight on extreme values being possible.

The right response to non-normality depends on what you are doing. For point estimation with reasonable sample sizes, the central limit theorem often saves you. For tail-risk estimation, value-at-risk calculations, or anywhere extreme events matter, Gaussian assumptions can be catastrophic and you should reach for a heavier-tailed distribution explicitly.

Generalizations and related distributions

Scroll sideways for more →

Most of these are connected to the normal distribution by simple transformations (Cholesky, exponentiation, summing of squares), which is one reason the Gaussian sits at the center of so much classical statistics.

Why does machine learning rely so heavily on Gaussians?

At the risk of belaboring the point, here is a short list of why Gaussians dominate machine learning specifically:

1. Tractability. Gaussian computations are closed-form. Marginals, conditionals, products, and sums of Gaussians are all Gaussian, and the parameters update in simple ways ^[11].
2. MLE = OLS. For linear regression with Gaussian noise, the maximum-likelihood estimate equals the ordinary least squares estimate. Every quadratic loss function in ML is implicitly a Gaussian likelihood ^[3].
3. Central limit theorem. Sums and averages of many small effects tend to be approximately Gaussian, which justifies Gaussian noise models for residuals and aggregates ^[12].
4. Maximum entropy. The Gaussian is the most uncertain distribution consistent with a known mean and variance, which is a principled default when you do not want to commit to more structure ^[4].
5. Rotational invariance. The isotropic Gaussian $N(0, \sigma^2 I)$ is invariant under rotations of the coordinate system, which is the right symmetry property for many high-dimensional problems where there is no preferred basis.
6. Reparameterizable. $N(\mu, \sigma^2)$ decomposes as $\mu + \sigma N(0, 1)$, enabling the reparameterization trick and gradient-based learning of distributional parameters ^[18].
7. Conjugate to itself. With known variance, a Gaussian prior plus a Gaussian likelihood gives a Gaussian posterior, with a one-line update rule for the parameters ^[11].

None of this means you should always use a Gaussian. It does mean that when you can use one without doing real damage, doing so will save you a lot of work.

See also

• Central limit theorem
• Gaussian process
• Variational autoencoder
• Diffusion models
• Batch normalization
• Linear regression
• Gaussian mixture model
• Bayesian inference
• Maximum likelihood
• Box-Muller transform
• Reparameterization trick
• Weight initialization
• Xavier initialization
• Mean squared error
• Markov chain Monte Carlo
• Information theory
• Entropy
• t-distribution
• Standard deviation
• Carl Friedrich Gauss
• Pierre-Simon Laplace

References

1. de Moivre, A. (1733). *Approximatio ad Summam Terminorum Binomii (a + b)^n in Seriem expansi*. Privately printed supplement, later incorporated into the 1738 second edition of *The Doctrine of Chances*. ↩
2. Laplace, P.-S. (1812). *Théorie analytique des probabilités*. Paris: Courcier. ↩
3. Gauss, C. F. (1809). *Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium*. Hamburg: Perthes & Besser. ↩
4. Quetelet, A. (1835). *Sur l'homme et le développement de ses facultés, ou Essai de physique sociale*. Paris: Bachelier. (Maximum-entropy and least-informative-distribution framing per Cover and Thomas, *Elements of Information Theory*.) ↩
5. "68-95-99.7 rule." *Wikipedia, The Free Encyclopedia*. Exact values 0.6826894921, 0.9544997361, 0.9973002039. https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule ↩
6. Stigler, S. M. (1986). *The History of Statistics: The Measurement of Uncertainty before 1900*. Harvard University Press.
7. Hald, A. (2007). "De Moivre's Normal Approximation to the Binomial, 1733, and Its Generalization." In *A History of Parametric Statistical Inference from Bernoulli to Fisher*. Springer. ↩
8. Stigler, S. M. (1999). *Statistics on the Table: The History of Statistical Concepts and Methods*. Harvard University Press. ↩
9. Pearson, K. (1893). "Asymmetrical Frequency Curves." *Nature* 48, 615-616. ↩
10. Pearson, K. (1920). "Notes on the History of Correlation." *Biometrika* 13(1), 25-45. ↩
11. Bishop, C. M. (2006). *Pattern Recognition and Machine Learning*. Springer. ↩
12. Murphy, K. P. (2012). *Machine Learning: A Probabilistic Perspective*. MIT Press. ↩
13. Box, G. E. P., and Muller, M. E. (1958). "A Note on the Generation of Random Normal Deviates." *Annals of Mathematical Statistics* 29(2), 610-611. ↩
14. Marsaglia, G., and Tsang, W. W. (2000). "The Ziggurat Method for Generating Random Variables." *Journal of Statistical Software* 5(8), 1-7. ↩
15. Glorot, X., and Bengio, Y. (2010). "Understanding the Difficulty of Training Deep Feedforward Neural Networks." *Proceedings of AISTATS* 9, 249-256. ↩
16. He, K., Zhang, X., Ren, S., and Sun, J. (2015). "Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification." *Proceedings of ICCV*. arXiv:1502.01852. ↩
17. Ioffe, S., and Szegedy, C. (2015). "Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift." *Proceedings of ICML*. arXiv:1502.03167. ↩
18. Kingma, D. P., and Welling, M. (2013). "Auto-Encoding Variational Bayes." arXiv:1312.6114. ↩
19. Ho, J., Jain, A., and Abbeel, P. (2020). "Denoising Diffusion Probabilistic Models." *Advances in Neural Information Processing Systems* 33. arXiv:2006.11239. ↩
20. Welling, M., and Teh, Y. W. (2011). "Bayesian Learning via Stochastic Gradient Langevin Dynamics." *Proceedings of ICML*. ↩
21. Dwork, C., and Roth, A. (2014). *The Algorithmic Foundations of Differential Privacy*. Foundations and Trends in Theoretical Computer Science 9(3-4). (Gaussian mechanism and the σ ≥ √(2 ln(1.25/δ)) · Δ_2 / ε bound; mechanism introduced in Dwork et al. 2006.) ↩
22. Rahimi, A., and Recht, B. (2007). "Random Features for Large-Scale Kernel Machines." *Advances in Neural Information Processing Systems* 20. ↩
23. Casella, G., and Berger, R. L. (2002). *Statistical Inference*, 2nd edition. Duxbury.
24. Wikipedia contributors. "Normal distribution." *Wikipedia, The Free Encyclopedia*. https://en.wikipedia.org/wiki/Normal_distribution