# How get a good numerical approximation for exponential quantile function?

I want to calculate the "scaled" quantiles for an exponential distribution and I have a function for the inverse CDF, $$iCDFExp(p,a,b)$$ (e.g. $$a=0, b=1$$). The argument $$p$$ is in my application close to $$1.0$$, like $$0.999\dots$$, and in this situation, I get a numerical overflow.

Actually, the argument $$p$$ is also the result of a calculation, on a normal distribution:
I generate a uniform random variable $$u$$, convert it to a standard normal variable, and scale it up with a factor $$s$$ (e.g. $$s=4$$). This way I get an "upscaled" normal variable, and next, I convert it to uniform by using the CDF of the normal distribution. This way I get $$p$$, which I pass to the exponential inverse CDF:

$$y=iCDFExp(CDFNorm(iCDFNorm(u)*s),a,b)$$

Already for moderate $$u$$, like $$0.99$$, I get quite large normal variable, like $$N=iCDFnorm(u)=2.5$$, and with upscaling even $$N*s=10$$, so the CDFnorm is very close to $$1.0$$.

Is there an option to approximate $$CDFNorm(iCDFNorm(u)*s)$$ for large $$u$$, or even the full term $$iCDFExp(CDFNorm(iCDFNorm(u)*s),a,b)$$?

With $$64$$-bits the direct method starts to fail at $$N=8.3$$, but I need correct values also up to approximately $$15$$.

From simple thinking in terms of PDF, I would expect that y should roughly follow $$s*s$$.

Among the Wiki approximations I think that Winitzki fits best. So I wrote two functions.
Function CDFnorm_Winitzki(x: Double) : Double;
Const a=0.140012;
Var erf, xu, y : Double;
Begin
xu:=x/sqrt(2);
y:=(4/Pi+a * xu * xu)/(1+a * xu * xu);
erf:=signum(xu) * sqrt(1-exp(-xu * xu * y));
CDFnorm_Winitzki:=0.5+0.5 * erf
end;

Function iCDFnorm_Winitzki(x: Double) : Double;
Const a=0.140012;
Var erfinv, xu, y, lnx : Double;
Begin
xu:=2 * x-1;
c:=2/Pi/a;
lnx:=ln(1-xu * xu);
y:=sqrt(sqr(c + lnx/2) - lnx/a) - (c + lnx/2);
erfinv:=signum(xu) * sqrt(y);
iCDFnorm_Winitzki:=erfinv * sqrt(2); end;
But still I get the overflow!

• Is this of any interest ? ncbi.nlm.nih.gov/pmc/articles/PMC6214622 – Claude Leibovici Jun 28 at 7:31
• My current workaround is to calculate the maximum possible scale value (around 2.2), and then use this and 2 further values like s=1 and an intermediate value to make a quadratic fit. This is at least better than clipping. – user32038 Jun 28 at 9:00
• Hi Claude, I found this article, but it is not addressing the tail problematic very well, I need to go much beyond 0.999999. Wikipedia has also nice approximations, but just using them does not help. – user32038 Jun 28 at 9:06

## 1 Answer

Quoting the Wikipedia page

"The quantile function of a distribution is the inverse of the cumulative distribution function. The quantile function of the standard normal distribution is called the probit function, and can be expressed in terms of the inverse error function"

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

If you are not too requiring, may be you could use what I wrote here that is to say $$\text{erf}(x)\approx \sqrt{1-\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 }$$ which can easily be inversed.

This is similar to the approximation proposed by Sergei Winitzki (have a look here

• Hi Claude, among the Wiki approximations I also thought that Winitzki fits best. So I wrote two functions, but still get the overflow. – user32038 Jul 2 at 13:27