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:
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;
Var erf, xu, y : Double;
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
Function iCDFnorm_Winitzki(x: Double) : Double;
Var erfinv, xu, y, lnx : Double;
xu:=2 * x-1;
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!