The Gaussian Integral and the Gaussian Probability Density Function

Abstract.

We review moment-generating functions to provide a context for understanding the utility of Gaussian functions and the Gaussian integral. We then explain Gaussian functions as a class of exponential function and demonstrate a common technique for integrating a Gaussian function over $(-\mathrm{\infty },\mathrm{\infty })$ by converting to polar coordinates, Finally, we use a moment-generating function and the Gaussian integral to construct the parameterized form of the Gaussian or normal probability density function.

1. Characterizing Probability Distributions

Although we could start by presenting a Gaussian function and proceed by evaluating its integral over the real numbers, that would not provide a context for the exercise. The relevance of the procedure would be lost and the reader would be left with an isolated mathematical process devoid of applicability. Given that the root of Gaussian functions lies in probability theory, where a specific instance defines the so-called normal distribution, we will review the necessary statistical principles to understand the utility of the Gaussian integral.

1.1. Moments

In probability theory, a parameter associated with a probability distribution, such as its mean or variance, is said to be a characteristic of the distribution. A probability distribution is a function that models a so-called population of data points. A population is the hypothetical collection of all possible measurements made under a given set of conditions. It is a theoretical set of data and not the actual set of measurements. To distinguish the characteristics of a sample from the parameters of a distribution, we use the term statistic. Therefore, the population mean is a parameter of the probability distribution and the sample mean is a statistic of the sample set.

To illustrate the difference between a parameter and a statistic, let’s compare the formula for the expected value of a discrete random variable $X$ with probability mass function $P(x)$,

$$E(X)={\mu }_X=\sum _{n=1}^N{x}_{n}P(x_n),$$ (1)

to the formula for the expected value of a continuous random variable $X$ with probability density function $P(x)$,

$$E(X)={\mu }_X={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}xP(x)dx.$$ (2)

Equation 1 represents a statistic and may have different values for different sample sets. Equation 2 represents a parameter and has a fixed value. Each is a characteristic of its respective discrete and continuous distribution.

A moment is an element of a set of characteristic values that collectively describe the distribution of a random variable, albeit not in a unique manner. Under the right conditions, two random variables with identical moments will have identical probability distributions. That makes it possible to approximate the probability distribution of a random variable using moments.

In the study of statistics, moments fall into two classes: raw and central. A raw moment is computed about the origin and a central moment is computed about the mean. Without proof, we state that if $X$ is a discrete random variable with probability mass function $P(x)$ and $g(X)$ is a function of $X$, then

$$E(g(X))={\mu }_{g(X)}=\sum _{n=1}^{N}g(x_n)\cdot P(x_n);$$ (3)

and if $X$ is a continuous random variable with probability density function $P(x)$ and $g(X)$ is a function of $X$, then

$$E(g(X))={\mu }_{g(X)}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}g(x)\cdot P(x)dx.$$ (4)

We set $g(X)=X^k$ to find the raw moments of a probability distribution.

The raw moments of a discrete probability distribution are defined as

$$E(X^k)={\mu }_k^{\prime }=\sum _{n=1}^N{x}_n^{k}P(x_n),$$ (5)

where $E(X^k)$ is the $k^{\mathrm{th}}$ moment of discrete random variable $X$ and written as ${\mu }_k^{\prime }$. Similarly, the raw moments of a continuous probability distribution for a continuous random variable $X$ are defined as

$$E(X^k)={\mu }_k^{\prime }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{x}^{k}P(x)dx.$$ (6)

The central moments are defined analogously, setting $g(X)={(X-\mu )}^k$ to center the moments about the mean $\mu$. The central moments of a discrete random variable are

$E({[X-\mu ]}^k)={\mu }_k$ $={\displaystyle \sum _{n=1}^N}{(x_n-\mu )}^{k}P(x_n)$ (7)

and for a continuous random variable they are

$E({[X-\mu ]}^k)={\mu }_k$ $={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{(x-\mu )}^{k}P(x)dx.$ (8)

Note that the first raw moment, ${\mu }_1^{\prime }$, is equal to the mean, $\mu$. The first central moment, ${\mu }_1$, will always equal zero because Equation 8 reduces to $\mu -\mu =0$ for $k=1$. The area under the curve to the left of the mean will always be equal to that to the right for a proper probability density function. The second central moment, ${\mu }_2$, is equal to the population variance. Knowing that the population variance, ${\sigma }^2$, of random variable $X$ is $E(X^2)-E^2(X)$, the variance can be expressed in terms of raw moments as

$${\sigma }^2={\mu }_2^{\prime }-{\mu }_1^{\prime 2}.$$ (9)

In addition, it is possible to define a set of moments about an arbitrary point, $a$, in which case you would replace $\mu$ in Equations 7 and 8 with the point $a$.

1.2. Moment-Generating Functions

Instead of calculating the moments of a probability distribution on an ad hoc basis, we can identify a function from which the moments can be immediately derived. This so-called moment-generating function exploits the series expansion of ${\mathrm{e}}^{tx}$ to create a function that, when successively differentiated at zero, produces the raw moment of a probability distribution corresponding to the level of differentiation. For example, the first derivative at zero produces the first raw moment, the second derivative the second raw moment, and so on. This works by applying Equation 4 to find the expected value of ${\mathrm{e}}^{tX}$, using $g(X)={\mathrm{e}}^{tX}$ and integrating to yield a moment-generating function,

$$M(t)=E({\mathrm{e}}^{tX})={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{tx}P(x)dx.$$ (10)

Expanding ${\mathrm{e}}^{tx}$ to its infinite series gives us

$M(t)$ $={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}\left({\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{\displaystyle \frac{t^n{x}^n}{n!}}\right)P(x)dx;$

$={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}(1+tx+{\displaystyle \frac{t^2{x}^2}{2!}}+{\displaystyle \frac{t^3{x}^3}{3!}}+\mathrm{\cdots })P(x)dx;$

$={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}P(x)dx+t{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}xP(x)dx+{\displaystyle \frac{t^2}{2!}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{x}^{2}P(x)dx+\mathrm{\cdots };$

$=1+t{\mu }_1^{\prime }+{\displaystyle \frac{t^2}{2!}}{\mu }_2^{\prime }+{\displaystyle \frac{t^3}{3!}}{\mu }_3^{\prime }+\mathrm{\cdots }.$ (11)

We arrive at Equation 11 by substituting terms using Equation 6.

If we differentiate $M(t)$ with respect to $t$ we get

$$M^{\prime }(t)=\frac{dM}{dt}={\mu }_1^{\prime }+\frac{2t}{2!}{\mu }_2^{\prime }+\frac{3t^2}{3!}{\mu }_3^{\prime }+\mathrm{\cdots }.$$

Evaluating $M^{\prime }(t)$ at $t=0$ eliminates all but the first term, yielding

$$M^{\prime }(0)={\frac{dM}{dt}|}_{t=0}={\mu }_1^{\prime }.$$

Taking the second derivative produces

$$M^{\prime \prime }(t)=\frac{d^{2}M}{d{t}^2}=\frac{2}{2!}{\mu }_2^{\prime }+\frac{6t}{3!}{\mu }_3^{\prime }+\frac{12t^2}{4!}+\mathrm{\cdots };$$

and evaluating it at $t=0$ results in

$$M^{\prime \prime }(0)={\frac{d^{2}M}{d{t}^2}|}_{t=0}={\mu }_2^{\prime }.$$

It should be apparent that the $k^{\mathrm{th}}$ derivative evaluated at $t=0$ produces the $k^{\mathrm{th}}$ raw moment according to

$${\frac{d^k}{d{t}^k}M(t)|}_{t=0}={\mu }_k^{\prime }\mathit{\hspace{1em}}\text{where}\mathit{\hspace{1em}}k>0.$$ (12)

As with Equation 4, if $X$ is a continuous random variable with probability density function $P(x)$ and $g(X)$ is a function of $X$, then

$$M(t)=E({\mathrm{e}}^{tg(X)})={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{tg(x)}P(x)dx.$$ (13)

We can use this formulation to find the moment-generating function of $X-\mu$ by setting $g(X)=X-\mu$. This can be useful to derive the variance directly from the second moment instead of using Equation 9. However, we will only be using raw moments.

1.3. Uniform Distribution Moments

So that we can understand how to derive the moment-generating function of the Gaussian distribution in Section 4, we will now work through the example of deriving the moment-generating function of the uniform probability distribution,

$$P(x)=\{\begin{array}{cc}\frac{1}{b-a}\hfill & \text{for}\hspace{1em}a\le x\le b,\hfill \\ 0\hfill & \text{for}\hspace{1em}x<a\hspace{1em}\text{or}\hspace{1em}x>b.\hfill \end{array}$$ (14)

Applying Equation 10 to the unifrom probability density function from Equation 14, we have

$$M(t)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{tx}\frac{1}{b-a}dx=\frac{1}{b-a}{\int }_a^b{\mathrm{e}}^{tx}dx.$$

Integrating gives us

$$M(t)={\frac{1}{b-a}\cdot \frac{{\mathrm{e}}^{tx}}{t}|}_a^b=\frac{1}{b-a}\left(\frac{{\mathrm{e}}^{bt}-{\mathrm{e}}^{at}}{t}\right),$$

which produces the indeterminate form $\frac{0}{0}$ at $t=0$, requiring us to apply l’Hôpital’s rule as follows:

$$\underset{t\to 0}{lim}\frac{1}{b-a}\left(\frac{{\mathrm{e}}^{bt}-{\mathrm{e}}^{at}}{t}\right)=\underset{t\to 0}{lim}\frac{b{\mathrm{e}}^{bt}-a{\mathrm{e}}^{at}}{b-a}=\frac{b-a}{b-a}=1.$$

This gives us the complete moment-generating function for the uniform probability distribution,

$$M(t)=\{\begin{array}{cc}\frac{{\mathrm{e}}^{bt}-{\mathrm{e}}^{at}}{(b-a)t}\hfill & \text{for}\hspace{1em}t\ne 0,\hfill \\ 1\hfill & \text{for}\hspace{1em}t=0.\hfill \end{array}$$ (15)

We can now find the first raw moment of the uniform probability distribution by finding the first derivative with respect to $t$ of its moment-generating function. Using the sum and product rules, we get

$M^{\prime }(t)={\displaystyle \frac{d}{dt}}M(t)$ $={\displaystyle \frac{1}{b-a}}\left({\displaystyle \frac{d}{dt}}{\displaystyle \frac{{\mathrm{e}}^{bt}}{t}}-{\displaystyle \frac{d}{dt}}{\displaystyle \frac{{\mathrm{e}}^{at}}{t}}\right);$

$={\displaystyle \frac{1}{b-a}}\left({\displaystyle \frac{b{\mathrm{e}}^{bt}}{t}}-{\displaystyle \frac{{\mathrm{e}}^{bt}}{t^2}}-{\displaystyle \frac{a{\mathrm{e}}^{at}}{t}}+{\displaystyle \frac{{\mathrm{e}}^{at}}{t^2}}\right);$

$={\displaystyle \frac{1}{b-a}}\left({\displaystyle \frac{b{\mathrm{e}}^{bt}-a{\mathrm{e}}^{at}}{t}}+{\displaystyle \frac{{\mathrm{e}}^{at}-{\mathrm{e}}^{bt}}{t^2}}\right).$

The resulting function is indeterminate at $t=0$, requiring us to rewrite it to produce a result of the form $\frac{0}{0}$, allowing us to apply l’Hôpital’s rule as follows:

${{\displaystyle \frac{d}{dt}}M(t)|}_{t=0}$ $=\underset{t\to 0}{lim}{\displaystyle \frac{1}{b-a}}\left({\displaystyle \frac{bt{\mathrm{e}}^{bt}-at{\mathrm{e}}^{at}+{\mathrm{e}}^{at}-{\mathrm{e}}^{bt}}{t^2}}\right);$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}{\displaystyle \frac{b^{2}t{\mathrm{e}}^{bt}+b{\mathrm{e}}^{bt}-a^{2}t{\mathrm{e}}^{at}-a{\mathrm{e}}^{at}+a{\mathrm{e}}^{at}-b{\mathrm{e}}^{bt}}{2t}};$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}{\displaystyle \frac{b^{2}t{\mathrm{e}}^{bt}-a^{2}t{\mathrm{e}}^{at}}{2t}};$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}{\displaystyle \frac{b^2{\mathrm{e}}^{bt}-a^2{\mathrm{e}}^{at}}{2}};$

$={\displaystyle \frac{1}{b-a}}\cdot {\displaystyle \frac{b^2-a^2}{2}}={\displaystyle \frac{b+a}{2}}.$

Swapping the terms from the last result, the first raw moment of the uniform probability distribution is

$$M^{\prime }(0)={\mu }_1^{\prime }=\frac{a+b}{2}.$$ (16)

We can verify the first raw moment is correct by using Equation 2 to find the mean,

$\mu$ $={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}x\cdot {\displaystyle \frac{1}{b-a}}dx;$

$={\displaystyle \frac{1}{b-a}}{\displaystyle {\int }_a^b}xdx;$

$={{\displaystyle \frac{1}{b-a}}\cdot {\displaystyle \frac{x^2}{2}}|}_a^b;$

$={\displaystyle \frac{1}{b-a}}\cdot {\displaystyle \frac{b^2-a^2}{2}}={\displaystyle \frac{a+b}{2}}.$

As a final exercise, let’s derive the second raw moment of the uniform probability distribution. Differentiating the first derivative from before, we find the second derivative is

$$M^{\prime \prime }(t)=\frac{1}{b-a}\left(\frac{b^2{\mathrm{e}}^{bt}-a^2{\mathrm{e}}^{at}}{t}-\frac{b{\mathrm{e}}^{bt}-a{\mathrm{e}}^{at}}{t^2}+\frac{a{\mathrm{e}}^{at}-b{\mathrm{e}}^{bt}}{t^2}-2\left(\frac{{\mathrm{e}}^{at}-{\mathrm{e}}^{bt}}{t^3}\right)\right).$$

Again, we have a function that cannot be readily evaluated at $t=0$. As before, we restructure it to allow the application of l’Hôpital’s rule, giving us

${{\displaystyle \frac{d^2}{d{t}^2}}M(t)|}_{t=0}$ $={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}({\displaystyle \frac{t^2(b^2{\mathrm{e}}^{bt}-a^2{\mathrm{e}}^{at})}{t^3}}-{\displaystyle \frac{t(b{\mathrm{e}}^{bt}-a{\mathrm{e}}^{at})}{t^3}}$

$+{\displaystyle \frac{t(a{\mathrm{e}}^{at}-b{\mathrm{e}}^{bt})}{t^3}}-{\displaystyle \frac{2({\mathrm{e}}^{at}-{\mathrm{e}}^{bt})}{t^3}});$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}{\displaystyle \frac{t^2(b^2{\mathrm{e}}^{bt}-a^2{\mathrm{e}}^{at})+2t(a{\mathrm{e}}^{at}-b{\mathrm{e}}^{bt})+2({\mathrm{e}}^{bt}-{\mathrm{e}}^{at})}{t^3}};$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}({\displaystyle \frac{t^2(b^3{\mathrm{e}}^{bt}-a^3{\mathrm{e}}^{at})}{3t^2}}$

$+{\displaystyle \frac{2t(b^2{\mathrm{e}}^{bt}-a^2{\mathrm{e}}^{at})+2t(a^2{\mathrm{e}}^{at}-b^2{\mathrm{e}}^{bt})}{3t^2}}$

$+{\displaystyle \frac{2(a{\mathrm{e}}^{at}-b{\mathrm{e}}^{bt})+2(b{\mathrm{e}}^{bt}-a{\mathrm{e}}^{at})}{3t^2}});$

$={\displaystyle \frac{1}{b-a}}\underset{t\to 0}{lim}{\displaystyle \frac{(b^3{\mathrm{e}}^{bt}-a^3{\mathrm{e}}^{at})}{3}};$

$={\displaystyle \frac{b^3-a^3}{3(b-a)}}={\displaystyle \frac{a^2+ab+b^2}{3}}.$

Using Equation 9, we can formulate the variance of the uniform probability distribution as

${\sigma }^2$ $={\displaystyle \frac{a^2+ab+b^2}{3}}-{\left({\displaystyle \frac{a+b}{2}}\right)}^2;$

$={\displaystyle \frac{4a^2+4ab+4b^2-3a^2-6ab-3b^2}{12}};$

$={\displaystyle \frac{a^2-2ab+b^2}{12}}={\displaystyle \frac{b^2-2ab+a^2}{12}};$

$={\displaystyle \frac{{(b-a)}^2}{12}}.$ (17)

2. Gaussian Functions

Mathematics literature uses the term Gaussian function either narrowly—using it exclusively to refer to the Gaussian probability density function—or broadly—using it to refer to a class of exponential functions. We will use the term in its broad sense, starting with the simplest instance, $f(x)={\mathrm{e}}^{-x^2}$, before examining more complex forms. Unfortunately, the author has failed to ascertain the historical origin of Gaussian functions and how their earliest applications eventually led to the identification of the Gaussian probability density function. Therefore, we will present a pedagogical analysis of how starting from $f(x)={\mathrm{e}}^{-x^2}$ eventually leads to the Gaussian probability density function.

Generally speaking, a Gaussian function is a function whose natural logarithm is a concave quadratic function. A quadratic function is concave if its second derivative is negative. In essence, it is a downward-growing parabola as in Figure 1.

Figure 1. A concave quadratic function.

Therefore, we can define a Gaussian function as having the form $f(x)=a{\mathrm{e}}^{-b{x}^2+cx+d}$, where $b$ is positive in order to ensure the second derivative of the quadratic function is negative. To simplify our discussion, we restrict our scrutiny to quadratic functions of the form $-a{(x-b)}^2$, giving $f(x)=a{\mathrm{e}}^{-b{(x-c)}^2}$ as our model Gaussian function. The simplest Gaussian function we can construct sets $a=1$, $b=1$, and $c=0$, leaving us with $f(x)={\mathrm{e}}^{-x^2}$. We plot this function in Figure 2 and show how it varies by changing the function parameters.

Figure 2. The Gaussian function, $f(x)={\mathrm{e}}^{-x^2}$, can be stretched, widened or narrowed, and shifted by changing $a$, $b$, and $c$ in $f(x)=a{\mathrm{e}}^{-b{(x-c)}^2}$.

Starting with $f(x)={\mathrm{e}}^{-x^2}$, we can visually note it is an even function, meaning that $f(x)=f(-x)$. This property means the function is symmetric about the $y$-axis. More generally, all Gaussian functions are symmetric about their midpoints, with the area under the curve to the left of the midpoint being equal to the area under the curve to the right of the midpoint. We will exploit this property in Section 3 to compute half of the Gaussian integral.

The midpoint of the function can be shifted by changing the value of $c$, as done in the plot of $f(x)={\mathrm{e}}^{-{(x-2)}^2}$ in Figure 2. Notice that the midpoint is shifted by an amount equal to $c$, shifting in the positive $x$ direction when $c$ is positive and in the negative $x$ direction when $c$ is negative. Changing the value of $b$ expands or contracts the width of the curve. A value greater than one narrows the curve and a value less than one and greater than zero widens the curve, as in the plot of $f(x)={\mathrm{e}}^{-x^2/4}$. Notice how the base of the curve expands by an amount equal to $\frac{1}{b}$. The function can be stretched upward or downward by multiplying it by a constant, as in the plot of $f(x)=2{\mathrm{e}}^{-x^2}$ where we see that the maximum of the function changes by a factor of $a$. The same effect can be achieved by adding a constant $k$ to the quadratic term, as in the plot of $f(x)={\mathrm{e}}^{-x^2+1}$, where $k=1$. But that is equivalent to multiplication by ${\mathrm{e}}^k$, which is one reason we chose not to use a general quadratic form for the exponent (another is that it simplifies integration by substitution as we will see in Section 4). It should be clear that $\mathrm{e}\cdot {\mathrm{e}}^{-x^2}={\mathrm{e}}^{-x^2+1}$ and $2{\mathrm{e}}^{-x^2}={\mathrm{e}}^{\ln 2}\cdot {\mathrm{e}}^{-x^2}={\mathrm{e}}^{-x^2+\ln 2}$, making an additive term redundant when parameterizing the function and studying the effects of changing parameter values. In Section 4, we will see that parameters $a$, $b$, and $c$ posses special meanings with respect to the Gaussian probability distribution.

3. The Gaussian Integral

Recognizing that $f(x)={\mathrm{e}}^{-x^2}$ ranges in value from $0$ to $1$, suppose we wanted to use a Gaussian function as a probability density function. Our first step would be to ensure that the integral over $(-\mathrm{\infty },\mathrm{\infty })$ is equal to $1$. We would have to integrate the function to determine the area, $A$, under the curve. If $A\ne 1$, we would normalize the function to $\frac{1}{A}f(x)$ so that its integral over $(-\mathrm{\infty },\mathrm{\infty })$ equaled $1$. We will, in fact, do all of this while deriving the Gaussian probability density function. However, embarking on this journey requires evaluating the Gaussian integral, for which we must first take a brief detour.

Solving problems in mathematics often requires what could be called tricks; like when we restructured the derivatives in Section 1.3 so we could apply l’Hôpital’s rule. The Gaussian integral can be solved in various ways, all of which require some trickery. We will be using two tricks, the first of which is based on the seemingly trivial identity $x=\sqrt{x^2}$, allowing us to observe that

${\displaystyle {\int }_a^b}{f}^{\prime }(x)dx$ $={\left({\displaystyle {\int }_a^b}{f}^{\prime }(x)dx\cdot {\displaystyle {\int }_a^b}{f}^{\prime }(x)dx\right)}^{\frac{1}{2}}$ (18)

implies

$f(b)-f(a)$ $=\sqrt{{(f(b)-f(a))}^2}.$ (19)

The second trick we will use is to convert from Cartesian to polar coordinates, for which the reader will have to refer to a calculus textbook if the process is unclear.

We wish to evaluate the improper integral

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx$$

to see if the integrand can be adapted for use as a probability density function. We cannot evaluate the integral based on an existing indefinite integral because the integrand has no elementary indefinite integral. To evaluate the improper integral will depend exactly on its improper nature, integrating from negative to positive infinity. First, we apply Equation 18 to rewrite the integral as

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx={\left({\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx\cdot {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx\right)}^{\frac{1}{2}}.$$

Next, we replace one of the $x$ dummy variables with another dummy variable, $y$, allowing us to rewrite the integral as the following double integral:

${\left({\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx\cdot {\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx\right)}^{\frac{1}{2}}$ $={\left({\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx\cdot {\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-y^2}dy\right)}^{\frac{1}{2}};$

$={\left({\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(x^2+y^2)}dxdy\right)}^{\frac{1}{2}}.$

For clarity, we’ll set aside the square root until the end and evaluate the double integral by converting to polar coordinates as follows:

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(x^2+y^2)}dxdy$ $={\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-r^2}drd\theta ;$

$={\displaystyle {\int }_0^{2\pi }}-{{\displaystyle \frac{1}{2}}{\mathrm{e}}^{-r^2}|}_0^{\mathrm{\infty }}d\theta ;$

$={\displaystyle {\int }_0^{2\pi }}{\displaystyle \frac{1}{2}}d\theta ;$

$={{\displaystyle \frac{1}{2}}\theta |}_0^{2\pi };$

$=\pi .$

Now we can apply the square root to arrive at our final result,

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx=\sqrt{\pi }.$$ (20)

The various transformations we performed were possible because $f(x)={\mathrm{e}}^{-x^2}$ is a continuous function. As noted in Section 2, $f(x)=f(-x)$ holds true for a Gaussian function. That allows us to infer that

$${\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx=\frac{\sqrt{\pi }}{2}.$$ (21)

More formally, if $f(x)$ is continuous on the interval $(-\mathrm{\infty },\mathrm{\infty })$ containing point $a$ then

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)dx={\int }_{-\mathrm{\infty }}^{a}f(x)dx+{\int }_a^{\mathrm{\infty }}f(x)dx.$$

This property of improper integrals allows us to write

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx={\int }_{-\mathrm{\infty }}^0{\mathrm{e}}^{-x^2}dx+{\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-x^2}dx.$$

Observing that $f(x)=f(-x)$ gives us

${\displaystyle {\int }_{-\mathrm{\infty }}^0}{\mathrm{e}}^{-x^2}dx$ $={\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx;$

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx$ $=2{\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx;$

${\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx$ $={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx;$

${\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx$ $={\displaystyle \frac{\sqrt{\pi }}{2}}.$

4. The Gaussian Probability Density Function

Having evaluated the most basic Gaussian integral as $\sqrt{\pi }$, we can start to build a probability density function derived from the simplest Gaussian function, $f(x)={\mathrm{e}}^{-x^2}$, normalizing it to $f(x)=\frac{1}{\sqrt{\pi }}{\mathrm{e}}^{-x^2}$ so that its integral over $(-\mathrm{\infty },\mathrm{\infty })$ evaluates to $1$. To investigate whether this is a useful construction, we can derive its moment-generating function and first two raw moments. Using Equation 10, we have

$M(t)$ $={\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{tx}{\mathrm{e}}^{-x^2}dx;$

$={\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(x^2-tx)}dx;$

$={\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(x^2-tx)}dx\cdot {\displaystyle \frac{{\mathrm{e}}^{-t^2/4}}{{\mathrm{e}}^{-t^2/4}}};$

$={\displaystyle \frac{{\mathrm{e}}^{t^2/4}}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(x^2-tx+t^2/4)}dx;$

$={\displaystyle \frac{{\mathrm{e}}^{t^2/4}}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-{(x-t/2)}^2}dx.$

Using $u=x-\frac{t}{2}$ and $du=1\cdot dx$, we find

$M(t)$ $={\displaystyle \frac{{\mathrm{e}}^{t^2/4}}{\sqrt{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-u^2}du;$

$={\displaystyle \frac{{\mathrm{e}}^{t^2/4}}{\sqrt{\pi }}}\cdot \sqrt{\pi };$

$={\mathrm{e}}^{t^2/4}.$

Note that we used Equation 20 to evaluate ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{e}}^{-{(x-t/2)}^2}dx$ via substitution.

We now use Equation 12 to find the first two raw moments. The first raw moment is

$${\frac{d}{dt}M(t)|}_{t=0}={\mu }_1^{\prime }=\underset{t\to 0}{lim}\frac{d}{dt}{\mathrm{e}}^{t^2/4}=\underset{t\to 0}{lim}\frac{t}{2}{\mathrm{e}}^{t^2/4}=0.$$

A mean of $0$ is exactly what we would expect for a probability density function centered about the origin. The second raw moment is

$${\frac{d^2}{d{t}^2}M(t)|}_{t=0}={\mu }_2^{\prime }=\underset{t\to 0}{lim}\frac{d^2}{d{t}^2}{\mathrm{e}}^{t^2/4}=\underset{t\to 0}{lim}\left(\frac{1}{2}{\mathrm{e}}^{t^2/4}+\frac{t^2}{4}{\mathrm{e}}^{t^2/4}\right)=\frac{1}{2},$$

giving us a variance of

$${\sigma }^2={\mu }_2^{\prime }-{\mu }_1^{\prime 2}=\frac{1}{2}-0=\frac{1}{2}$$

after using Equation 9.

The values of the raw moments we found are constants that don’t really help us toward formulating a general Gaussian probability density function. But they do give us clues about how such a function might be formulated. In Section 2, we saw how the center of a Gaussian function could be shifted by changing the value of the parameter $c$ in $f(x)=a{\mathrm{e}}^{-b{(x-c)}^2}$. We now have a clue that $c$ is the mean of the distribution because it lies at the midpoint about which the function has symmetric areas. The parameter $b$ may be related in some way to the variance, given how it affects the width of the curve. The parameter $a$ must be a normalizing constant. To test our theory, we can hold these parameters constant, with $b$ being positive, and try to find a moment-generating function parameterized in terms of $a$, $b$, and $c$.

We already know that $c$ does not affect the area under the curve. Therefore we will ignore it and evaluate the integral of $f(x)=a{\mathrm{e}}^{-b{x}^2}$ over $(-\mathrm{\infty },\mathrm{\infty })$ to find a value of $a$ dependent on $b$ that normalizes the function. That will enable us to find a moment-generating function for the function including $c$ using substitution as we did earlier. Using our earlier work, we can readily evaluate the integral as follows:

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-b{x}^2}dx$ $=a{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle \frac{{\mathrm{e}}^{-u^2}}{\sqrt{b}}}du$

where $u=\sqrt{b}x$ and $du=\sqrt{b}dx$. Substituting Equation 20, we get

${\displaystyle \frac{a}{\sqrt{b}}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-u^2}du$ $={\displaystyle \frac{a}{\sqrt{b}}}\cdot \sqrt{\pi }=a\sqrt{{\displaystyle \frac{\pi }{b}}}.$

Solving for $a$ in

$a\sqrt{{\displaystyle \frac{\pi }{b}}}$ $=1$

to normalize the function leaves us with

$a$ $=\sqrt{{\displaystyle \frac{b}{\pi }}}.$ (22)

Now we can use $f(x)=\sqrt{\frac{b}{\pi }}{\mathrm{e}}^{-b{(x-c)}^2}$ as a probability density function and find its moment-generating function. Starting with

$M(t)$ $={\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{tx}\cdot \sqrt{{\displaystyle \frac{b}{\pi }}}{\mathrm{e}}^{-b{(x-c)}^2}dx;$

$=\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{tx}{\mathrm{e}}^{-b{(x-c)}^2}dx;$

$=\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-(2bc+t)x+b{c}^2)}dx;$

$=\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-(2bc+t)x+b{c}^2)}dx\cdot {\displaystyle \frac{{\mathrm{e}}^{-(ct+t^2/4b)}}{{\mathrm{e}}^{-(ct+t^2/4b)}}};$

$={\mathrm{e}}^{ct+t^2/4b}\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-(2bc+t)x+b{c}^2+ct+t^2/4b)}dx;$

$={\mathrm{e}}^{ct+t^2/4b}\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-\left(b{x}^2-(2bc+t)x+{\left(\frac{bc+t/2}{\sqrt{b}}\right)}^2\right)}dx;$

$={\mathrm{e}}^{ct+t^2/4b}\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-{\left(\sqrt{b}x-\left(\frac{bc+t/2}{\sqrt{b}}\right)\right)}^2}dx;$

we can substitute $u=\sqrt{b}x-\left(\frac{bc+t/2}{\sqrt{b}}\right)$ and $du=\sqrt{b}dx$ to get

$M(t)$ $={\mathrm{e}}^{ct+t^2/4b}\sqrt{{\displaystyle \frac{b}{\pi }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle \frac{{\mathrm{e}}^{-u^2}}{\sqrt{b}}}du;$

$={\mathrm{e}}^{ct+t^2/4b}\sqrt{{\displaystyle \frac{b}{\pi }}}\cdot \sqrt{{\displaystyle \frac{\pi }{b}}}={\mathrm{e}}^{ct+t^2/4b}.$

The first moment is

${{\displaystyle \frac{d}{dt}}M(t)|}_{t=0}=\mu$ $=\underset{t\to 0}{lim}{\displaystyle \frac{d}{dt}}{\mathrm{e}}^{ct+t^2/4b};$

$=\underset{t\to 0}{lim}\left((c+{\displaystyle \frac{t}{2b}})\cdot {\mathrm{e}}^{ct+t^2/4b}\right);$

$=c.$

As we suspected, $c$ is equal to the mean, $\mu$. The second moment is

${{\displaystyle \frac{d^2}{d{t}^2}}M(t)|}_{t=0}$ $=\underset{t\to 0}{lim}{\displaystyle \frac{d^2}{d{t}^2}}{\mathrm{e}}^{ct+t^2/4b};$

$=\underset{t\to 0}{lim}{\displaystyle \frac{d}{dt}}\left((c+{\displaystyle \frac{t}{2b}})\cdot {\mathrm{e}}^{ct+t^2/4b}\right);$

$=\underset{t\to 0}{lim}{\displaystyle \frac{d}{dt}}(c{\mathrm{e}}^{ct+t^2/4b})+\underset{t\to 0}{lim}{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{t}{2b}}{\mathrm{e}}^{ct+t^2/4b}\right);$

$=\underset{t\to 0}{lim}\left(c(c+{\displaystyle \frac{t}{2b}})\cdot {\mathrm{e}}^{ct+t^2/4b}\right)+\underset{t\to 0}{lim}\left({\displaystyle \frac{{\mathrm{e}}^{ct+t^2/4b}}{2b}}+\left({\displaystyle \frac{t}{2b}}(c+{\displaystyle \frac{t}{2b}})\cdot {\mathrm{e}}^{ct+t^2/4b}\right)\right);$

$=c^2+{\displaystyle \frac{1}{2b}}.$

We can now find the variance using Equation 9 as

${\sigma }^2$ $={\mu }_2^{\prime }-{\mu }^2={\mu }_2^{\prime }-c^2;$

$=c^2+{\displaystyle \frac{1}{2b}}-c^2;$

$={\displaystyle \frac{1}{2b}}.$

Solving for $b$ in terms of $\sigma$ produces

$$b=\frac{1}{2{\sigma }^2}.$$

Substituting this value into Equation 22 gives us

$$a=\frac{1}{\sigma \sqrt{2\pi }}.$$

Replacing $a$, $b$, and $c$ in $f(x)=a{\mathrm{e}}^{-b{(x-c)}^2}$ with their parameterized values produces the probability density function for what is known as the Gaussian or normal distribution,

$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{\mathrm{e}}^{-\frac{{(x-\mu )}^2}{2{\sigma }^2}}.$$ (23)

Notice how Equation 23 is dependent on the mean, $\mu$, and the standard deviation, $\sigma$. As mentioned previously, changing $\mu$ just moves the center of the distribution left or right. But $\sigma$ controls both the width and the height of the distribution. When the distribution widens, it gets shorter. When the distribution narrows, it gets taller. The role of $\sigma$ is to constrain the area under the curve to remain constant independent of the value of $\sigma$. Therefore, the Gaussian probability density function represents a family of functions of unit area.

When $\mu =0$ and $\sigma =1$, the distribution corresponding to the probability density function is called the standard normal distribution,

$$f(x)=\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{-\frac{x^2}{2}}.$$ (24)

The standard normal distribution’s probability density function is used primarily as a reference function for numerically integrating probabilities that can be scaled to calculate probabilities for other normal distributions based on the value of $\sigma$. The probability within $n$ standard deviations of the mean is the same for all normal distributions regardless of the value of $\sigma$. Therefore, calculating probabilities for the standard normal distribution allows you to determine them for any normal distribution as long as you express your ranges in relative terms of standard deviations from the mean instead of using absolute numbers.

5. Commentary

The Gaussian probability density function is usually presented as a formula to be used, but not ncessarily understood. Although we attempted to show a step-by-step process from which one can get from $f(x)={\mathrm{e}}^{-x^2}$ to Equation 23, we did not explain the origin of $f(x)={\mathrm{e}}^{-x^2}$. Also, we cheated and chose $-b{(x-c)}^2$ as an exponent instead of $-b{x}^2+cx+d$ with little explanation. At this point, it should be apparent that we needed the exponent to be in the form of a square paralleling $-x^2$ so that we could integrate by substitution.

Ultimately, our goal was to show that the Gaussian probability density function did not sprout out of thin air fully formed. Instead, it evolved from a series of observations about the family of functions derived from $f(x)={\mathrm{e}}^{-x^2}$. The Gaussian integral formed the core of our exercise because we only needed to evaluate it once and were subsequently able to use it multiple times via substitution, demonstrating its utility.

The Gaussian probability distribution occurs frequently in many contexts as a result of the central limit theorem.

The Central Limit Theorem.

Given a random sample $X_{\mathrm{1}}\mathrm{,}{X}_{\mathrm{2}}\mathrm{,}\mathrm{\dots }\mathrm{,}{X}_N$ of independent and identically distributed random variables from a distribution with finite mean $\mu$ and finite variance ${\sigma }^{\mathrm{2}}$, let

$D_N={\displaystyle \frac{{\sum }_{n=1}^N{X}_n-N\mu }{\sigma \sqrt{N}}}$ $={\displaystyle \frac{\sqrt{N}}{\sigma }}(\overline{X}-\mu )\mathit{\hspace{1em}}\text{𝑤ℎ𝑒𝑟𝑒}\mathit{\hspace{1em}}\overline{X}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{n=1}^N}{X}_n.$

Then

$\underset{N\to \mathrm{\infty }}{lim}P(D_N\le d)$ $={\displaystyle {\int }_{-\mathrm{\infty }}^d}{\displaystyle \frac{1}{\sqrt{2\pi }}}{\mathrm{e}}^{-\frac{x^2}{2}}dx\mathit{\hspace{1em}}\mathit{\text{for all d.}}$

In other words, the distribution function $D_N$ converges to the standard normal distribution from Equation 24 as $N\mathrm{\to }\mathrm{\infty }$.

For large $N$—usually greater than $30$—one can pretend that $\overline{X}$ is distributed according to the normal distribution when calculating probabilities such as $P(a\le \overline{X}\le b)$. This should not be done recklessly, but does work in many situations. We will not prove the central limit theorem, but familiarity with it goes a long way to understanding why it is so common for textbooks to assume a Gaussian distribution in various contexts.

As a final comment, we explored Gaussian functions and the Gaussian probability distribution function in one dimension. The functions can can be extended to multiple dimensions, where analysis becomes considerably more complicated than what we have explored.

6. Addendum

For completeness, we have decided to demonstrate how to evaluate the integral of $f(x)=a{\mathrm{e}}^{-b{x}^2+cx+d}$ over the real numbers and summarize the integrals we covered in this paper for easy reference.

6.1. General Gaussian Integral

After evaluating the integral of the exponential function of a general concave quadratic function, it should be possible to evaluate the integral of any Gaussian function by simple parameter substitution. As before, the integral will depend on already having arrived at the result of Equation 20 and completing the square to allow integration by substitution as follows:

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-b{x}^2+cx+d}dx$ $=a{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-cx-d)}dx;$

$=a{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-cx-d)}dx\cdot {\displaystyle \frac{{\mathrm{e}}^{c^2/4b}}{{\mathrm{e}}^{c^2/4b}}};$

$=a{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-(b{x}^2-cx+c^2/4b)}dx\cdot {\mathrm{e}}^{\frac{c^2}{4b}+d};$

$=a{\mathrm{e}}^{\frac{c^2}{4b}+d}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-{(\sqrt{b}x-\frac{c}{2\sqrt{b}})}^2}dx;$

substituting $u=\sqrt{b}x-\frac{c}{2\sqrt{b}}$ and $du=\sqrt{b}dx$ to get

$=a{\mathrm{e}}^{\frac{c^2}{4b}+d}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle \frac{{\mathrm{e}}^{-u^2}}{\sqrt{b}}}du;$

$=a\sqrt{{\displaystyle \frac{\pi }{b}}}{\mathrm{e}}^{\frac{c^2}{4b}+d}.$

6.2. Summary of Gaussian Integrals

We now summarize the results of evaluating the Gaussian integrals in this paper:

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx$ $=\sqrt{\pi };$

${\displaystyle {\int }_0^{\mathrm{\infty }}}{\mathrm{e}}^{-x^2}dx$ $={\displaystyle \frac{\sqrt{\pi }}{2}};$

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-b{x}^2}dx$ $=a\sqrt{{\displaystyle \frac{\pi }{b}}};$

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-b{(x-c)}^2}dx$ $=a\sqrt{{\displaystyle \frac{\pi }{b}}};$

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-\frac{{(x-\mu )}^2}{2{\sigma }^2}}dx$ $=a\sigma \sqrt{2\pi };$

and

${\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}a{\mathrm{e}}^{-b{x}^2+cx+d}dx$ $=a\sqrt{{\displaystyle \frac{\pi }{b}}}{\mathrm{e}}^{\frac{c^2}{4b}+d}.$

All you really need is the final identity, from which all the others can be derived by substituting the appropriate values for the equation’s parameters. Note, however, that we needed to evaluate the first integral in order to evaluate all of the others.