Why Sums Become Gaussian
12 min read
Sums of independent random quantities are everywhere: measurement errors accumulate, sample means average, noise composes across the steps of an algorithm. Two facts about the Gaussian govern all of them. Sums of independent Gaussians are exactly Gaussian. And sums of almost anything else become Gaussian. The textbook proves the first with a convolution integral and a completed square, which certifies the identity without ever showing why the Gaussian is the shape that survives. There is a proof that shows why, which I learned from Grant Sanderson’s 3Blue1Brown lessons [Sanderson, ’23]:
Everything on this side of the story is making those sentences true: summing is slicing, and round means no direction is special. The second fact starts from the first: a shape addition cannot change is a shape sums can rest in, and the real work is showing they are pulled there.
Sums Collect Anti-Diagonals
Start with two dice. The pair of faces takes each of its \(6 \times 6 = 36\) values with equal probability, one cell each in a square grid. The sum does not care which cell you land in, only which anti-diagonal the cell lies on: the rolls with total seven are exactly the six cells of the main anti-diagonal, while a total of two or twelve claims a single corner cell. Counting cells per anti-diagonal gives the triangle in the plot: adding two flat distributions has already built a peak in the middle.
The continuous version is the same picture with the cells shrunk to points [Sanderson, ’23]. For independent draws \(x\) and \(y\), the joint density is the surface \(p_x(t)\,p_y(u)\) over the plane, and the pairs that total \(s\) form the line \(t + u = s\), a continuous anti-diagonal. The density of the sum at \(s\) collects everything on that line: it is the integral of the joint density along it. Summing is slicing: sweep the line across the plane and record what each slice holds. So everything a sum can do is written in the shape of one surface.
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\), the squared distance from the origin, so every contour is a circle: the surface is a perfectly round bell, unchanged by any rotation of the plane. The cause is mechanical. Multiplying densities adds exponents, the Gaussian’s exponent is a square, and a sum of squares is a squared radius: independence hands the Gaussian the Pythagorean theorem.
Roundness is rarer than it looks. Two independent Laplace draws, each with density \(\frac{1}{2}e^{-|t|}\), have joint density \(\frac{1}{4}e^{-(|t| + |u|)}\), whose contours \(|t| + |u| = c\) are diamonds. In the plot, both surfaces are sliced two ways: vertically (dotted), which reads off the density of the first draw alone, and along the anti-diagonals (dashed), which reads off the density of the sum. On the diamonds the two slicings meet genuinely different profiles, and indeed the sum of two Laplace draws is no longer Laplace. On the circles the two slicings are the same slices, rotated.
Roundness is in fact the Gaussian’s signature: Herschel, and later Maxwell, showed that a pair of independent identically distributed draws has circular contours only when the draws are Gaussian [Maxwell, Phil. Mag. 1860].
A Rotation Instead of an Integral
With the surface in hand, we can trade the convolution integral for a rotation. 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}\) is what makes this a rotation rather than a stretch: \(w_1^2 + w_2^2 = z_1^2 + z_2^2\).) A rotation changes neither the surface, whose height depends only on the radius, nor areas, so the rotated pair \((w_1, w_2)\) has exactly the joint density that \((z_1, z_2)\) had: \(w_1\) and \(w_2\) are again independent standard Gaussians. The anti-diagonal slices that build the sum have become the vertical slices that build a marginal, and the marginal of a product density is just its first factor. Reading off the sum, a standard Gaussian stretched by \(\sqrt{2}\):
\[ z_1 + z_2 = \sqrt{2}\, w_1 \sim \mathcal{N}(0, 2). \]
No convolution integral and no completed square: the sum is Gaussian because a round surface cannot tell the diagonal from an axis. Unequal spreads only tilt the angle, since the same argument with the rotation whose first axis points along \((\sigma_1, \sigma_2)\) shows that \(\sigma_1 z_1 + \sigma_2 z_2 \sim \mathcal{N}(0, \sigma_1^2 + \sigma_2^2)\), and deterministic means just shift the center. That is the closure rule in the form it gets used:
\[ \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) \qquad \text{(independent summands)}. \]
Means add and variances add, which any independent sum must obey; the content of the identity is the shape surviving.
The Fixed Point
Closure has a consequence one step further on. The operation “add an independent copy, then rescale to restore the variance” maps distributions to distributions, and the Gaussian passes through untouched: \((z_1 + z_2)/\sqrt{2}\) is standard Gaussian again. The Gaussian is a fixed point of averaging. The central limit theorem says this fixed point attracts everything. Take independent draws \(x_1, \dots, x_n\) from any distribution with variance \(\sigma^2\) (subtract the mean once so it is zero) and standardize the sum:
\[ y_n = \frac{x_1 + \cdots + x_n}{\sigma \sqrt{n}}. \]
The distribution of \(y_n\) then approaches the standard Gaussian, whatever the start, so long as the variance is finite. Sanderson’s lesson plays this as an animation, convolving a density with itself over and over and watching the bumps melt [Sanderson, ’23]; the panels below are four frames of that animation for a deliberately awkward start, a lopsided two-bump density. At \(n = 2\) the density has three bumps, since the bumps convolve pairwise; by \(n = 4\) the bumps survive only as ripples; at \(n = 16\) the density sits on the Gaussian to within the width of the line.
The panels show it happening. What they do not show is why the destination is the Gaussian and never anything else, and that has a headline of its own:
Addition Keeps Only Two Numbers
To make “everything else” precise, we need coordinates for the shape of a distribution, chosen so that adding independent draws does something simple to them. Exponentials are the natural tool, since an exponential turns a sum into a product. For a random variable \(x\) and a fixed number \(t\), consider \(\mathbb{E}[e^{tx}]\), a function of \(t\) called the moment generating function. Expanding \(e^{tx} = 1 + tx + \frac{(tx)^2}{2!} + \cdots\) inside the expectation and averaging term by term shows what the function stores:
\[ \mathbb{E}\big[e^{tx}\big] = 1 + t\,\mathbb{E}[x] + \frac{t^2}{2!}\,\mathbb{E}[x^2] + \frac{t^3}{3!}\,\mathbb{E}[x^3] + \cdots. \]
One function packages every moment of the distribution. (The standing assumption that \(\mathbb{E}[e^{tx}]\) is finite near \(t = 0\) is exactly what licenses the term-by-term averaging.)
The reason to prefer this function to the density itself is what it does to sums. For independent \(x\) and \(y\), the exponential of a sum splits into a product, and the expectation of a product of independent factors is the product of the expectations:
\[ \mathbb{E}\big[e^{t(x+y)}\big] = \mathbb{E}\big[e^{tx}\, e^{ty}\big] = \mathbb{E}\big[e^{tx}\big]\, \mathbb{E}\big[e^{ty}\big]. \]
Adding random variables multiplies their generating functions, and the logarithm turns that multiplication into addition:
\[ \ln \mathbb{E}\big[e^{t(x+y)}\big] = \ln \mathbb{E}\big[e^{tx}\big] + \ln \mathbb{E}\big[e^{ty}\big]. \]
Convolution of densities has become addition of functions, and this one line is the entire reason the transform exists.
The cumulants \(\kappa_j\) of \(x\) are the Taylor coefficients of that logarithm:
\[ \ln \mathbb{E}\big[e^{t x}\big] = \kappa_1 t + \kappa_2 \frac{t^2}{2!} + \kappa_3 \frac{t^3}{3!} + \cdots \]
(no constant term, because at \(t = 0\) the expectation is \(\mathbb{E}[1] = 1\) and its logarithm is zero). The first two coefficients are old friends. Substitute the moment expansion into \(\ln(1 + u) = u - \frac{u^2}{2} + O(u^3)\), taking \(u = t\,\mathbb{E}[x] + \frac{t^2}{2}\,\mathbb{E}[x^2] + O(t^3)\), and keep everything up to \(t^2\):
\[\begin{align} \ln \mathbb{E}\big[e^{tx}\big] &= t\,\mathbb{E}[x] + \frac{t^2}{2}\,\mathbb{E}[x^2] - \frac{1}{2}\big(t\,\mathbb{E}[x]\big)^2 + O(t^3) \\ &= t\,\mathbb{E}[x] + \frac{t^2}{2}\Big(\mathbb{E}[x^2] - \mathbb{E}[x]^2\Big) + O(t^3), \end{align}\]
where the first line keeps the two terms of \(u\) and the square of its leading term (every other product is \(O(t^3)\)), and the second collects the powers of \(t\). Matching coefficients against the definition: \(\kappa_1\) is the mean and \(\kappa_2\) is the variance. Two more rounds of the same expansion give the shape features: \(\kappa_3 = \mathbb{E}[(x - \kappa_1)^3]\), the third central moment, and \(\kappa_4 = \mathbb{E}[(x - \kappa_1)^4] - 3\kappa_2^2\). Dividing out scale gives the classical summaries, the skewness \(\kappa_3/\sigma^3\) (asymmetry) and the excess kurtosis \(\kappa_4/\sigma^4\) (tail weight). The \(-3\kappa_2^2\) correction in \(\kappa_4\) looks unmotivated here; in a moment it will be exactly what the Gaussian requires.
Two rules now fall out of the definition.
Sums. The logarithms add for independent summands, and the Taylor coefficients of a sum of functions are the sums of the coefficients, so every cumulant adds:
\[ \kappa_j(x + y) = \kappa_j(x) + \kappa_j(y). \]
Scale. Replacing \(x\) by \(c\,x\) substitutes \(ct\) for \(t\) in the generating function, since \(\mathbb{E}[e^{t(cx)}] = \mathbb{E}[e^{(ct)x}]\), and re-expanding in powers of \(t\), the \(j\)-th power of \(t\) drags along \(c^j\):
\[ \kappa_j(c\, x) = c^j\, \kappa_j(x). \]
High-order cumulants carry many powers of the scale, which is about to matter.
Where does the Gaussian sit in these coordinates? One completed square answers:
\[\begin{align} \mathbb{E}\big[e^{t z}\big] &= \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}\, e^{t\tau - \tau^2/2}\, d\tau \\ &= e^{t^2/2} \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}\, e^{-(\tau - t)^2/2}\, d\tau \\ &= e^{t^2/2}, \end{align}\]
where the second line rewrites the exponent as \(t\tau - \tau^2/2 = \frac{t^2}{2} - \frac{(\tau - t)^2}{2}\) (expand \((\tau - t)^2\) to check) and pulls the constant factor out of the integral, and the third recognizes the integrand as a Gaussian density recentered at \(t\), which integrates to one like any density. So for the standard Gaussian:
\[ \ln \mathbb{E}\big[e^{t z}\big] = \frac{t^2}{2}. \]
The series does not merely start with a quadratic; it is a quadratic. Matching coefficients one more time: \(\kappa_2 = 1\), and every cumulant from \(\kappa_3\) on is exactly zero. (This is what the \(-3\kappa_2^2\) correction was for: the Gaussian’s fourth central moment is \(3\sigma^4\), and the correction cancels it precisely.) A general Gaussian is \(\mu + \sigma z\), and \(\mathbb{E}[e^{t(\mu + \sigma z)}] = e^{\mu t}\, \mathbb{E}[e^{(\sigma t) z}]\) gives the logarithm \(\mu t + \sigma^2 t^2/2\). A mean, a variance, and nothing else.
Now push the standardized sum \(y_n\) through the two rules, one move per line:
\[\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\sqrt{n})^{\,j}} \\ &= \frac{\kappa_j(x_1)}{\sigma^{\,j}} \cdot \frac{n}{n^{\,j/2}} \\ &= \frac{\kappa_j(x_1)}{\sigma^{\,j}}\; n^{\,1 - j/2}, \end{align}\]
where the first line applies the scale rule with \(c = 1/(\sigma\sqrt{n})\), the second applies the sum rule to \(n\) summands with identical cumulants, the third splits \((\sigma\sqrt{n})^j = \sigma^j\, n^{j/2}\), and the fourth collects the powers of \(n\). The competition on the third line is the whole theorem: the numerator grows like \(n\), one deposit per summand, while the rescaling raises \(\sqrt{n}\) to the \(j\)-th power in the denominator. At \(j = 2\) the two exactly tie. For \(j \ge 3\) the denominator wins.
Read the exponent \(1 - j/2\) across \(j\). The mean stays zero, since we centered. The variance gets exponent zero, exactly \(1\) for every \(n\): no surprise, because we divided by the standard deviation on purpose, and the algebra can only confirm the bookkeeping. The skewness shrinks like \(n^{-1/2}\), the excess kurtosis like \(n^{-1}\), the fifth cumulant faster still: sum a hundred draws and the skewness has dropped \(10\times\), the excess kurtosis \(100\times\). Every feature past the first two drains away at a known rate, and the distribution left holding a bare mean and variance is, by the completed square above, the Gaussian.
A worked check makes the machine concrete. The exponential distribution with rate one has \(\mathbb{E}[e^{tx}] = \frac{1}{1-t}\) for \(t < 1\), and its logarithm has a series every calculus course computes:
\[ \ln \frac{1}{1 - t} = t + \frac{t^2}{2} + \frac{t^3}{3} + \cdots, \qquad \text{so} \qquad \kappa_j = \frac{j!}{j} = (j - 1)!, \]
matching the coefficient of \(t^j\) against \(\kappa_j/j!\). So the mean and variance are both \(1\), the skewness is \(\kappa_3/\sigma^3 = 2\), and the excess kurtosis is \(6\), the textbook values; after summing \(n\) draws, the standardized skewness is exactly \(2/\sqrt{n}\). In the plot, the same drain runs on the two-bump density from the panels above: its two shape features fall on exact power laws.
(One repair for full rigor: heavy-tailed distributions can have finite variance and no moment generating function, so the standard proof runs the same drain on \(\mathbb{E}[e^{\mathrm{i}tx}]\), which is finite for every distribution because \(|e^{\mathrm{i}tx}| = 1\); finite variance alone suffices [Lindeberg, Math. Z. ’22].)
The Closure at Work
The two facts work as a pair: the limit theorem delivers sums to the neighborhood of the Gaussian, and closure keeps them there, exactly. An instrument’s error accumulates many small independent disturbances, which is why the Gaussian error model is the default and why the sample mean’s interval \(\bar{x} \pm 1.96\,\sigma/\sqrt{n}\) ends in a Gaussian quantile. In differential privacy, the evidence leaked by each noised release is Gaussian-distributed and the leaks of successive releases add, so the privacy of a thousand-release pipeline is one Gaussian tail bound rather than a thousand compounding ones. In randomized linear algebra, a matrix of independent Gaussian entries looks the same in every orthogonal basis (roundness, run in \(n\) dimensions), which is exactly what the sketches inside private regression via iterative Hessian mixing rely on.
Means add and variances add for any independent sum; the Gaussian is the distribution that adds nothing else, because it has nothing else. Model accumulated noise as Gaussian and the whole distribution rides on two running totals, and the central limit theorem is the reason so much noise arrives already in that form.