# 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

$$\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.