\documentclass{article}
\usepackage{fullpage}[]
\usepackage{cancel}[]
\usepackage{amssymb,amsmath}
%\usepackage[margin=2cm]{geometry}
\begin{document}
  
\noindent  Homework 1 -- February 1, 2006 \hfill \textbf{Solution} prepared by Tobin Fricke\\
Physics 402 -- Probability  \hfill  Prof. S.G. Rajeev, Univ. of Rochester \\

\noindent\textbf{1. Moments of the Poisson distribution.} The mean of a discrete random variable $\xi$ is calculated as $\langle \xi\rangle = \sum k P[\xi=k]$ where the sum is over all possible values $k$ the random variable may attain.  For a random variable $\xi$ with Poisson distribution with parameter $\lambda$, the probability mass function is $P[\xi=k]=(\lambda^k/k!)\exp(k)$.

$$
\langle \xi\rangle 
= \sum_{k=0}^\infty k \frac{\lambda^k}{k!} e^{-\lambda}
= e^{-\lambda} \sum_{k=1}^\infty k \frac{\lambda^k}{k!}
= e^{-\lambda} \lambda \sum_{k=1}^\infty \frac{\lambda^{k-1}}{(k-1)!}
= e^{-\lambda} \lambda \sum_{k'=0}^\infty \frac{\lambda^{k'}}{k'!}
= e^{-\lambda}e^\lambda \lambda = \boxed{\lambda}
$$

Instead of directly calculating $\sigma^2=\left\langle\left(\xi-\langle\xi\rangle\right)^2\right\rangle = \langle\xi^2\rangle-\langle\xi\rangle^2$ or even $\langle\xi^2\rangle$, it is easiest to first compute $\langle\xi(\xi-1)\rangle = \langle\xi^2\rangle-\langle\xi\rangle$, following a trick similar to the above:

$$
\langle \xi(\xi-1)\rangle 
= \sum_{k=0}^\infty k(k-1) \frac{\lambda^k}{k!} e^{-\lambda}
%= e^{-\lambda} \sum_{k=2}^\infty k(k-1) \frac{\lambda^k}{k!}
= e^{-\lambda} \lambda^2 \sum_{k=2}^\infty \frac{\lambda^{k-2}}{(k-2)!}
= e^{-\lambda} \lambda^2 \sum_{k'=0}^\infty \frac{\lambda^{k'}}{k'!}
= e^{-\lambda}e^\lambda \lambda^2 = \lambda^2
$$

Using $\langle\xi(\xi-1)\rangle=\lambda^2$ and $\langle\xi\rangle=\lambda$ we can find the variance $\sigma^2$:

$$
\sigma^2 = \langle\xi^2\rangle-\langle\xi\rangle^2 = \langle\xi(\xi-1)\rangle + \langle\xi\rangle - \langle\xi\rangle^2 = \lambda^2 + \lambda - \lambda^2 = \lambda
\implies\boxed{\sigma = \sqrt{\lambda}}$$

The same trick works for the computation of $\langle\xi^3\rangle$ (the ``third moment'' of the distribution):

$$
\langle \xi(\xi-1)(\xi-2)\rangle 
= \sum_{k=0}^\infty k(k-1)(k-2) \frac{\lambda^k}{k!} e^{-\lambda}
%= e^{-\lambda} \sum_{k=3}^\infty k(k-1)(k-2) \frac{\lambda^k}{k!}
= e^{-\lambda} \lambda^3 \sum_{k=2}^\infty \frac{\lambda^{k-3}}{(k-3)!}
= e^{-\lambda} \lambda^3 \sum_{k'=0}^\infty \frac{\lambda^{k'}}{k'!}
%= e^{-\lambda}e^\lambda \lambda^3 
= \lambda^3
$$
%\langle\xi(\xi-1)(\xi-2)\rangle
% = (\xi^2 - \xi)(\xi - 2) 
% = (\xi^3 - \xi^2 - 2\xi^2 + 2\xi)
%= \langle\xi^3 - 3\xi^2 + 2\xi\rangle$
$$
\left\langle\xi^3\right\rangle 
= \langle\xi(\xi-1)(\xi-2)\rangle + 3\langle\xi^2\rangle - 2\langle\xi\rangle
= \lambda^3 + 3\left(\sigma^2 + \langle\xi\rangle^2\right) - 2\langle\xi\rangle
= \lambda^3 + 3\left(\lambda + \lambda^2\right) - 2\lambda
= \boxed{\lambda^3 + 3\lambda^2 + \lambda}
$$

\noindent\textbf{Aside:} We will see later that, for a random variable $\xi$, we can define a function $f_\xi(t)$ called $\xi$'s \emph{characteristic function} that is the Fourier transform of its probability distribution, $f_\xi(t)=\langle\exp\{i t\xi\}\rangle$. Given this function, it turns out that we may easily compute the $n$th moment $\langle\xi^n\rangle$ of $\xi$ as $\langle\xi^n\rangle=(1/i^n)(\partial^n/\partial t^n)f_\xi(t)|_{t=0}$.   For a random variable $\xi$ with Poisson distribution, one can find $f_\xi(t)=\exp\left\{\lambda\left(\exp(i t)-1\right)\right\}$.  The following Mathematica code computes the characteristic function and uses it to generate a table of moments:

\begin{verbatim}
p[k_] := Exp[-\[Lambda]] \[Lambda]^k / Factorial[k] 

f[t_] = Sum[Exp[I k t] p[k], {k, 0, Infinity}]

Table[{\[LeftAngleBracket]\[Xi]^n\[RightAngleBracket],
      I^(-n) D[f[t], {t, n}] /. t -> 0 // Expand}, {n, 1, 4}] // TableForm
\end{verbatim}

\vspace{1cm}\pagebreak\noindent \textbf{2.  Convergence of the binomial distribution to the Poisson distribution.} 
For a random variable $\xi$ with Binomial distribution, given $N$ trials and probabilities $p$ and $q=(1-p)$ for success and failure, the probability of exactly $k$ successes is:

$$ P[\xi=k] = {N \choose k} p^k q^{N-k} $$

$$ P[\xi=k] = \frac{N!}{(N-k)! k!} p^k (1-p)^{N-k} $$

$$ P[\xi=k] =\frac{N (N-1) (N-2) \cdots (N-k+1)}{k!}p^k (1-p)^{N-k} $$

Notice that there are $k$ terms in the denominator of the leading fraction.  Divide the top and bottom by $N^k$:

$$ P[\xi=k] = N^k \frac{N}{N} \frac{N-1}{N} \frac{N-2}{N} \cdots \frac{N-k+1}{N}\frac{1}{k!}p^k (1-p)^{N-k} $$

$$ P[\xi=k] = \frac{N^k p^k}{k!} \left(1-\frac{1}{N}\right) \left(1-\frac{2}{N}\right) \cdots \left(1-\frac{k-1}{N}\right) (1-p)^{N-k} $$

We write $p$ in terms of the number of expected successes, $p = \lambda/N$:

$$ P[\xi=k] = \frac{\lambda^k}{k!} \left(1-\frac{1}{N}\right) \left(1-\frac{2}{N}\right) \cdots \left(1-\frac{k-1}{N}\right) \left(1-\frac{\lambda}{N}\right)^{N-k} $$

To show convergence to the Poisson distribution, we take the limit $N\to\infty$.  We'll need the identity $\lim_{N\to\infty}\left(1-\lambda/N\right)^N\to\exp(-\lambda)$, which we may derive:

$$ 
\lim_{N\to\infty}\left(1+\frac{\lambda}{N}\right)^N = 
\lim_{N\to\infty}\sum_{k=0}^{N}{N \choose k}\left(\frac{\lambda}{N}\right)^k 
=\lim_{N\to\infty}\sum_{k=0}^N\left(\prod_{i=0}^{k-1}\frac{N-i}{N}\right)\frac{\lambda^k}{k!}
=\sum_{k=0}^\infty \frac{\lambda^k}{k!} = e^\lambda$$

Our expression for $P[\xi=k]$ was:

$$ P[\xi=k] = \frac{\lambda^k e^{-\lambda}}{k!} \left(1-\frac{1}{N}\right) \left(1-\frac{2}{N}\right) \cdots \left(1-\frac{k-1}{N}\right) \left(1-\frac{\lambda}{N}\right)^{N-k}$$

In the limit $N\to\infty$ all of the $(1-\cdot/N)$ terms will go to unity and the final term will go to $\exp(-\lambda)$ as we just saw, and we obtain the familiar Poisson distribution:

$$ \lim_{N\to\infty,Np\to\lambda}P[\xi=k] = \frac{\lambda^k e^{-\lambda}}{k!} $$

\noindent We are, however, interested in quantifying the error when we use the Poisson distribution to approximate a Binomial distribution with ``large'' $N$.  First we must also expand the $(1-\lambda/N)^N$ term in powers of N:

$$
\left(1-\frac{\lambda}{N}\right)^{N-k}
=\exp\left\{(N-k)\log\left[1-\frac{\lambda}{N}\right]\right\}
$$

Expand $\log(1-x)$ in a power series about $x_0=1$:

$$ \log(1-x) = - \sum_{j=1}^{\infty} \frac{x^j}{j} = - x - \frac{x^2}{2} - \frac{x^3}{3} - \cdots$$


Using that expansion:

$$
\left(1-\frac{\lambda}{N}\right)^{N-k}
% =\exp\left\{(N-k)\log\left[1-{\lambda\over N}\right]\right\}
=\exp\left\{-(N-k)\sum_{j=1}^\infty\frac{1}{j}\left(\frac{\lambda}{N}\right)^j\right\}
%=\exp\left\{-(N-k)\left( \frac{\lambda}{N} + \frac{1}{2} \left(\frac{\lambda}{N}\right)^2 + \frac{1}{3} \left(\frac{\lambda}{N}\right)^3 + \cdots \right)\right\}
%=\exp\left\{\lambda - \frac
%-(N-k)\left( \frac{\lambda}{N} - \frac{1}{2} \left(\frac{\lambda}{N}\right)^2 + \frac{1}{3} \left(\frac{\lambda}{N}\right)^3 + \cdots \right)\right\}
$$

Working out the term in curly braces:
$$
-(N-k)\sum_{j=1}^\infty\frac{1}{j}\left(\frac{\lambda}{N}\right)^j
= -\sum_{j'=0}^\infty \frac{1}{j'+1} \frac{\lambda^{j'+1}}{N^{j'}} +\sum_{j=1}^\infty \frac{k}{j}\frac{\lambda^j}{N^j}
= -\lambda - \sum_{j=1}^\infty \left(\frac{\lambda}{N}\right)^j \left(\frac{\lambda}{j+1} - \frac{k}{j}\right)
$$

$$
%-(N-k)\sum_{j=1}^\infty\frac{1}{j}\left(\frac{\lambda}{N}\right)^j
= -\lambda  + \frac{k\,\lambda }{N} + \frac{k\,{\lambda }^2}{2\,N^2} - \frac{{\lambda }^2}{2\,N} - \frac{{\lambda }^3}{3\,N^2} + \cdots
%= -\lambda \left( -1 + \frac{k}{N} + \frac{k\,\lambda }{2\,N^2} - \frac{\lambda }{2\,N} -   \frac{{\lambda }^2}{3\,N^2} + \cdots\right)
$$

$$
\left(1-\frac{\lambda}{N}\right)^{N-k}
= \exp(-\lambda) \exp\left\{  \frac{k\,\lambda }{N} + \frac{k\,{\lambda }^2}{2\,N^2} - \frac{{\lambda }^2}{2\,N} - \frac{{\lambda }^3}{3\,N^2} + \cdots\right\}
$$

Putting everything together, we have:

$$ P[\xi=k] = \frac{\lambda^k e^{-\lambda}}{k!} 
%\left(1-\frac{1}{N}\right) \left(1-\frac{2}{N}\right) \cdots \left(1-\frac{k-1}{N}\right) 
\left( \prod_{i=0}^{k-1}\left(1-\frac{i}{N}\right)\right)
\exp\left\{  \frac{k\,\lambda }{N} + \frac{k\,{\lambda }^2}{2\,N^2} - \frac{{\lambda }^2}{2\,N} - \frac{{\lambda }^3}{3\,N^2} + \cdots\right\}
$$

We can expand the exponential, too, as a power series about 0, and we may group the series by powers of $1/N$.  Expanding the product in the middle expression and multiplying by this power series will result in another series in powers of $1/N$.  To first order, we will see that the estimate is correct.   The next term, in $1/N$ is the first error term.  Therefore the error of our estimate is of order $1/N$.

\vspace{1cm}\pagebreak
\noindent \textbf{3. The Stirling Approximation.}  
The Stirling approximation (in one form) is:

$$\boxed{n! \approx n^n e^{-n} \sqrt{n} \sqrt{2 \pi}}$$

This bit of Mathematica produces a decent table:

\begin{verbatim}
stirling[n_] := n^n Exp[-n] Sqrt[n] Sqrt[2 Pi]

Table[{n, N[Factorial[n]], N[stirling[n]], 
  N[stirling[n]/Factorial[n] - 1]}, {n, 1, 20}] // TableForm
\end{verbatim}

$$
\begin{tabular}{|*{1}{*{4}{|l}|}|}
\hline 
$n$ & $n!$ & $n^n e^{-n} \sqrt{n} \sqrt{2 \pi}$ & relative error \\
\hline
$1$ & $1.$ & $0.922137$ & $-0.077863$ \\
$2$ & $2.$ & $1.919$ & $-0.0404978$ \\
$3$ & $6.$ & $5.83621$ & $-0.0272984$ \\
$4$ & $24.$ & $23.5062$ & $-0.020576$ \\
$5$ & $120.$ & $118.019$ & $-0.0165069$ \\
$6$ & $720.$ & $710.078$ & $-0.0137803$ \\
$7$ & $5040.$ & $4980.4$ & $-0.0118262$ \\
$8$ & $40320.$ & $39902.4$ & $-0.0103573$ \\
$9$ & $362880.$ & $359537.$ & $-0.00921276$ \\
$10$ & $3.6288\times{10}^6$ & $3.5987\times{10}^6$ & $-0.00829596$ \\
$11$ & $3.99168\times{10}^7$ & $3.96156\times{10}^7$ & $-0.00754507$ \\
$12$ & $4.79002\times{10}^8$ & $4.75687\times{10}^8$ & $-0.00691879$ \\
$13$ & $6.22702\times{10}^9$ & $6.18724\times{10}^9$ & $-0.0063885$ \\
$14$ & $8.71783\times{10}^{10}$ & $8.6661\times{10}^{10}$ & $-0.0059337$ \\
$15$ & $1.30767\times{10}^{12}$ & $1.30043\times{10}^{12}$ & $-0.00553933$ \\
$16$ & $2.09228\times{10}^{13}$ & $2.08141\times{10}^{13}$ & $-0.00519412$ \\
$17$ & $3.55687\times{10}^{14}$ & $3.53948\times{10}^{14}$ & $-0.0048894$ \\
$18$ & $6.40237\times{10}^{15}$ & $6.3728\times{10}^{15}$ & $-0.00461846$ \\
$19$ & $1.21645\times{10}^{17}$ & $1.21113\times{10}^{17}$ & $-0.00437596$ \\
$20$ & $2.4329\times{10}^{18}$ & $2.42279\times{10}^{18}$ & $-0.00415765$ \\
\hline 
\end{tabular}
$$

This form of the Stirling approximation is very accurate.  Its accuracy is better than $10\%$ for $n\ge 1$ and better than $1\%$ for $n\ge 9$.  The worst error, at $n=1$ is about $7.8\%$.  In this form, the approximation is indeterminate for $n=0$, since the approximation involves an $n^n$ term.

In thermodynamics and statistical mechanics, it is common to use a `Stirling' approximation for $\log (n!)$:

$$\boxed{\log(n!)\approx (n\log n)-n}$$.
\end{document}
