Normal Distribution formula for Percentile point function

Question

I am looking for a formula that allows me to calculate the Z value of the normal distribution acumulative for example:

I have a value ${\bf \alpha} = 0.975$

and in the table the ${\bf Z} = 1.960 $

in a nutshell I have the percentage value of $\alpha$ and my goal is to find Z. In python there is a library that allows me to do this.

from scipy.stats import norm as zeta
alpha = 0.95
rs = zeta.ppf(vara)
print(rs)

Information scipy

Answer 1

According to the documentation, help(zeta):

Percent point function (inverse of cdf) at q of the given RV.

So you want to invert the cumulative distribution function

$$ F(a) = \int_{-\infty}^a \frac{1}{\sqrt{2\pi}} \, \exp\left(-\dfrac{x^2}{2} \right) dx.$$

Say you want to find $a$ with $F(a) = 0.95$:

$$ 0.95 = \int_{-\infty}^a \frac{1}{\sqrt{2\pi}} \, \exp\left(-\dfrac{x^2}{2} \right) dx.$$

Unfortunately, there is not an analytical solution, only a numerical solution.

Answer 2

As already said, for a rigorous calculation of $a$, you will need some numerical method.

However, we can have quite good approximation since $$F(a) = \int_{-\infty}^a \frac{1}{\sqrt{2\pi}} \, \exp\left(-\dfrac{x^2}{2} \right) dx=\frac{1}{2} \left(1+\text{erf}\left(\frac{a}{\sqrt{2}}\right)\right)$$

Rewrite is as $$(2F(a)-1)^2=\left(\text{erf}\left(\frac{a}{\sqrt{2}}\right)\right)^2$$ and have a look here where you will see good approximations $$\left(\text{erf}\left(x)\right)\right)^2=1-e^{-k x^2}$$ where $k=\frac{\pi^2} 8$ or, (supposed to be slightly better) $k=(1+\pi )^{2/3} \log ^2(2)$. This gives as an approximation $$a=\sqrt{-\frac{2 \log [4 (1-F(a)) F(a)]}{k} }$$ Appplied to $F(a)=0.95$, the first $k$ would give $$a=\frac{4 \sqrt{\log \left(\frac{100}{19}\right)}}{\pi } \approx 1.64082$$ and the second $$a=\frac{\sqrt{2 \log \left(\frac{100}{19}\right)}}{\sqrt[3]{1+\pi } \log (2)} \approx 1.63726$$ while the "exact" solution would be $1.64485$

A bit more complex (but this is just a quadratic equation in $x^2$) would use $$\mathrm{erf}\!\left(x\right)^2\approx1-\exp\Big(-\frac 4 {\pi}\,\frac{1+\alpha x^2}{1+\beta x^2}\,x^2 \Big)$$ where $$\alpha=\frac{10-\pi ^2}{5 (\pi -3) \pi }\qquad \text{and} \qquad \beta=\frac{120-60 \pi +7 \pi ^2}{15 (\pi -3) \pi }$$

Applied to the worked example, this last formula would give $a=1.64528$.

Edit

After comments, the last equation was used for the range $0.90 \leq F(a) \leq 0.99$. The table below gives the results. $$\left( \begin{array}{ccc} F(a) & \text{approximation} & \text{exact} \\ 0.90 & 1.28164 & 1.28155\\ 0.91 & 1.34087 & 1.34076\\ 0.92 & 1.40523 & 1.40507\\ 0.93 & 1.47600 & 1.47579\\ 0.94 & 1.55507 & 1.55477\\ 0.95 & 1.64528 & 1.64485\\ 0.96 & 1.75133 & 1.75069\\ 0.97 & 1.88180 & 1.88079\\ 0.98 & 2.05548 & 2.05375\\ 0.99 & 2.32999 & 2.32635 \end{array} \right)$$