The Normal Distribution: A Derivation through MultiVariable Calculus
Posted 2023/08/29
The normal distribution or Gaussian distribution (informally known as a bell curve) is the most used continuous probability distribution with many applications in statistics and science. Its importance comes from the central limit theorem, which establishes that the distribution of a sum (or sample mean) of many independent and identically distributed random variables will tend to look like a standard normal distribution with a large enough sample size. This property makes the distribution commonly used to represent random variables whose distributions are not known.
As such, the distribution is commonly taught as a tool in many courses without a full explanation of all of its properties and derivation. For me, the distribution was recently used in my SYDE162 Human Factors in Design course for calculating the distribution of anthropometric dimensions using zscores, but no explanation of what the distribution is and its properties were presented in the course. However, the distribution and many of its properties are easy to understand and the HershelMaxwell derivation can be understood with some basic understanding of multivariable differentiation and integration.
Basic Properties
Continuous probability distributions like the normal distribution have a cumulative distribution (CDF) and a probability density function (PDF). However, in this post I will be focusing on the probability density function which is the “bell curve” function most people use; the cumulative distribution function can be found by taking the antiderivative of the probability density function we will be deriving.
The probability density function describes how likely a random outcome is with outcomes near the mean being more likely and outcomes further away from the mean being less likely. This creates a function that looks like a “bell curve” as seen in the image on the right.
An example of this distribution being used is the sum of die rolls. In a two die roll, rolling a sum of 7 is more likely than the more extreme values of 2 and 12, graphing the probabilities of each outcome from 212 will give you a graph that vaguely looks like a normal distribution. However, if you increased the number of die and redrew the graph, the distribution will continue to look closer to the normal “bell curve” distribution as you increase the number of die. This tendency to approach the probability density function of the normal distribution as sample size is increased is called the central limit theorem, and the theorem can be seen in many applications which is why this distribution and its probability density function is so important.
Probability density functions (in general for any probability distribution) have a few properties that are important:

The probability density function is an integrable over its domain (or else our distribution would not have a cumulative distribution function and property number 4 would not work)

The probability density function is positive everywhere in its domain, $p(x)>0$ (the probability of an outcome is never negative)

The area under the curve of the probability density function (area between the curve and the $x$axis) is $1$ ($\int_{\infty}^\infty p(x) \mathit{dx} = 1$)

The probability that a random value $X$ lies in the interval $[a,b]$ is given by the area under the curve over that interval, $P(X \in [a,b]) = \int_a^b p(x) \mathit{dx}$
This last property is what probability density functions are most commonly used for, they help calculate the probability that a random variable will be a certain value or fall within a certain range. The normal distribution specifically has a couple more properties that are useful and will be used in the derivation:
 The probability density function is symmetric around the point $x=\mu$ which is the mean of the distribution (and also the median and mode).
 The PDF is unimodal (has a single peak), this means that the first derivative (instantaneous rate of change) is positive at $x<\mu$, negative at $x>\mu$, and zero at $x=\mu$
The HerschelMaxwell Theorem and the Cartesian Dartboard Analogy
Earlier, I explained that we will derive the normal probability density function by following the HerschelMaxwell derivation. This derivation uses the HerschelMaxwell’s probability theory, which states that if the probability distribution of a vector in $\mathbb R^n$ is unchanged by (or independent of) rotation and the individual vector components are independently randomly distributed, then the components of the vector are all normally distributed in an identical manner. Following this theorem, we will derive the probability distribution in a basic twodimensional case before simplifying to onedimension as it makes some of the later integration easier.
A good analogy of this theorem and a twodimensional normal distribution is throwing darts at the origin of a twodimensional Cartesian plane. The darts are aimed at the origin, but random errors in the throw produce different results. Under this analogy, we can assume that small errors are more likely than large errors, errors of equal distance from the origin are equally likely (errors are independent of orientation/rotation), and errors in the two axes are independent and normally distributed.
Finding the Basic Shape of the Probability Density Function
Under the assumptions of independent axes and independence from rotation, we can start with this equation, which equates probabilities in Cartesian coordinates to probabilities in polar coordinates:
$g(r)=p(x)p(y)$
Here, $p(x)$ represents the probability density function on the $x$axis, $p(y)$ represents the function on the $y$axis, and $g(r)$ represents the probability density with respect to distance from the origin (or radius). We can then partially differentiate the equation with respect to $\theta$ or rotation around the origin to get the following:$0=p(x)\frac{dp(y)}{d\theta}+p(y)\frac{dp(x)}{d\theta}$
Observe that the probability with respect to radius ($g(r)$) is independent of orientation and as such differentiates to $0$ and the use product rule on the right side. We can then substitute in $x=r\cos\theta$ and $y=r\sin\theta$ to convert Cartesian coordinates into polar coordinates on the right side and then separate the variables:$\begin{align*} 0 &= p(x)\frac{dp(r\sin\theta)}{d\theta} + p(y)\frac{dp(r\cos\theta)}{d\theta} \\ 0 &= p(x)p'(y)(r\cos\theta) + p(y)p'(x)(r\sin\theta) \\ p(y)\,p'(x)\,(r\sin\theta) &= p(x)\,p'(y)\,(r\cos\theta) \\ y\,p(y)\,p'(x) &= x\,p(x)\,p'(y) \\ \frac{p'(x)}{x\,p(x)} &= \frac{p'(y)}{y\,p(y)} \end{align*}$
For this differential equation to be true for any independent $x$ and $y$ value, the ratio of the differential equation must be constant. This unknown constant can be written as $C$:$\frac{p'(x)}{x\,p(x)} = \frac{p'(y)}{y\,p(y)} = C$
We can then use this equation to solve for the basic form of the normal distribution’s probability density function:$\begin{align*} \frac{p'(x)}{x\,p(x)} &= C \\ \frac{p'(x)}{p(x)} &= Cx \\ \int \frac{p'(x)}{p(x)} \mathit{dx} &= \int Cx \mathit{dx} \\ \ln(p(x)) &= \frac{Cx^2}{2}+c \\ p(x) &= e^{\frac{Cx^2}{2}+c} \\ p(x) &= Ae^{\frac{C}{2}x^2} \end{align*}$
However, there are two conditions that need to be met for this to be a valid probability density function. First, we are assuming large deviations from the mean (or large errors in the dartboard analogy) are less likely than small errors and the area under the curve must converge to $1$. This condition makes the unknown constant $C$ negative, as if it was positive, large deviations would be more likely and the area would diverge to $\infty$. The second condition is that all probability density functions must be positive, this makes the unknown coefficient $A$ positive (the exponential function it scales is always positive). This allows us to rewrite the basic form of the function into the following:
$p(x) = Ae^{\frac{k}{2}x^2}\qquad A,k>0$
Finding the Coefficient $A$
The value of $A$ in this basic function form can be found by using the area property of probability density functions (area under the curve is equal to $1$):
$\begin{align*} \int_{\infty}^\infty Ae^{\frac{k}{2}x^2} \mathit{dx} &= 1 \\ A\int_{\infty}^\infty e^{\frac{k}{2}x^2} \mathit{dx} &= 1 \\ \int_{\infty}^\infty e^{\frac{k}{2}x^2} \mathit{dx} &= \frac{1}{A} \end{align*}$
The function is symmetric and in its current form is an even function, so the integral can be rewritten like this:$\begin{align*} 2\int_0^\infty e^{\frac{k}{2}x^2} \mathit{dx} &= \frac{1}{A} \\ \int_0^\infty e^{\frac{k}{2}x^2} \mathit{dx} &= \frac{1}{2A} \end{align*}$
It is difficult to integrate a Gaussian function ($e^{x^2}$) so we will bring back the seconddimension from earlier (which looks identical due to symmetry) to make integration easier by rewriting the equation into a double integral:$\begin{align*} \left(\int_0^\infty e^{\frac{k}{2}x^2} \mathit{dx}\right)\left(\int_0^\infty e^{\frac{k}{2}y^2} \mathit{dy}\right) &= \left(\frac{1}{2A}\right)\left(\frac{1}{2A}\right) \\ \int_0^\infty \int_0^\infty e^{\frac{k}{2}x^2} \cdot e^{\frac{k}{2}y^2} \mathit{dy}\,\mathit{dx} &= \frac{1}{4A^2} \\ \int_0^\infty \int_0^\infty e^{\frac{k}{2}(x^2+y^2)} \mathit{dy}\,\mathit{dx} &= \frac{1}{4A^2} \\ \end{align*}$
This double integral is symmetrical around the origin, so we can convert it from Cartesian coordinates to polar coordinates. Note the appearance of $r$ as the Jacobian determinant from this change of coordinates. Following this change of coordinates, we can use standard integration techniques to solve for $A$:$\begin{align*} \int_0^\frac{\pi}{2} \int_0^\infty e^{\frac{k}{2}r^2}r \,\mathit{dr}\,\mathit{d\theta} &= \frac{1}{4A^2} \\ \int_0^\frac{\pi}{2}\mathit{d\theta} \cdot \int_0^\infty e^{\frac{k}{2}r^2}r \,\mathit{dr} &= \frac{1}{4A^2} \\ \Big[\theta\Big]^{\theta=\frac{\pi}{2}}_{\theta=0} \cdot \int_0^\infty e^{\frac{k}{2}r^2}r \,\mathit{dr} &= \frac{1}{4A^2} \\ \frac{\pi}{2} \cdot \int_0^\infty e^{\frac{k}{2}r^2}r \,\mathit{dr} &= \frac{1}{4A^2} \\ \int_0^\infty e^{\frac{k}{2}r^2}r \,\mathit{dr} &= \frac{1}{2\pi A^2} \end{align*}$
This single improper integral can be evaluated with an $u$substitution of $u=\frac{k}{2}r^2$ ($\mathit{du} = kr \mathit{dr}$), allowing the rest of the equation to be solved:$\begin{align*} \frac{1}{k}\int_0^{\infty} e^{u} \,\mathit{du} &= \frac{1}{2\pi A^2} \qquad\text{($u$substitution)} \\ \int_0^{\infty} e^{u} \,\mathit{du} &= \frac{k}{2\pi A^2} \\ \lim_{t\rightarrow\infty}\int_0^t e^{u} \,\mathit{du} &= \frac{k}{2\pi A^2} \qquad\text{(rewriting the improper integral as a limit)} \\ \lim_{t\rightarrow\infty}\Big[e^u\Big]^{u=t}_{u=0} &= \frac{k}{2\pi A^2} \\ \lim_{t\rightarrow\infty}\left(e^te^0\right) &= \frac{k}{2\pi A^2} \\ \lim_{t\rightarrow\infty}\left(\frac{1}{e^t}1\right) &= \frac{k}{2\pi A^2} \\ \left(01\right) &= \frac{k}{2\pi A^2} \\ 1 &= \frac{k}{2\pi A^2} \\ A^2 &= \frac{k}{2\pi} \\ A &= \sqrt{\frac{k}{2\pi}} \end{align*}$
Note that $A\neq\sqrt{\frac{k}{2\pi}}$ due to the condition that $A$ must be positive that we deduced above. This gives us a probability function of the form $p(x)=\sqrt{\frac{k}{2\pi}} e^{\frac{k}{2}x^2}$ with $k$ being the last step we need to solve for.Finding the Value $k$
To solve for the value $k$, we have to bring in the variance of the probability distribution. The variance is the square of the perfect deviation and is defined as $\sigma^2 = \int_{\infty}^{\infty}(x\mu)^2\,p(x)\,\mathit{dx}$. The mean of the function ($\mu$) is defined by the integral $\int_{\infty}^\infty x\,p(x)\,\mathit{dx}$. This integral is odd given the even probability function we found above, so our mean in this case is $0$. We can substitute this value and our probability density function into the integral for variance (squared deviation) to solve for $k$:
$\begin{align*} \int_{\infty}^{\infty}(x\mu)^2\,p(x)\,\mathit{dx} &= \sigma^2 \\ \int_{\infty}^{\infty}(x0)^2\,p(x)\,\mathit{dx} &= \sigma^2 \\ \int_{\infty}^{\infty}x^2\cdot\sqrt{\frac{k}{2\pi}} e^{\frac{k}{2}x^2}\,\mathit{dx} &= \sigma^2 \\ 2\sqrt{\frac{k}{2\pi}}\int_0^{\infty}x^2\cdot e^{\frac{k}{2}x^2}\,\mathit{dx} &= \sigma^2 \qquad\text{(simplify even integral)} \\ 2\sqrt{\frac{k}{2\pi}}\int_0^{\infty}x\cdot xe^{\frac{k}{2}x^2}\,\mathit{dx} &= \sigma^2 \end{align*}$
This integral can then be solved with an integration by parts with $u=x$ and $dv=xe^{\frac{k}{2}x^2}$:$\begin{align*} 2\sqrt{\frac{k}{2\pi}}\left(\lim_{t\rightarrow\infty}\left[x\cdot\frac{1}{k}\,e^{\frac{k}{2}x^2}\right]_{x=0}^{x=t}\int_0^{\infty}\frac{1}{k}e^{\frac{k}{2}x^2}\,\mathit{dx}\right) &= \sigma^2 \\ 2\sqrt{\frac{k}{2\pi}}\left(\frac{1}{k}\lim_{t\rightarrow\infty}\left[\frac{x}{e^{\frac{k}{2}x^2}}\right]_{x=0}^{x=t}+\frac{1}{k}\int_0^{\infty}e^{\frac{k}{2}x^2}\,\mathit{dx}\right) &= \sigma^2 \\ 2\sqrt{\frac{k}{2\pi}}\left(\frac{1}{k}\lim_{t\rightarrow\infty}\left(\frac{t}{e^{\frac{k}{2}t^2}}\frac{0}{e^{\frac{k}{2}0^2}}\right) + \frac{1}{k}\int_0^{\infty}e^{\frac{k}{2}x^2}\,\mathit{dx}\right) &= \sigma^2 \\ 2\sqrt{\frac{k}{2\pi}}\left(\frac{1}{k}\lim_{t\rightarrow\infty}\left(\frac{1}{e^{\frac{k}{2}t^2}\cdot{k}t}\right) + \frac{1}{k}\int_0^{\infty}e^{\frac{k}{2}x^2}\,\mathit{dx}\right) &= \sigma^2 \qquad\text{(L'Hôpital's rule)} \\ 2\sqrt{\frac{k}{2\pi}}\cdot\frac{1}{k}\int_0^{\infty}e^{\frac{k}{2}x^2}\,\mathit{dx} &= \sigma^2 \end{align*}$
This final integral we have already solved above when finding the value of $A$:$\begin{align*} 2\sqrt{\frac{k}{2\pi}}\cdot\frac{1}{2k}\int_{\infty}^{\infty}e^{\frac{k}{2}x^2}\,\mathit{dx} &= \sigma^2 \qquad(\mathbb R^+\rightarrow\mathbb R) \\ 2\sqrt{\frac{k}{2\pi}}\cdot\frac{1}{2k}\cdot\frac{1}{A} &= \sigma^2 \\ 2\frac{\sqrt{k}}{\sqrt{2\pi}}\cdot\frac{1}{2k}\cdot\frac{\sqrt{2\pi}}{\sqrt{k}} &= \sigma^2 \\ \frac{1}{k} &= \sigma^2 \\ k &= \frac{1}{\sigma^2} \end{align*}$
The Normal Probability Density Function
We can plug in this value of $k$ to find the normal probability distribution function:
$\begin{align*} p(x)=\sqrt{\frac{k}{2\pi}}\,e^{\frac{k}{2}x^2} \\ p(x)=\sqrt{\frac{1}{2\pi\sigma^2}}\,e^{\frac{1}{2\sigma^2}x^2} \\ p(x)=\frac{1}{\sigma\sqrt{2\pi}}\,e^{\frac{1}{2}\left(\frac{x}{\sigma}\right)^2} \\ \end{align*}$
This form of the normal probability density function has a single constant $\sigma$, the desired standard deviation of the distribution, and has a mean ($\mu$) of $0$. However, most applications of this distribution don’t have a mean of $0$, so we can horizontally translate the function so that the peak of the “bell curve” lies on the desired mean:$p(x)=\frac{1}{\sigma\sqrt{2\pi}}\,e^{\frac{1}{2}\left(\frac{x\mu}{\sigma}\right)^2}$
Personal Thoughts
I first worked on a presentation about this derivation of the normal distribution with a group of friends back in my Grade 12 Data Management course in high school as a final project. That presentation and this post were loosely based on a “dartboard” based derivation written by Dan Teague which I did not fully understand back in high school due to some of the derivation being beyond my limited high school level calculus knowledge. After finishing my Calculus 1 and 2 courses (SYDE111 & SYDE112) and seeing more applications of the distribution during my studies at the University of Waterloo, I was inspired to revisit the derivation and hopefully write it out in a more comprehensive form than the original derivation.
I am not a mathematician or educator, so there may be errors post above or potential improvements I can make for clarity. If you spot any errors or improvements that can be made, feel free to send me an email!
3Blue1Brown also recently covered the Central Limit Theorem in a series of videos that I found very interesting available here which demonstrate why the normal probability distribution is so useful.
I haven’t posted anything to my blog in a while so hopefully this marks the beginning of more technical writing posts related to math, programming, and information security!