% LaTeX source for the Probability Integral

\documentclass{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{ifthen}
\usepackage{times}

\renewcommand{\Pr}{\mbox{$\mathsf P$}}
\newcommand{\E}{\mbox{$\mathsf E$}}
\newcommand{\Var}{\text{Var\ }}
\newcommand{\Cov}{\text{Cov\ }}
\newcommand{\half}{\mbox{$\frac{1}{2}$}}
\newcommand{\quarter}{\mbox{$\frac{1}{4}$}}
\newcommand{\B}{\text{B}}
\newcommand{\blankline}{\vspace{\baselineskip}\noindent}
\renewcommand{\d}{\text{d}}
\newcommand{\e}{\text{e}}
\newcommand{\gap}{\phantom{.}\hspace{3mm}}
\newcommand{\gapplus}{\phantom{.}\hspace{3mm}\phantom{1}}

\setlength{\hoffset}{-2.4cm} % 2.5cm from the left (A4 is 21.0 x 29.7cm)
\setlength{\voffset}{-2.4cm} % 2cm from the top
\setlength{\textwidth}{17cm}
\setlength{\textheight}{24cm}

\newcounter{refno}
\newcommand{\refer}{\stepcounter{refno}\newline(\therefno)\ \hangindent=1em}

\begin{document}
\thispagestyle{empty}
\begin{center}
{\large{\textbf{THE PROBABILITY INTEGRAL}}}
\end{center}

\noindent
On a donc fait un hypoth\ese, et cette hypoth\ese a \'et\'e appel\'ee
loi des erreurs.  Elle ne s'obtient pas des d\'eductions rigoreuses
\dots Tout le monde y croit cependent,'' me disait un jour M Lippman,
car les exp\'erimenteurs s'imaginent que c'est un th\'eorems de
math\'ematiques, et les math\'ematiciens que c'est un fait
exp\'erimental", H~Poincar\'e, \textit {Calcul des Probabilit\'es},
Paris: Gauthier-Villars  1896 and 1912.

\bigskip

\noindent
We say that $Z$ has a standard normal distribution if it has the
probability density function
$f_Z(z)=\phi(z)$
where $\phi(z)$ is the function
$\phi(z)=\frac{1}{\sqrt{2\pi}}\exp(-\half z^2).$
According to Gnedenko, \S22, the integral
$\int_{-\infty}^{+\infty}\phi(z)\,dz$ is called the Poisson integral.
Although this function is clearly non-negative, it is by no means clear
that it integrates to unity.  There are a number of methods of showing
that
$I=\int_{0}^{\infty}\exp(-\half z^2)\,dz=\sqrt{\frac{\pi}{2}}$
none of which is obvious.

\begin{enumerate}

\item
De Moivre showed (see A De Moivre, \textit{Approximatio ad Summam
Terminorum Binomii $\overline{a+b}\backslash^n$ in Seriem expansi}
1733, reprinted in
R~C~Archibald, A rare pamphlet of Moivre and some or his discoveries,
\textit{Isis} \textbf{8} (1926), 671--683; translated with some
A~De~Moivre, \textit{The Doctrine of Chances} (2nd edn), London:
H~Woodfall  1738, reprinted London: Cass 1967, A~De~Moivre,
\textit{The Doctrine of Chances} (3rd edn), London: A~Millar 1756,
reprinted New York, NY: Chelsea
1967 with a biographical article from \textit{Scripta Mathematica}
\textbf{2}(4) (1934), 316--333 by H~M~Walker) that if $n=2m$ and
$b(x)=\binom{2m}{x}\left(\frac{1}{2}\right)^x\left(\frac{1}{2}\right)^{2m-x}= \binom{2m}{x}2^{-2m}$
then
$b(x)\sim\frac{0.7976}{\sqrt{n}}$
and the exact value of the constant was shown by James Stirling to be
$\sqrt{2/\pi}$ (see A~Hald, \textit{A History of the Theory of
Probability and Statistics and Their Application before 1750}, New
York, NY, etc: Wiley 1990, \S24.4).  He then went on to show that
$b(m+l)\sim b(m)\exp(-2l^2/n)\sim\sqrt{\frac{2}{\pi n}}\exp(-2l^2/n).$
With $\sigma^2$ taking its binomial value $pqn=n/4$ this is of the form
$b(m+l)\sim\frac{1}{\sqrt{2}\pi\sigma^2}\exp(-\half l^2/\sigma^2) =\sigma^{-1}\phi(l/\sigma)$
but De Moivre did not use the concept of variance and did not express it
that way.  Using the binomial theorem it could be concluded that
\begin{align*}
1=\sum_{l=-m}^{+m} b(m+l) &\sim
\int_{-m}^{+m}\sigma^{-1}\phi(l/\sigma)\,dl\\
&\sim\int_{-\infty}^{+\infty}\phi(x)\,dx.
\end{align*}
It is not entirely simple to justify the limiting processes.  De Moivre
did not in fact make any remark about the integral as such.

\item
The first method by which the integral was explicitly calculated
appears to have been given by P~S~Laplace in his 1774 paper M\'emoire
sur la probabilit\'e des causes par les \'evenments.  An extract from
\noindent\dots From this we can easily conclude
$E=\frac{(p+1)\cdots(p+q+1)}{1\cdot 2\cdot 3\cdot \cdots q} .\frac{p^pq^q}{(p+q)^{p+q}}\int 2\,dz\,\cdot\exp \left(-\frac{(p+q)^3}{2pq}zz\right).$
Let $-[(p+q)^3/2pq]zz=\ln\mu$, and we will have
$\int 2\,dz\,\cdot\exp\left(-\frac{(p+q)^3}{2pq}zz\right)= -\frac{\sqrt{2qp}}{(p+q)^2}\int\frac{d\mu}{\sqrt{-\ln\mu}}.$
The number $\mu$ can here have any value between 0 and 1, and, supposing
the integral begins at $\mu=1$, we need its value at $\mu=0$.  This may
be determined using the following theorem (see M.\ Euler's Calcul
int\'egral).
Supposing the integral goes from $\mu=0$ to $\mu=1$ we have
\footnote{Cf.\ the equation
$\int_0^{\pi/2} \sin^\alpha dt=\frac{\sqrt{\pi}}{\alpha} \frac{\Gamma((1+\alpha)/2)}{\Gamma(\alpha/2)}$
(R~Courant, \textit{Differential and Integral Calculus} (2 volumes),
London and Glasgow: Blackie 1934--6, Volume II, Chapter IV, \S6,
page 338).  We use  the substitution $\mu^n=\cos\theta$ to reduce
the integrals to
$\int \frac{1}{i}\cos^{(n-i+1)/i}\theta\,d\theta\qquad\hbox{and}\qquad \int \frac{1}{i}\cos^{(n+1)/i}\theta\,d\theta.$
Essentially what it needs is that if
$I_k=\int_0^{\pi/2}\cos^k\theta\,d\theta$
then
$I_kI_{k+1}=\frac{\pi}{2(k+1)},$
which as $I_k=B((k+1)/2,1/2)/2$ is easily seen to be equivalent to
$\Gamma(1/2)^2=\pi$.
}
$\int \frac{\mu^n\,d\mu}{\sqrt{(1-\mu^{2i})}}\cdot \int \frac{\mu^{n+i}\,d\mu}{\sqrt{(1-\mu^{2i})}}= \frac{1}{i(n+1)}\cdot\frac{\pi}{2},$
whatever be $n$ and $i$.  Supposing $n=0$ and $i$ is infinitely small,
we will have $(1-\mu^{2i})/(2i)=-\ln\mu$, because the numerator and
the denominator oof this quantity become zero when $i=0$, and if we
differentiate them both, regarding $i$ alone as variable, we will have
$(1-\mu^{2i})/(2i)=\ln\mu$, therefore $1-\mu^{2i}=-2i\ln\mu$.  Under
these conditions we will thus have
$\int \frac{\mu^n\,d\mu}{\sqrt{(1-\mu^{2i})}}\cdot \int \frac{\mu^{n+i}\,d\mu}{\sqrt{(1-\mu^{2i})}}= \int \frac{d\mu}{\sqrt{2i}\sqrt{-\ln\mu}} \int \frac{d\mu}{\sqrt{2i}\sqrt{-\ln\mu}}= \frac{1}{i}\frac{\pi}{2};$
Therefore
$\int \frac{d\mu}{\sqrt{-\ln\mu}}=\sqrt{\pi}.$
Thus
$\int 2\,dz\,\cdot\exp\left(-\frac{(p+q)^3}{2pq}zz\right)= \frac{\sqrt{\phantom{2}pq}\sqrt{2\pi}}{(p+q)^{3/2}},$
from which we obtain $E=1$.''\quad Taking $p=q=1/2$ in Laplace's
result we get
$\int_{0}^{\infty}\exp(-2z^2)\,dz=\sqrt{\pi/8}$
and so on taking $x=2z$
$\int_{0}^{\infty}\exp\left(-\half x^2\right)\,dz=\sqrt{\pi/2}$
or in other words $I=\sqrt{\pi/2}$.

\item
The integral is commonly evaluated using a double integral.
The first method based on a double integral depends on noting that
$I=\int_{0}^{\infty}\exp(-\half z^2)\,dz= \int_{0}^{\infty}\exp(-\half(xy)^2)\,y\,dx$
for any $y$ (on setting $z=xy$).  Putting $z$ in place of $y$, it
follows that for any $z$
$I = \int_{0}^{\infty}\exp(-\half(zx)^2)\,z\,dx$
so that
$I^2=\left(\int_{0}^{\infty}\exp(-\half z^2)\,dz\right) \left(\int_{0}^{\infty} \exp(-\half(zx)^2)z\,dx\right) =\int_{0}^{\infty}\int_{0}^{\infty} \exp\{-\half(x^2+1)z^2\}\,z\,dz\,dx.$
Now set $(1+x^2)z^2=2t$ so that $z\,dz=dt/(1+x^2)$ to get
\begin{align*}
I^2&=\int_{0}^{\infty}\int_{0}^{\infty}
\exp(-t)\,\frac{dt}{(1+x^2)}\,dx
=\left(\int_{0}^{\infty}\exp(-t)\,dt\right)
\left(\int_{0}^{\infty}\frac{dx}{(1+x^2)}\right) \\
&=\left[-\exp(-t)\right]_0^{\infty}\,
\left[\tan^{-1}x\right]_0^{\infty}
=\bigl[1\bigr]\bigl[\half\pi\bigr] \\
&=\frac{\pi}{2}
\end{align*}
and hence $I=\sqrt{\pi/2}$ so that the integral of $\phi$ from $-\infty$
to
$\infty$ is 1, and hence $\phi$ \textit{is} a probability density
function.
This method is apparently due to P.S.~Laplace (1749--1827),
\textit{Th\'eorie Analytiques des Probabilit\'es}, \S24, pages 94--95
in the first edition.; cf.\ I~Todhunter, \textit{A History of the
Mathematical Theory of Probability from the time of Pascal to that
of Laplace}, Cambridge and London: Macmillan 1865, reprinted New York,
NY: Chelsea 1949, art.\ 899.  See, e.g., G~Valiron,  \textit{Cours
d'Analyse Math\'ematique} (2 volumes), Paris: Masson 1947, Volume I,
page 152.

\item
The usual double integral'' method is based on defining $I$ as above
and noting that
$I^2=\left(\int_0^{\infty}\exp(-\half x^2)\,dx\right) \left(\int_0^{\infty}\exp(-\half y^2)\,dy\right) =\int_0^{\infty}\!\int_0^{\infty} \exp\{-\half(x^2+y^2)\}\,dx\,dy.$
We then change to polar co-ordinates $(r, \theta)$ in which
$dx\,dy=r\,dr\,d\theta$, so that
\begin{align*}
I^2&=\int_0^{\pi/2}\!\int_0^{\infty}
\exp\left(-\half r^2\right)\,r\,dr\,d\theta
=\left(\int_0^{\pi/2}d\theta\right)
\left(\int_0^{\infty}\exp\left(-\half
r^2\right)\,r\,dr\right)
=\bigl[\theta\bigr]_0^{\pi/2}\,
\bigl[-\exp\left(-\half r^2\right)\bigr]_{0}^{\infty}
\\
&=\frac{\pi}{2}
\end{align*}
from which it follows that $I=\sqrt{\pi/2}$ so that the integral of
$\phi$ from
$-\infty$ to $\infty$ is 1, and hence $\phi$ \textit{is} a probability
density function.  This method is apparently due to Sim\'eon Denis
Poisson (1781--1840) and was popularized by Jacob Karl Franz Sturm
(1803--1855)---see his \textit{Cours d'Analyse de l'\'ecole
polytechnique}, Paris: Mallet-Bachelier, Volume 2, pages 16--17

466.  L'integrale
$A=\int_0^{\infty} e^{-x^2}dx$
a \'et\'e d\'etermin\'ee par M.\ Poisson \a l'aide d'un proc\'ed\'e
tr\es-remarquable.  Si l'on change $x$ en $y$, on aura encore
$A=\int_0^{\infty} e^{-y^2}dy,$
et, par suite,
$A^2=\int_0^{\infty} e^{-x^2}dx.\int_0^{\infty} e^{-y^2}dy =\int_0^{\infty}\int_0^{\infty}e^{-x^2-y^2}dxdy.$

Soient mantainent trois axes rectangulaire O$x$, O$y$, O$z$ et
$y=0,\quad z=e^{-x^2},$
les equations d'une courbe situ\'ee dans le plan $z$O$x$.  Si cette
courbe tourne autour de l'axe O$z$, elle engendrera une surface ayant
pour \'equation
$z=e^{-x^2-y^2},$
et l'int\'egrale double
$\int_0^{\infty}\int_0^{\infty}e^{-x^2-y^2}dxdy,$
r\'epr\'esentera le quart du volume compris entre le surface et le plan
$x$O$y$.  On peut \'evaluer ce volume par le partageant en une
infinit\'e de tranches cylindriques dont O$z$ soit l'axe commun.  La
tranche dont les surfaces ext\'erieurs ont pour rayone $r$ et $r+dr$ est
\'egale \a sa base $2\pi rdr$ multipli\'ee par sa hauteur $z$ on
$e^{-r^2}$: on a donc
$A^2=\frac{1}{4}\int_0^{\infty}e^{-r^2}\times 2\pi rdr=\frac{1}{4}\pi$
d'o\u
$A=\frac{1}{2}\sqrt{\pi}.\text{''}$

\item
Another method comes from the fact that
$\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin\pi z}$
with $z=\half$---see E~T~Whittaker and G~N~Watson, \textit{A Course of
Modern Analysis}, Cambridge University Press 1902, 1915, 1920 and 1927.
\S12.14, J~C~Burkill and H~Burkill, \textit{A Second Course in
Mathematical Analysis}, Cambridge University Press 1970, \S14.6, or
E~T~Copson, \textit{Theory of Functions of a Complex Variable}, Oxford:
Clarendon Press 1935, \S9.22.

\item
Yet another method results from substituting $u=\exp(-\half z^2)$,
giving
$I = \int_0^1 \frac{du}{\sqrt{-2\ln u}}.$
Now note that, for $x > 0$,
$\frac{x-1}{x} \leqslant \ln x \leqslant x - 1$
which follows geometrically from the convexity of the logarithmic
function, or can be easily established using calculus to show, for
example, that $x - \ln x$ has smallest value 1.

For $0 < v < 1$ and for any positive integer $n$, write $v_n=v^{1/n}$,
so that $\ln v = n\ln v_n$.  From the above inequalities with $x=v_n$,
$n\left(\frac{v_n-1}{v_n}\right) \leqslant \ln v \leqslant n(v_n-1)$
from which
$\frac{1}{\sqrt{2n}}\sqrt{\frac{v_n}{1-v_n}} \leqslant \frac{1}{\sqrt{-2\ln v}} \leqslant \frac{1}{\sqrt{2n}}\frac{1}{\sqrt{1-v_n}}.$
Integrating these inequalities between 0 and 1, we obtain
$J_n \leqslant I \leqslant K_n$
where
$J_n = \frac{1}{\sqrt{2n}} \int_0^1 \sqrt{\frac{v_n}{1-v_n}}dv \qquad\text{and}\qquad K_n = \frac{1}{\sqrt{2n}} \int_0^1 \frac{1}{1-v_n} dv.$
Now substitute $v_n=\sin^2\phi$, i.e.\ $v=\sin^{2n}\phi$.  Then
$J_n = \sqrt{2n}\int_0^{\pi/2} \sin^{2n}\phi d\phi \qquad\text{and}\qquad K_n = \sqrt{2n}\int_0^{\pi/2} \sin^{2n-1}\phi d\phi.$
It is thus clear (as integrals of powers of $\sin\phi$ must decrease
with the power involved) that
$0 \leqslant \sqrt{\frac{n}{n+1}}K_{n+1} \leqslant J_n \leqslant K_n.$
Furthermore, by the usual reduction method
$\frac{J_{n+1}}{J_n} = \frac{K_n}{K_{n+1}} = \frac{2n+1}{2\sqrt{n(n+1)}} > 1$
so that
$J_{n+1}K_{n+1} = J_nK_n = \dots = J_1K_1 = \pi/2.$
It follows that, as $n\to\infty$, $K_n$ decreases and $J_n$ increases
to a common limit $\sqrt{\pi/2}$.  It follows that as $J_n\leqslant I\leqslant K_n$, we have $I=\sqrt{\pi/2}$.  This method can be found in
N~Gauthier, Note 72.22 Evaluating the probability integral,
\textit{Mathematical Gazette}, \textbf{72} (1988), 124--125, and
D~Desbrow, Note 74.28 Evaluating the probability integral,
\textit{Mathematical Gazette} \textbf{74} (1990), 169--170,
but I am not sure whether it originated there.

\item
It was long supposed that the integral could not be evaluated by the
Cauchy method of residues, but it turns out that it can be (see, e.g.,
J~C~Burkill and H~Burkill, \textit{A Second Course in Mathematical
Analysis}, Cambridge University Press 1970, Exercises 14(a), no.\ 15).
This method depends on setting
$f(z) = \exp\left(\pi i z^2\right)/\sin(\pi z)$
which function has residue $1/\pi$ at $z=0$.  Then since
$\sin\{\pi(z-1)\} = -\sin\{\pi z\}$
and
$z^2=z(z-1)+z\qquad\text{and}\qquad(z-1)^2=z(z-1)-z+1$
we see that
$f(z) - f(z-1) = 2i\exp\{\pi i z(z-1)\}.$
By integrating $f$ round a parallelogram with vertices
$\pm\half\pm R\exp(\quarter\pi i)$, where $R$ is large, (putting
$z=t\exp(i\pi/4)$ $\pm\half$ so $dz=dt\exp(i\pi/4)$ along long sides) we
see that
$\int_{-\infty}^{\infty} \exp(-\pi t^2) dt = 1.$
This method was due to G.\ P\'olya (1949).  A previous method using
contour integrals due to L.J.\ Mordell (1920) contains [the
probability integral] as a special case, [but] the methods used by
Mordell are too complicated and it is not really worthwhile applying
them to [this case]''.  Another method is due to J.H.\ Cadwell (1947).
For more details, see D~S~Mitrinovi\'c and J~D~Ke\v cki\'c,
\textit{The Cauchy Method of Residues: Theory and Applications},
Dordrecht, etc: Reidel 1984, \S5.3.4.10, pp.\ 158--168.

\item Recently T~P~Jameson (1994) (at the age of 16!) has suggested yet
another method (subsequently suggested independently by S~P~Eveson (2005)).
Consider the volume under the surface $z=\text{e}^{-(x^2+y^2)}$,
which is clearly given by
\begin{align*}
V &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
\text{e}^{-(x^2+y^2)}\text{d}x\text{d}y \\
&= \left(\int_{-\infty}^{\infty}
\text{e}^{-x^2}\text{d}x\right)^2.
\end{align*}
This can, however, also be thought of as a volume of revolution about
the $z$ axis where as $z=\text{e}^{-x^2}$ we have $x=\sqrt{-\log z}$.
Using the standard formula for a volume of revolution
$V = \pi\int_0^1 x^2\text{d}z = \pi\int_0^1\{-\log z\}\,\text{d}z = \pi[-z\log z+z]_0^1 = \pi.$
and hence
$\int_{-\infty}^{\infty}\text{e}^{-x^2}\text{d}x = \sqrt{\pi}.$

\end{enumerate}

\bigskip
\noindent {\Large{\textit{References}}}
\normalsize
\ \\
\refer R~C~Archibald, A rare pamphlet of Moivre and some or his
discoveries, \textit{Isis} \textbf{8} (1926), 671--683.
\refer J C Burkill and H Burkill, \textit{A Second Course in
Mathematical Analysis}, Cambridge: University Press 1970.
\refer R Courant, \textit{Differential and Integral Calculus} (2
volumes), London and Glasgow: Blackie 1934--6.
\refer Copson, \textit{Theory of Functions of a Complex Variable},
Oxford: Clarendon Press 1935.
\refer A De Moivre, \textit{Approximatio ad Summam Terminorum Binomii
$\overline{a+b}\backslash^n$ in Seriem expansi} 1733, reprint\-ed
in Archibald (1929); translated with some additions in De Moivre
(1738), De Moivre (1756) and Smith (1929).
\refer A De Moivre, \textit{The Doctrine of Chances} (2nd ed.), London:
H Woodfall 1738, reprinted London: Cass 1967.
\refer A De Moivre, \textit{The Doctrine of Chances} (3rd ed.),
London: A Millar 1756, reprinted New York, NY: Chelsea 1967
with a biographical article from \textit{Scripta Mathematica}
\textbf{2}(4) (1934), 316--333 by H M Walker.
\refer D~Desbrow, Note 74.28 Evaluating the probability integral,
\textit{Mathematical Gazette} \textbf{74} (1990), 169--170,
\refer L Euler, \textit{Calcul int\'egral}, translation of
\textit{ Institvtionvm calcvli integralis}, Petropoli: Impensis
\refer S P Eveson, Private communication (2005).
\refer T P Jameson, \textit{Mathematical Gazette} \textbf{78} (1994),
note 78.16, pp. 339--340.
\refer N~Gauthier, Note 72.22 Evaluating the probability integral,
\textit{Mathematical Gazette}, \textbf{72} (1988), 124--125.
\refer B V Gnedenko, \textit{Theory of Probability}, Moscow: Nauka 1954
and 1961, English translation Moscow: Mir 1969 and New York, NY:
Chelsea 1967.
\refer A Hald, \textit{A History of the Theory of Probability and
Statistics and Their Application before 1750}, New York, NY, etc:
Wiley 1990.
\refer P S Laplace, M\'emoire sur la probabilit\'e des causes par
les \'evenments, \textit{M\'emoires de mathe\'ematique et de
physique present\'es \a l'Acad\'emie royale des sciences, par
divers savans, \& lus dans ses assembl\'ees} \textbf{6}, 621--626,
reprinted in Laplace's \textit{Oeuvres compl\etes} \textbf{8},
27--65, translated with an introduction by S M Stigler,
\textit{Statistical Science} \textbf{1} (1986), 359--378.
\refer P S Laplace, \textit{Th\'eorie Analytiques des Probabilit\'es},
Paris: Courcier 1812, reprinted Bruxelles: Culture et Civilisation
1967.
\refer D S Mitrinovi\'c and J D Ke\v cki\'c, \textit{The Cauchy Method
of Residues: Theory and Applications}, Dordrecht, etc: D Reidel
1984.
\refer D E Smith, \textit{A Source Book in Mathematics} (2 volumes),
New York, NY: McGraw-Hill 1929, reprinted New York, NY: Dover 1959.
\refer J K F Sturm, \textit{Cours d'Analyse de l'\'ecole
polytechnique}, Paris: Mallet-Bachelier, 1857.
\refer E C Titchmarsh, \textit{Theory of Functions}, Oxford: University
Press 1932 and 1939.
\refer I Todhunter, \textit{A History of the Mathematical Theory of
Probability from the time of Pascal to that of Laplace}, Cambridge
and London: Macmillan 1865, reprinted New York, NY: Chelsea 1949.
\refer G Valiron, \textit{Cours d'Analyse Math\'ematique} (2 volumes),
Paris: Masson 1947.
\refer E T Whittaker and G N Watson, \textit{A Course of Modern
Analysis}, Cambridge: University Press 1902, 1915, 1920 and 1927.

\begin{flushright}
P.M.L.
\end{flushright}

\end{document}

% `