# Adequate Root Finder To Compute The Quantile Function

The cumulative distribution function of the standard normal distribution $$\Phi(z)=\displaystyle\frac{1}{\sqrt{2\pi}}\int_{-\infty}^z e^{-t^2/2}dt$$ cannot be expressed in terms of elementary functions, thus computing values is subject of numerical approximations (or iterating converging series up to a reasonable limit). Ditto for the inverse of it, quantile function or probit, $$\text{probit}(\Phi(z))=\Phi^{-1}(\Phi(z))=z$$. While there exist approximations for the later it may also be set up as root finding procedure for $$\Phi^{-1}(\Phi(z))-z=0$$.
A look at the plot of the probit function

shows from $$\Phi(z)=0.1$$ to $$0.9$$ roughly a straight line. In this section Newton's method works perfectly. In contrast for the lower and upper ends an alternate root finder is more apt.
"Algorithm 39", Areas Under The Normal Curve, uses as approximation for the tails $$\Phi(z)\approx\displaystyle\frac{1}{\sqrt{2\pi}}\exp\left (-\frac{z^2}{2}\right )\cdot f(z)$$ where $$f$$ is a rational function. (Note, "Algorithm 39" computes actually the Q-function, $$Q(z)=1-\Phi(z)$$, a minor detail me think with no impact on the choice for an adequate root finder.)

Question: Which root finder would be the best for the outer sections of this approximation?

Many standard math libraries (e.g. C99, C++11,Fortran 2008) offer implementations of the complementary error function $$\mathrm{erfc}$$ making it worthwhile to take a detour through that. If we let $$\mathrm{normcdfinv}()$$ represent the inverse of the CDF of the normal distribution, and $$\mathrm{erfcinv}()$$ the inverse of the complementary error function, then $$\mathrm{normcdfinv}(a) = -\sqrt{2} \ \mathrm{erfcinv} (2a)$$.

The inverse of the complementary error function lends itself to root finding by Halley iteration:

$$x_{i+1}=x_{i}-\frac{f(x_{i})}{f'(x_{i})}\left[1-\frac{f(x_{i})}{f'(x_{i})} \cdot\frac{f''(x_{i})}{2f'(x_{i})} \right]^{-1}$$

With $$f(x)=\mathrm{erfc}(x)-y$$, we have $$\frac{f(x_{i})}{f'(x_{i})} = -\frac{1}{2}\sqrt{\pi}\exp(x^2)(\mathrm{erfc}(x_{i}) - y)$$ and $$\frac{f''(x_{i})}{2f'(x_{i})}=-x_{i}$$. Given $$y$$, a simple starting guess would be $$x_{0} = \sqrt{-\ln{y}}$$. The complete setup becomes:

$$c := -\frac{1}{2}\sqrt {\pi}$$
$$f := y \gt 1$$
$$\mathrm{if} \ \ (f) \ \ \mathrm{then} \ \ y := 2 - y$$
$$x_{0} := \sqrt{-\ln{y}}$$

Now iterate for as long as necessary:

$$u_{i} := c \cdot \exp (x_{i}^{2}) \cdot (\mathrm{erfc}(x_{i}) - y)$$
$$d_{i} := \frac{u_{i}}{1 + x_{i} \cdot u_{i}}$$
$$x_{i+1} := x_{i} - d_{i}$$

We need to apply a final fix-up based on the original argument transformation:

$$\mathrm{erfcinv} := \mathrm{if} \ \ (f) \ \ \mathrm{then} \ \ -x_{n} \ \ \mathrm{else} \ \ x_{n}$$

A quick check using IEEE-754 double precision suggests that five iterations are sufficient for fully accurate $$\mathrm{erfcinv}$$ results.

• Hi njuffa -- TY for this fast converging algorithm. Since I'm not a C programmer (and FORTRAN 77 under VM/CMS was several decades ago) I testet it in a CAS. With default accuracy after 3..4 iterations $d_i$ is $0$. Nice :) -- But my query is about the ideal root finder for the outer sections of "Algorithm 39". I'll take a closer look at Halley's method as you suggest, also to Euler-Tschebyschow, before I tick off your reply to answer my question. It may take few days, sorry. Commented yesterday