Why Sums Become Gaussian
10 min read
Measurement errors add up, and sample averages combine many random observations. Gaussian distributions appear repeatedly in both settings, for two different reasons. Adding independent Gaussians gives another Gaussian exactly; averaging many other kinds of independent observations can give something approximately Gaussian.
The geometric proof, which I learned from Grant Sanderson’s 3Blue1Brown lessons (one, two, and three), explains why addition preserves the Gaussian shape.
Sums Collect Anti-Diagonals
How does the density of a sum come from the densities of the two summands? Start with two dice. Their \(6 \times 6 = 36\) pairs of faces are equally likely and can be arranged in a square grid. The sum depends only on the anti-diagonal containing a pair: six cells total seven, while only one totals two and one totals twelve. Counting the cells on each anti-diagonal gives the triangle below, so adding two uniform distributions already produces a peak in the middle.
The continuous case is the same picture with the cells shrunk to points. For independent draws \(x\) and \(y\), the joint density is the surface \(p_x(t)\,p_y(u)\), and all pairs with sum \(s\) lie on the line \(t + u = s\). Integrating the joint density along that line gives the density of the sum at \(s\). Whether that slice has the same shape as a marginal depends on the shape of the joint surface. For independent Gaussians the surface is round, which is the whole proof; cumulants then explain why other standardized sums drift toward the same shape.
Circles and Diamonds
What is special about the joint density of two independent Gaussians? For two independent standard Gaussians \(z_1\) and \(z_2\), the joint density is a product, and the exponents merge:
\[ \frac{1}{\sqrt{2\pi}}\, e^{-t^2/2} \cdot \frac{1}{\sqrt{2\pi}}\, e^{-u^2/2} = \frac{1}{2\pi}\, e^{-(t^2 + u^2)/2}. \]
The pair \((t, u)\) enters only through \(t^2 + u^2\), its squared distance from the origin. Every contour is therefore a circle, and rotating the plane leaves the surface unchanged. This roundness comes directly from independence: multiplying the two Gaussian densities adds their squared exponents, and \(t^2 + u^2\) is a squared radius.
Most product densities are not round. Two independent Laplace draws, each with density \(\frac{1}{2}e^{-|t|}\), have joint density \(\frac{1}{4}e^{-(|t| + |u|)}\), whose contours, the sets on which \(|t| + |u|\) is constant, are diamonds. The plot slices both surfaces vertically (dotted), to obtain a marginal density, and along anti-diagonals (dashed), to obtain the density of the sum. For the Laplace density, the two directions produce different profiles, and the sum of two Laplace draws is not Laplace. For the Gaussian, one set of slices is just a rotation of the other.
This roundness characterizes the Gaussian: Herschel, and later Maxwell, showed that the joint density of two independent identically distributed draws has circular contours only when the draws are Gaussian [Maxwell, Phil. Mag. 1860].
A Rotation Instead of an Integral
How does roundness make the sum Gaussian without integrating? We can replace the convolution integral with a rotation, and what comes out is the usual closure rule.
Claim: for independent summands with means \(\mu_1, \mu_2\) and variances \(\sigma_1^2, \sigma_2^2\),
\[ \mathcal{N}(\mu_1, \sigma_1^2) + \mathcal{N}(\mu_2, \sigma_2^2) = \mathcal{N}\big(\mu_1 + \mu_2,\; \sigma_1^2 + \sigma_2^2\big). \]
Proof
Start with the standard Gaussians \(z_1\) and \(z_2\), and rotate the plane by \(45°\) so that the diagonal becomes an axis:
\[ w_1 = \frac{z_1 + z_2}{\sqrt{2}}, \qquad w_2 = \frac{z_1 - z_2}{\sqrt{2}}. \]
(The \(\sqrt{2}\) makes this a rotation rather than a stretch: \(w_1^2 + w_2^2 = z_1^2 + z_2^2\).) Rotation preserves both radius and area, so \((w_1, w_2)\) has the same joint density as \((z_1, z_2)\). Thus \(w_1\) and \(w_2\) are again independent standard Gaussians. The anti-diagonal slices that define the sum are now the vertical slices that define a marginal, giving
\[ z_1 + z_2 = \sqrt{2}\, w_1 \sim \mathcal{N}(0, 2). \]
For unequal standard deviations, rotate the first axis to point along \((\sigma_1, \sigma_2)\); then \(\sigma_1 z_1 + \sigma_2 z_2 \sim \mathcal{N}(0, \sigma_1^2 + \sigma_2^2)\). Adding deterministic means shifts the center, which gives the closure rule.
The sum is Gaussian because the round joint density cannot distinguish a diagonal from an axis. Means and variances add for any independent sum; what is special here is that the Gaussian shape also survives.
The Fixed Point
Why do sums of non-Gaussian draws also approach a Gaussian? Closure also makes the Gaussian a fixed point of averaging: add an independent copy and rescale to restore the variance, and \((z_1 + z_2)/\sqrt{2}\) is standard Gaussian again. The central limit theorem says that, under suitable conditions, repeated averaging draws other distributions toward this fixed point. Take independent draws \(x_1, \dots, x_n\) from a distribution with mean zero and finite variance \(\sigma^2\), and standardize their sum:
\[ y_n = \frac{x_1 + \cdots + x_n}{\sigma \sqrt{n}}. \]
As \(n\) grows, the distribution of \(y_n\) approaches the standard Gaussian. The panels below show this for a deliberately awkward starting point: a lopsided, two-bump density. Pairwise convolution produces three bumps at \(n = 2\); by \(n = 4\) they have become ripples, and by \(n = 16\) the density is visually indistinguishable from the Gaussian.
The plots show the convergence but not why the Gaussian is the limit. Cumulants provide one explanation: after standardization, repeated addition suppresses higher-order shape features while preserving variance.
The Cumulant Decay Rates
What gets washed out with each addition, and how fast? We can see the first two shape changes without defining cumulants. Recall that each \(x_i\) has mean zero and variance \(\sigma^2\), and that \(y_n\) is their sum divided by \(\sigma\sqrt{n}\). Because \(y_n\) has mean zero and variance one, its skewness is simply \(\mathbb{E}[y_n^3]\).
When we expand the cube of the sum, every mixed term has mean zero. For example, independence gives \(\mathbb{E}[x_i^2x_k] = \mathbb{E}[x_i^2]\mathbb{E}[x_k] = 0\) when \(i \ne k\). Only the \(n\) pure cubes remain, so the skewness is:
\[\begin{align} \mathbb{E}[y_n^3] &= \frac{\mathbb{E}[(x_1 + \cdots + x_n)^3]}{\sigma^3 n^{3/2}} \\ &= \frac{n\,\mathbb{E}[x_1^3]}{\sigma^3 n^{3/2}} \\ &= \frac{\mathbb{E}[x_1^3]}{\sigma^3}\,\frac{1}{\sqrt{n}}. \end{align}\]
The first line substitutes the standardized sum, the second keeps the \(n\) pure cubes, and the third cancels the powers of \(n\). Thus adding a hundred draws and restoring variance reduces the original skewness by a factor of ten.
The fourth power has one extra kind of surviving term. Besides the \(n\) pure fourth powers, products such as \(x_i^2x_k^2\) remain because neither factor has mean zero. There are \(6\binom{n}{2}\) of these pairings, which gives:
\[\begin{align} \mathbb{E}[(x_1 + \cdots + x_n)^4] &= n\,\mathbb{E}[x_1^4] + 6\binom{n}{2}\sigma^4 \\ &= n\,\mathbb{E}[x_1^4] + 3n(n-1)\sigma^4. \end{align}\]
The first line lists the surviving terms, and the second simplifies the number of pairings.
For a variance-one variable, excess kurtosis is its fourth moment minus \(3\). Dividing the last display by the fourth power of \(\sigma\sqrt{n}\) shows what remains:
\[\begin{align} \mathbb{E}[y_n^4] - 3 &= \frac{\mathbb{E}[x_1^4]}{\sigma^4 n} + \frac{3(n-1)}{n} - 3 \\ &= \left(\frac{\mathbb{E}[x_1^4]}{\sigma^4} - 3\right)\frac{1}{n}. \end{align}\]
The pairings produce the Gaussian value \(3\). Any excess above or below it is divided by \(n\), so a hundred draws reduce the original excess kurtosis by a factor of one hundred.
Cumulants extend these two calculations to every order. The first cumulant is the mean and the second is the variance. The next two are the quantities we just tracked:
\[ \kappa_3 = \mathbb{E}\big[(x-\mathbb{E}[x])^3\big], \qquad \kappa_4 = \mathbb{E}\big[(x-\mathbb{E}[x])^4\big] - 3\mathrm{Var}(x)^2. \]
After division by \(\sigma^3\) and \(\sigma^4\), these are skewness and excess kurtosis. Higher cumulants subtract combinations of lower-order moments, generalizing the subtraction of \(3\mathrm{Var}(x)^2\) at order four.
Their formal definition packages this cancellation into two useful rules. Provided the expectation is finite for \(t\) near zero, define the cumulant generating function of \(x\) by:
\[ K_x(t) = \ln \mathbb{E}\big[e^{tx}\big] = \kappa_1 t + \kappa_2\frac{t^2}{2!} + \kappa_3\frac{t^3}{3!} + \cdots. \]
The logarithm is there because independence turns an exponential average of a sum into a product. For independent \(x\) and \(y\), it turns that product back into addition:
\[\begin{align} K_{x+y}(t) &= \ln \left(\mathbb{E}\big[e^{tx}\big]\mathbb{E}\big[e^{ty}\big]\right) \\ &= K_x(t) + K_y(t). \end{align}\]
Scaling \(x\) by a constant \(c\) instead replaces \(t\) by \(ct\):
\[ K_{cx}(t) = K_x(ct). \]
Comparing the Taylor coefficients gives the two rules we need:
\[ \kappa_j(x+y) = \kappa_j(x) + \kappa_j(y), \qquad \kappa_j(cx) = c^j\kappa_j(x). \]
Now apply addition and scaling to the standardized sum:
\[\begin{align} \kappa_j(y_n) &= \frac{\kappa_j(x_1 + \cdots + x_n)}{(\sigma\sqrt{n})^{\,j}} \\ &= \frac{n\,\kappa_j(x_1)}{\sigma^j n^{j/2}} \\ &= \frac{\kappa_j(x_1)}{\sigma^j}\,n^{1-j/2}. \end{align}\]
The lines apply the scale rule, add \(n\) identical cumulants, and then collect the powers of \(n\). Addition supplies \(n\) copies of every cumulant, while restoring the variance divides an order-\(j\) cumulant by \(n^{j/2}\). At order two the factors cancel, so the variance stays one. At every order above two, the division is stronger: skewness falls as \(n^{-1/2}\), excess kurtosis falls as \(n^{-1}\), and the higher orders fall faster.
The Gaussian is what remains because its cumulant generating function stops after the variance term. For a standard Gaussian \(z\), we can compute it directly:
\[\begin{align} K_z(t) &= \ln \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}e^{-\tau^2/2}e^{t\tau}\,d\tau \\ &= \ln \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}e^{-(\tau^2-2t\tau)/2}\,d\tau \\ &= \ln \left(e^{t^2/2}\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}e^{-(\tau-t)^2/2}\,d\tau\right) \\ &= \frac{t^2}{2}. \end{align}\]
The lines write the expectation as an integral, combine the exponents, and complete the square. The remaining integral is one because it is the total density of a Gaussian centered at \(t\). There are no third- or higher-order terms, so every higher Gaussian cumulant is zero. The plot shows the exact decay rates for the two-bump density used above.
(This cumulant argument uses a moment generating function that exists near zero. A finite-variance distribution need not have one, so general proofs work with the characteristic function \(\mathbb{E}[e^{\mathrm{i}tx}]\), which always exists because \(|e^{\mathrm{i}tx}| = 1\), or with Lindeberg’s argument that swaps the summands for Gaussians one at a time [Lindeberg, Math. Z. ’22].)
The Closure at Work
When an application invokes the Gaussian, which of the two facts is it using? The limit theorem and the closure rule serve different roles: the first makes many standardized sums approximately Gaussian, and the second makes sums of Gaussians exactly Gaussian. An instrument’s error often combines many small independent disturbances, which motivates the Gaussian error model and the Gaussian quantile in the sample-mean interval \(\bar{x} \pm 1.96\,\sigma/\sqrt{n}\). In differential privacy, privacy loss from Gaussian mechanisms adds across releases, allowing a long pipeline to be analyzed through a single Gaussian tail bound. In randomized linear algebra, a matrix with independent Gaussian entries has the same distribution in every orthogonal basis, the higher-dimensional version of roundness used by the sketches in private regression via iterative Hessian mixing.
For independent sums, means and variances always add. For Gaussians, those two totals determine the entire distribution; the central limit theorem explains why they also come to dominate many non-Gaussian sums.