$\mathcal{N}(0,1)$ c.d.f. and quantile

Question

It is often useful to estimate the quantile for $\Phi^{-1}(d)$ for some real number $d$. Here $\Phi\sim \mathcal{N}(0,1)$. We cannot use for some unknown reasons the Taylor expansion as for estimating say $e^x$, and we have to use a table occurring at the end of some books as an appendix. I would like to be sure that there is no way around this.

EDIT

This doesn't look like $$\Phi^{-1}(d):$$

It's written in Maple. What's wrong here ? Even this is not OK

Answer 1

$$\Phi^{-1}(d) = \sqrt2\,\operatorname{erf}^{-1}(2d - 1)\qquad \text{for} \qquad d\in(0,1)$$

What we can do is to build a $[2n+1,2n]$ Padé approximant $P_n$ around $d=\frac 12$ and obtain for example $$\color{blue}{\Phi^{-1}(d) \sim x \,\, \frac{111547 x^4-1846740 x^2+3643920}{307995 x^4-2454060 x^2+3643920}}\qquad \color{red}{x=\sqrt{2 \pi } \left(d-\frac{1}{2}\right)}$$

Computing the norms, $$\Psi_n(\epsilon)=\int_{\frac \epsilon {100}}^{1-\frac \epsilon {100}} \Big[\Phi^{-1}(d) -P_n\Big]^2\, dd$$

$$\left( \begin{array}{ccccccc} n & \epsilon=0& \epsilon=1 & \epsilon=2 & \epsilon=3 & \epsilon=4 & \epsilon=5 \\ 1 & 1.7\times 10^{-2} & 4.2\times 10^{-3} & 1.8\times 10^{-3} & 8.7\times 10^{-4} & 4.5\times 10^{-4} & 2.4\times 10^{-4} \\ 2 & 5.8\times 10^{-3} & 4.9\times 10^{-4} & 1.2\times 10^{-4} & 3.4\times 10^{-5} & 1.1\times 10^{-5} & 4.0\times 10^{-6} \\ 3 & 2.8\times 10^{-3} & 7.2\times 10^{-5} & 8.7\times 10^{-6} & 1.5\times 10^{-6} & 3.1\times 10^{-7} & 7.2\times 10^{-8} \\ 4 & 1.6\times 10^{-3} & 1.2\times 10^{-5} & 7.1\times 10^{-7} & 7.1\times 10^{-8} & 9.1\times 10^{-9} & 1.4\times 10^{-9} \\ 5 & 1.1\times 10^{-3} & 1.9\times 10^{-6} & 6.0\times 10^{-8} & 3.4\times 10^{-9} & 2.8\times 10^{-10} & 2.8\times 10^{-11} \\ \end{array} \right)$$

Edit

For the generation of $P_n$, $$x=\sqrt{2 \pi } \left(d-\frac{1}{2}\right)\implies \Phi^{-1}(d) = \sqrt2\,\text{erf}^{-1}\left(x \sqrt{\frac{2}{\pi }} \right)$$ P[n_]:=FullSimplify[PadeApproximant[Sqrt[2]*InverseErf[x*Sqrt[2/Pi]],{x, 0, {2*n+1, 2*n}}]]

For the remaining

Q[n_]:=Limit[P[n],x->Sqrt[2*Pi]*(d-1/2)]

Psi[n_]:=NIntegrate[(Sqrt[2]*InverseErf[2d-1]-Q[n])^2,{d,0.05,0.95}]