|A fast and compact algorithm for the normal quantile|
Message #1 Posted by Dieter on 21 Apr 2011, 4:58 p.m.
In earlier threads in this forum various improvements have been discussed concerning the most efficient way to evaluate the inverse normal distribution (quantile) on the 34s. Since this might be of some interest for other calculators and/or implementations I think it's okay to start a new topic here.
As usual, it's a tradeoff between memory usage, execution speed and numeric accuracy. The current implementation starts with two quite good initial guesses - one value slightly high, the other a bit low. Then the solver is called to determine the exact result.
We can do better. ;-)
I have tried the one or other approach, and finally it turned out that there is a simple method that converges very fast and might even use less memory than the current method. The very sophisticated 34s solver is not required here. There is a much simpler, yet effective way to determine the quantile with a different algorithm:
First we guess a good estimate for the quantile. Only one single value is required, for instance this one:
p > 0,2 p < 0,2
u = sqrt(2*pi)*(0,5-p) u = -2 * ln p
x = u + u^3 / 6 x = sqrt(-2 * ln(p * sqrt((u-1)*2*pi))) + 0,2/u
Yes, the second term is a bit more complex than before but it's worth it.
After this, a dedicated though simple solver is used. The preferred
algorithm here is the Halley method, which uses both the first and the second derivative of a function. In our case, the first derivative simply is the PDF, and the second derivative is a simple function thereof. This leads to a very compact form of the Halley equation:
Assuming x > 0 and p < 0.5, first the well-known Newton-Raphson quotient f / f' is evaluated:
cdf(x) - p
t = ----------
The pdf already has been evaluated within the cdf routine, so its value is known and no additional calculation is required. Now the new and improved approximation is:
x_new = x + ---------
1 - t*x/2
Look, Pauli - no slow logs required here. ;-)
Okay, how does it perform? I tried this method in Excel/VBA with 15
significant digits. After just two (!) iterations the final result was achieved. The method seems to converge at least quadratic, theoretically it's even cubic. So a third iteration probably should return something like 30 valid digits.
Pauli, what do you think? This should be faster than before and maybe it even requires less memory since no equation for the solver has to be set up. And only three (hard-coded) iterations should return a result with far more valid digits than actually required.
(Edit: corrected error in first guess formula)
Edited: 22 Apr 2011, 6:06 a.m. after one or more responses were posted