Calculating Probability (P) given Z bounds.

Question

I need to program a simple Probability calculation function for any given Z boundaries (Area P under the normal distribution curve):

I know we can use the The Z table, but I want to actually calculate it - I found that the actual calculation is:

And for speeding up we can use Taylor Expansion:

Now I implemented it in my own code but I came across a function that looks extremely simple and produce accurate results:

 function calc_Area_Between(_zActual) {
   var p1 = 1*0;
   var p2 = 1*0;
   var pmid = 1*0;
   var z = Math.abs(_zActual);
   p2 = (((((.000005383*z+.0000488906)*z+.0000380036)*z+.0032776263)*z+.0211410061)*z+.049867347)*z+1;
   p2 = Math.pow(p2, ‐16);
   p1 = p2/2;
   pmid = 1‐p2;
   pmid = Math.round(pmid*10000)/10000;
   p2 = Math.round(p2*10000)/10000; //To sides
   p1 = Math.round(p1*10000)/10000; //One sided
   if(z>3.75) return 100;
   return pmid;
 }

Can anyone Explain What is going on there? especially the p2 calculation with those float numbers?

Answer 1

Note: No explanation but an excerpt from Appendix VII of Uncertainty, Calibration, Probability by C.F. Dietrich (see also this link)

The algorithm implements an approximation formula $P(k)$ of the Normal Distribution between the limits $-k$ to $k$.

\begin{align*} P(k)\simeq \frac{1}{\sqrt{2\pi}}\int_{-k}^{k}e^\frac{-z^2}{2}dz \end{align*}

with

\begin{align*} P(k)=1-\frac{1}{\left(1+a_1k+a_2k^2+a_3k^3+a_4k^4+a_5k^5+a_6k^6\right)^{16}} \end{align*} The coefficients $a_1$ to $a_6$ are \begin{align*} a_1&=0.04986734697\\ a_2&=0.02114100615\\ a_3&=0.00327762632\\ a_4&=0.000038{\color{blue}{{\bf 6}}}03575\\ a_5&=0.00004889{\color{blue}{\bf{6}}}635\\ a_6&=0.000005382975 \end{align*}

The maximum error $\varepsilon=\pm0.0000003$ for $0\leq k\leq 5.7$

Note the different digits in $a_4$ and $a_5$.