Why is binomial probability calculation wrong in pgfplots?

Question

During my investigation of this nice question by @Sentient, I found a rather bizarre problem: Given the binomial parameters n=15 and p=0.7, @Sentient’s code produces some probability values greater than two! This can be illustrated with other parameters as well:

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
    \addplot[cyan, domain=0:1, samples=26, smooth]{normd(x, 25, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 25, 0.72)};
\end{axis}
\end{tikzpicture}\quad
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
    \addplot[cyan, domain=0:1, samples=27, smooth]{normd(x, 26, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 26, 0.72)};
\end{axis}
\end{tikzpicture}
\end{document}

On the left, I set n=25 and p=0.72; On the right, I set n=26 and p=0.72. The orange curves are the Gaussian approximations generated by the function normaldensity (the formula can be found in my previous answer), which should be close to the cyan curves (generated by normd, “normalized binomial”) in both graphs.

However, the calculations of binomial probabilities by normd are way off on the right, producing probability values bigger than 3! On the contrary, the calculations on the left are perfectly fine. I’m sure @marmot has already noticed this problem when writing this excellent answer, because the gamma function approach produces curves very different from @Sentient’s (but are very close to the Gaussian curves).

My question is: Why does this calculation error occur (but only sometimes, e.g., n = 23, 24 or 26 but not for n = 25), and how can I fix it?

Answer 1

I forgot my own take on how TeX stores real numbers (how embarrassing!). Long story short, in the original function normd by @Sentient, \x*\n is passed to be put into factorial calculations. Even after I told PGFplots to use

\x = 0, 1/\n, 2/\n, ... , \n/\n

it is very likely that the resulting \x*\n’s are none of the integers from 1 to n-1. This leads to ridiculous results of the factorials, and then leads to probabilities greater than one.

So, a clear fix is to use integer \x, then normalizing the x-coordinate as @marmot tried in the comments. This method at least guarantees the factorials are somewhat accurate. The problem of overflow is another beast which I don’t want to bring up.

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\section*{The correct left is by luck}
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
    \addplot[cyan, domain=0:1, samples=26, smooth]{normd(x, 25, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 25, 0.72)};
    \addplot[red,only marks,samples at={0,1,...,25}](x/25,{binom(x, 25, 0.72)});
\end{axis}
\end{tikzpicture}\quad
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
    \addplot[cyan, domain=0:1, samples=27, smooth]{normd(x, 26, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 26, 0.72)};
    \addplot[red,only marks,samples at={0,1,...,26}](x/26,{binom(x, 26, 0.72)});
\end{axis}
\end{tikzpicture}
\section*{One should use integers for factorial calculations}
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
%    \addplot[cyan, domain=0:1, samples=26, smooth]{normd(x, 25, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 25, 0.72)};
    \addplot[red,only marks,samples at={0,1,...,25}](x/25,{binom(x, 25, 0.72)});
\end{axis}
\end{tikzpicture}\quad
\begin{tikzpicture}[
    scale=0.88,
    declare function={binom(\x,\n,\p)=\n!/(\x!*(\n-\x)!)*\p^\x*(1-\p)^(\n-\x);},
    declare function={normd(\x,\n,\p)=binom(\x*\n,\n,\p);},
    declare function={normaldensity(\x,\n,\p)=exp(-\n*(\x-\p)^2/(2*\p*(1-\p)))/sqrt(2*pi*\n*\p*(1-\p));}
]
\begin{axis}
%    \addplot[cyan, domain=0:1, samples=27, smooth]{normd(x, 26, 0.72)};
    \addplot[orange, domain=0:1, smooth]{normaldensity(x, 26, 0.72)};
    \addplot[red,only marks,samples at={0,1,...,26}](x/26,{binom(x, 26, 0.72)});
\end{axis}
\end{tikzpicture}
\end{document}

As someone who works with probabilities a lot, I barely used PGFplots to program my distribution functions. I found it easier (and more accurate) to use data generated by statistics software to draw in PGFplots.