# Introduction and Motivation

Binomial distribution plays an important role in real world applications. It measures the probability of having $k$ successes out of $n$ i.i.d. trials with each trial having a success rate $p$. The PDF of binomial distribution with success rate $p$, total number of trials $n$ and the number of successes $k$ is
\begin{equation}
f(k,n,p)={n\choose k}p^{k}(1-p)^{n-k}
\label{bino_pdf}
\end{equation}

However, it may be difficult to directly use formula because it may contain large and small terms. Typically, ${n\choose k}$ is very large, whereas $p^k$ and $(1-p)^{n-k}$ are very small. These may cause overflow or underflow when computing them. It would be more convenient if we could reformulate the formula that avoids this issue. Is it possible to do? Maybe we can observe an example of this distribution to get some ideas.

Given the success rate $p$ of i.i.d. trials and the total number of trials $n$, the above equation gives the relationship between $k$ and $f(k,n,p)$. Below is an example of binomial distribution PDF with $n=3000$ and $p=0.3$.



The above figure looks very similar to a normal distribution, which makes intuitive sense. The peak of the distribution should coorespond to $k=np$. Furthermore, it is also very likely that $k$ can deviate around $np$ because it's very possible that we may get some more or fewer successes than expected. However, it is unlikely that the $k$ we get deviates too far from $np$. All of these intuitions can be visualized on the graph. Therefore, one may ask, is it possible to have a normal distribution to approximate the binomial distribution? It turns out the answer is yes. In this project, we will study the relationships between binomial distributions their corresponding normal approximation.

# Derivation

We begin by presenting the derivation from this paper[1].
Before we start, it is worth noting an approximation for factorials called Stirling's formula[2].
$$
\begin{aligned}
n!&=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{1}{51840n^3}+...\right)\\
&=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\left(1+O\left(\frac{1}{n}\right)\right)
\end{aligned}
$$

Equipped with this formula, we can proceed on making the approximation.
$$
\begin{aligned}
f(k,n,p)&={n\choose k}p^{k}(1-p)^{n-k}\\
&=\frac{n!}{k!(n-k)!}p^{k}(1-p)^{n-k}\\
&=\sqrt{\frac{n}{2\pi k(n-k)}}\left(\frac{np}{k}\right)^{k}\left(\frac{n(1-p)}{n-k}\right)^{n-k}\left[1+O\left(\frac1n\right)\right]
\end{aligned}
$$

Defining $\delta=k-np$, we have $k=\delta+np$ and $n-k=n(1-p)-\delta$. Expanding $\textrm{ln}(1+k)=k-\frac{1}{2}k^2+O(k^3)$, we have,
$$
\textrm{ln}\left[\left(\frac{np}{k}\right)^k\left(\frac{n(1-p)}{n-k}\right)^{n-k}\right]=-\frac{\delta^2}{2np(1-p)}+O\left(\frac{\delta^3}{n^2}\right)
$$

Exponentiating the result, we have
$$
\begin{aligned}
\left(\frac{np}{k}\right)^k\left(\frac{n(1-p)}{n-k}\right)^{n-k}&=\textrm{exp}\left(-\frac{\delta^2}{2np(1-p)}+O\left(\frac{\delta^3}{n^2}\right)\right)\\
&=\textrm{exp}\left(-\frac{\delta^2}{2np(1-p)}\right)\textrm{exp}\left(O\left(\frac{\delta^3}{n^2}\right)\right)\\
\end{aligned}
$$

Because $\delta$ is just the difference between observed $k$ and the expected value of $k$ which is $np$, $\delta$ should be on the scale of a few standard deviations, which is $\sqrt{np(1-p)}$. Therefore, we argue that $\delta\sim\sqrt{n}$.
Therefore,
$$
O\left(\frac{\delta^3}{n^2}\right)=O\left(\frac{n^{\frac32}}{n^2}\right)=O\left(\frac{1}{\sqrt{n}}\right).
$$
As $n$ gets large, $\frac{1}{\sqrt{n}}$ gets small.
Using $\textrm{exp}(x)=1+x+\frac{x^2}{2}+...=1+x+O(x^2)$, the above equation becomes
$$
\begin{aligned}
\left(\frac{np}{k}\right)^k\left(\frac{n(1-p)}{n-k}\right)^{n-k}=\textrm{exp}\left(-\frac{\delta^2}{2np(1-p)}\right)\left[1+O\left(\frac{1}{\sqrt{n}}\right)\right].
\end{aligned}
$$

For the square root part, we have
$$
\begin{aligned}
\sqrt{\frac{n}{2\pi k(n-k)}}&=\sqrt{\frac{n}{2\pi(np+\delta)(n(1-p)-\delta)}}\\
&=\sqrt{\frac{1}{2\pi np(1-p)}}\left[1+O\left(\frac{1}{\sqrt{n}}\right)\right]
\end{aligned}
$$

Multiplying them together, we obtain
$$
\begin{aligned}
f(k,n,p)=\sqrt{\frac{1}{2\pi np(1-p)}}\textrm{exp}\left(-\frac{(k-np)^2}{2np(1-p)}\right)\left[1+O\left(\frac{1}{\sqrt{n}}\right)\right]
\end{aligned}
$$

We know that for binomial distribution $f(k)$ with fixed $n$ and $p$, its mean of $k$ is $\mu=np$, and its variance is $\sigma^2=np(1-p)$. Therefore, we can rewrite the formula above as
$$
f(k)=\frac{1}{\sqrt{2\pi} \sigma}\textrm{exp}\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)\left[1+O\left(\frac{1}{\sqrt{n}}\right)\right]
$$

This is exactly a normal distribution with mean $\mu$ and variance $\sigma^2$ with a correction term $O\left(\frac{1}{\sqrt{n}}\right)$. As $n\rightarrow\infty$, $O\left(\frac{1}{\sqrt{n}}\right)$ becomes a small value compared to $1$.

There's a caveat here. we arrived at this result by assuming that $\delta\sim\sqrt{n}$, which means that this formula is applicable only when $k$ is within a few standard deviations to $np$. This formula may fail for significantly deviated $k$. This makes intuitive sense. After all, the domain of $f(k)$ is $[0,n]$ whereas the domain of normal distribution is $(-\infty,\infty)$. For $k$ smaller than $0$ or larger than $n$, our approximation returns a positive value whereas the true binomial formula returns $0$, and the ratio between them will be infinity. This is not to be worried about because our approximation disregards those exceptional cases where $\delta\sim n$ rather than $\delta\sim\sqrt{n}$.

We can now define the approximation $g(k)$ of $f(k)$ by neglecting the $O\left(\frac{1}{\sqrt{n}}\right)$ term.
$$
g(k)=\frac{1}{\sqrt{2\pi} \sigma}\textrm{exp}\left(-\frac{(k-\mu)^2}{2\sigma^2}\right).
$$

# Approximation Error Visualizations

As previously discussed, our approximation works well when $n$ is large $k$ is within a few standard deviations to $\mu$. We can visualize both the approximation and the exact function with varying $n$. Here are some examples. The x-axis limits are determined by $\mu\pm 5\sigma$ for all of the plots. The blue and red curves represent the approximations $g(k)$ and the exact functions $f(k)$, respectively.



As $n$ gets larger, the differences between two plots becomes more similar. At $n=400$, they are essentially indistinguishable. However, we may still be interested in seeing how the differences decrease with increasing $n$. Defining $d(k)=g(k)-f(k)$, we can visualize $d(k)$.



We can see that as $n$ increases the magnitudes of $d(k)$ becomes smaller, but it gets wider. Moreover, its general structures are more or less preserved. Take $n = 400$ and $n=1000$ as a comparison, if we scale y-coordinates the $n=400$ plot down, but scale up its x-coordinates, we will get some results similar to the $n=1000$ plot. In other words, if we ignore the x and y scales, the $n=400$ plot and $n=1000$ plot look pretty much identical.

We should be aware that the correction is also a function of $k$, because otherwise $d(k)$ will be a constant function. Furthermore, the observations above suggest that on the x-axis, the pattern induced by the correction lengthens in accordance to $\sigma$. For example, if $\sigma$ doubles, we should see the distance between the two local minima around the center peak of $d(k)$ also doubles. That is, the plot becomes twice as wide. For a normal distribution, if $\sigma$ doubles, we need to go $\sqrt{2}$ as far from $\mu$ to achieve the same percentile as before. Because $\sigma$ scales with $\sqrt{n}$, it implies that when $n$ is large, the contribution from the term $\frac{k}{\sqrt{n}}$ dominates those from other terms as a function $y\left(\frac{x}{\sqrt{n}}\right)$ is $\sqrt{n}$ times as wide as $y(x)$, and the disproportionalities of the graphs caused by other terms, if they exist, are visually negligible.

We may also visualize the correction terms itself by defining a new function $r(k)$ such that
$$
\begin{aligned}
r(k)&=-\frac{d(k)}{g(k)}\\
&=-\frac{\frac{-1}{\sqrt{2\pi} \sigma}\textrm{exp}\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)}{\frac{1}{\sqrt{2\pi} \sigma}\textrm{exp}\left(-\frac{(k-\mu)^2}{2\sigma^2}\right)}O\left(\frac{1}{\sqrt{n}}\right)\\
&=O\left(\frac{1}{\sqrt{n}}\right)
\end{aligned}
$$



The plots above confirms the claims that for large $n$ and $k$ not too far from $\mu$, the correction terms are small and therefore the approximations are very accurate.

An interesting observation of $r(k)$ is that its local maxima doesn't seem to be at $\mu$, but rather at some distances away from $\mu$. We can zoom in around $k=\mu$ to see more details.



It looks like at $k=\mu$, the correction $r(k)$ is negative, making $g(\mu)>f(\mu)$. As $k$ deviates from $\mu$, the difference becomes smaller and eventually $g(k)<f(k)$ for positive corrections. Then, the corrections starts decreasing which will make $g(k)>f(k)$ again.

Scale wise, while we already know that $r(k)$ roughly widens as $\sqrt{n}$, its amplitudes decrease roughly as $n^{-1}$. Take $n=100$ and $n=1000$ as a comparison, $r(\mu)$ for $n=1000$ is roughly $10$ times smaller than that for $n=100$. 

# Conclusion

We have constructed a normal approximation to a binomial distribution. The error analysis indicates that the approximation errors tends to be insignificant if $n$ is large and $k$ is not too far from $\mu$. More specifically, the corrections stays "in sync" with the approximation in a sense that the final graphs corresponding to different $n$'s are roughly proportional. Moreover, the amplitudes of the corrections decreases  roughly as $n^{-1}$.

# Generalizations

We may generalize our discussions about binomial distributions to multinomial distributions, using the same technique. In fact, the generalization has already been discussed. For a multinomial $f(\mathbf{k},n,\mathbf{p})$ such that
$$
\begin{aligned}
\mathbf{k}&=[k_1,k_2,...k_m]\\
\mathbf{p}&=[p_1,p_2,...p_m]\\
n&=\sum_{i=1}^{m}k_i\\
0&<=p_1,p_2,...,p_m<=1
\end{aligned}
$$
for varibles $\mathbf{X}=[X_1,X_2,...X_M]$, where $X_1,X_2,...,X_m$ are i.i.d. only within themselves. We can approximate this multinomial distribution PDF as
$$
\mathcal{N}(n\mathbf{p},n\mathbf{M}),
$$
where $\mathbf{M}=\mathbf{P}-\mathbf{p}\mathbf{p}^\top$
where $\mathbf{P}$ is a diagonal matrix with diagonal entries being elements of $\mathbf{p}$.

# References

[1]The Normal Approximation to the Binomial Distribution. http://scipp.ucsc.edu/~haber/ph116C/NormalApprox.pdf

[2]Stirling's approximation. In Wikipedia. https://en.wikipedia.org/wiki/Stirling%27s_approximation

[3]Geyer, C. Stat 5151 Notes: Brand Name Distributions, http://www.stat.umn.edu/geyer/5102/notes/brand.pdf
    


```julia

```

a261f86ae5426c8f532bed63c7952196849475c3
ipynb
Jupyter Notebook
_notebooks/2020-11-22-Normal-Approximation.ipynb
TianqiT/MyFastpages
e70f481e5df3e3619d890f5e57c7e985ce13dc30
Apache-2.0
_notebooks/2020-11-22-Normal-Approximation.ipynb
TianqiT/MyFastpages
e70f481e5df3e3619d890f5e57c7e985ce13dc30
Apache-2.0
2021-09-28T05:43:58.000Z
2022-02-26T10:15:46.000Z
_notebooks/2020-11-22-Normal-Approximation.ipynb
TianqiT/MyFastpages
e70f481e5df3e3619d890f5e57c7e985ce13dc30
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES