Thought Experiment that gave us the most famous probability distribution function ever.

The classic Herschel-Maxwell derivation (originally conceptualized by John Herschel in 1850 and formalized by James Clerk Maxwell in 1860). It is one of the most elegant and intuitive proofs in all of statistics because it builds the entire bell curve out of just two logical assumptions.

Here is how Herschel first built that argument from scratch.

The Setup

Imagine you are throwing a dart at a 2D plane, aiming for the bullseye at the origin (0, 0). The coordinates of where the dart actually lands are $(x, y)$.

To figure out the probability distribution of where the darts land, we only need to make two reasonable, common-sense postulates.

The Two Golden Assumptions

1. Orthogonal Independence Your error in the horizontal direction has absolutely no effect on your error in the vertical direction. They are completely independent. Because of this, the joint probability of landing at a specific $(x, y)$ coordinate is simply the product of their individual probabilities:

$P(x, y) = p(x)p(y)$

2. Rotational Symmetry (Isotropy) The probability of the dart landing at any point depends only on its total distance from the bullseye, not the angle. A dart 2 inches strictly right of the bullseye has the exact same probability of occurring as a dart 2 inches strictly above it. Therefore, the joint probability can also be written as a function of the radius $r$:

$P(x, y) = g(r) = g(\sqrt{x^2 + y^2})$

---

Building the Argument

Step 1: Equate the two perspectives Because both of our assumptions describe the exact same joint probability $P(x, y)$, we can set them equal to each other:

$p(x)p(y) = g(\sqrt{x^2 + y^2})$

Step 2: Logarithmic transformation To make the math easier to manipulate, take the natural log of both sides.

$\ln(p(x)) + \ln(p(y)) = \ln(g(\sqrt{x^2 + y^2}))$

To clean up the notation, let's invent two new functions. Let $f(x) = \ln(p(x))$ and let $h(r^2) = \ln(g(r))$. The equation simplifies to:

$f(x) + f(y) = h(x^2 + y^2)$

Step 3: The calculus trick Now we take the partial derivative of both sides with respect to $x$. Since $y$ is independent, the $f(y)$ term treats $x$ as a constant and becomes 0. Using the chain rule on the right side:

$f'(x) = 2x \cdot h'(x^2 + y^2)$

Divide both sides by $x$ to isolate the functions:

$\frac{f'(x)}{2x} = h'(x^2 + y^2)$

Here is the "aha!" moment. If we had taken the partial derivative with respect to $y$ instead, we would have found:

$\frac{f'(y)}{2y} = h'(x^2 + y^2)$

Since both expressions equal the same right-hand side, we can set them equal to each other:

$\frac{f'(x)}{x} = \frac{f'(y)}{y}$

Step 4: The constant of proportionality Look closely at that last equation. The left side depends entirely on $x$. The right side depends entirely on $y$. Since $x$ and $y$ are completely independent variables, the only mathematical way a function of purely $x$ can always equal a function of purely $y$ is if they are both equal to the exact same constant. Let's call that constant $c$.

$\frac{f'(x)}{x} = c$

Step 5: Solving the differential equation Multiply by $x$:

$f'(x) = cx$

Integrate both sides with respect to $x$:

$f(x) = \frac{1}{2}cx^2 + C$

Remembering our earlier substitution, $f(x) = \ln(p(x))$:

$\ln(p(x)) = \frac{1}{2}cx^2 + C$

Exponentiate both sides to solve for our final probability distribution $p(x)$:

$p(x) = e^{\frac{1}{2}cx^2 + C} = e^C \cdot e^{\frac{1}{2}cx^2}$

Are we clueless about these constants though? Think of more logical arguments, because they might help us constraint these currently arbitrary constant factors.

The Final Form

We can clean up those constants. Let's call $e^C$ a new constant $A$.

We also know that the total probability of the dart landing somewhere must be 1. For the area under the curve to be finite, the probability must decay toward 0 as $x$ gets infinitely far away from the bullseye. This forces our constant $c$ to be a negative number. We can replace $\frac{1}{2}c$ with a negative constant, $-k$.

This gives us the foundational equation for the normal curve:

$p(x) = A e^{-k x^2}$

To complete it, welcome one of the most famous and clever tricks in calculus: solving the Gaussian integral.

Step 1: Solving for $A$ (The Gaussian Integral)

We know the total probability must equal 1. So, the area under our foundational curve from $-\infty$ to $\infty$ must be 1:

$A \int_{-\infty}^{\infty} e^{-kx^2} dx = 1$

The integral $\int_{-\infty}^{\infty} e^{-kx^2} dx$ is notoriously impossible to solve using standard antiderivatives. The trick is to square the integral and solve it in 2D space. Let's call the integral $I$:

$I = \int_{-\infty}^{\infty} e^{-kx^2} dx$

If we multiply it by the exact same integral (just using $y$ as a dummy variable), we get $I^2$:

$I^2 = \left(\int_{-\infty}^{\infty} e^{-kx^2} dx\right) \left(\int_{-\infty}^{\infty} e^{-ky^2} dy\right)$

Combine them into a double integral:

$I^2 = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-k(x^2 + y^2)} dx dy$

Now we use the rotational symmetry we established at the very beginning. We convert this from Cartesian coordinates $(x, y)$ to polar coordinates $(r, \theta)$.

• $x^2 + y^2 = r^2$
• The area element $dx dy$ becomes $r dr d\theta$
• The limits change: an infinite plane means the radius $r$ goes from $0$ to $\infty$, and the angle $\theta$ goes from $0$ to $2\pi$.

$I^2 = \int_{0}^{2\pi} \int_{0}^{\infty} e^{-kr^2} r dr d\theta$

Now it's easily solvable using standard u-substitution. Let $u = -kr^2$, which means $du = -2kr dr$. Evaluating the inner integral for $r$:

$\int_{0}^{\infty} e^{-kr^2} r dr = \left[ -\frac{1}{2k} e^{-kr^2} \right]_{0}^{\infty} = (0) - \left(-\frac{1}{2k}\right) = \frac{1}{2k}$

Now plug that into the outer integral for $\theta$:

$I^2 = \int_{0}^{2\pi} \frac{1}{2k} d\theta = \frac{1}{2k} (2\pi) = \frac{\pi}{k}$

Since $I^2 = \frac{\pi}{k}$, we take the square root to find our original integral $I$:

$I = \sqrt{\frac{\pi}{k}}$

Plug this back into our very first equation:

$A \sqrt{\frac{\pi}{k}} = 1 \implies A = \sqrt{\frac{k}{\pi}}$

So our probability density function is now:

$p(x) = \sqrt{\frac{k}{\pi}} e^{-kx^2}$

---

Step 2: Swapping $k$ for Variance ($\sigma^2$)

In statistics, variance ($\sigma^2$) is the expected value of the squared deviation from the mean. Since our dartboard bullseye is at 0, the mean is 0. The formula for variance is:

$\sigma^2 = \int_{-\infty}^{\infty} x^2 p(x) dx$

Plug in our new function:

$\sigma^2 = \int_{-\infty}^{\infty} x^2 \left( \sqrt{\frac{k}{\pi}} e^{-kx^2} \right) dx$

You solve this using integration by parts. Let $u = x$ and $dv = x e^{-kx^2} dx$. After applying the integration by parts formula $\int u dv = uv - \int v du$ and calculating the boundaries, it simplifies beautifully to:

$\sigma^2 = \frac{1}{2k}$

Solving for $k$, we get:

$k = \frac{1}{2\sigma^2}$

---

Step 3: The Final Textbook Formula

Now we just plug $k = \frac{1}{2\sigma^2}$ into everything.

First, update our constant $A$:

$A = \sqrt{\frac{k}{\pi}} = \sqrt{\frac{1/(2\sigma^2)}{\pi}} = \sqrt{\frac{1}{2\pi\sigma^2}} = \frac{1}{\sigma\sqrt{2\pi}}$

Next, update the exponent:

$-kx^2 = -\left(\frac{1}{2\sigma^2}\right)x^2 = -\frac{x^2}{2\sigma^2}$

Put it all together, and the completed derivation yields the exact modern textbook definition of the normal distribution:

$p(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{x^2}{2\sigma^2}}$