Why Sums Become Gaussian
8 min read
The sum of two fair dice is much more likely to be seven than two. Addition changes the shape of the distribution, even though both dice were uniform. Independent Gaussians behave differently: their sum has the same bell shape, just shifted or rescaled.
A beautiful geometric proof explains why, without having to complete the square in an integral. I learned it from Grant Sanderson’s 3Blue1Brown lessons (one, two, and three). The same geometry gives us a starting point for understanding why many other sums approach a Gaussian.
Sums Collect Anti-Diagonals
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\). Parameterize that line by its first coordinate \(t\), so the second is \(s-t\). The density of the sum is the convolution:
\[ p_{x+y}(s)=\int_{-\infty}^{\infty}p_x(t)\,p_y(s-t)\,dt. \]
Whether that slice has the same shape as a marginal depends on the shape of the joint surface. For independent standard Gaussians, the surface is round, so changing the direction of a slice changes only its scale.
Circles and Diamonds
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
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
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, positive 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
We can see the first two shape changes without defining cumulants. For these calculations, assume the third and fourth moments are finite; the central limit theorem above needs only a finite, positive variance. 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.
The generating-function argument for all orders
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
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.