Re: 32E's Normal & Inverse Normal Distribution Message #16 Posted by Dieter on 9 Apr 2012, 5:31 p.m., in response to message #12 by Les Wright
Les, thank you for your work and the error plot. However, it seems you used different formulas for the initial quantile guess than those that actually got implemented. ;-)
After a first approach with the solver that required two guesses (one a bit high, the other a bit low) the whole thing was changed to an algorithm that requires just *one* initual guess. This guess is different for 0<p<0,2 and 0,2<p<0,5. The details can be found here in the original discussion.
In VBA it looks like this:
const root2pi = 2.50662827463100
If q >= 0.2 Then
u = (0.5 - q) * root2pi
x = u + u * u * u / 6
Else
u = -2 * Log(q)
x = Sqr(-2 * Log(q * Sqr(u - 1) * root2pi)) + 0.2 / u
End If
Edit: the actual implementation may slightly differ - the last term for small p may be 1/u2 instead of the original 0,2/u. But that will not change much.
The switch at 0,2 was chosen because at this point the absolute (!) errors match quite closely. They are 0,015 resp. 0,019, in both cases a bit low in absolute terms. You are right in that a switch at p=0,25 would be preferable in terms of the relative (!) error. And 0,22 looks best for both. ;-)
Here's a graph of the absolute error. Sorry for the quality, I had to squeeze it below 125 K.
This initial guess is then improved by a Newton-Halley iteration. When evaluated with sufficient precision, three iterations usually return 35...40-digit accuracy - please see the examples in the thread mentioned above.
BTW, the guess for larger p actually uses the first two terms of the series expansion for the Normal quantile. That's why it gets better and finally exact as p approaches 0,5. ;-)
Dieter
Edited: 9 Apr 2012, 6:29 p.m.
|