Accurate Bernoulli numbers on the 41C, or "how close can you get"?
03-19-2014, 01:37 PM (This post was last modified: 03-19-2014 06:54 PM by Dieter.)
Post: #18
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-19-2014 06:47 AM)Ángel Martin Wrote:  Hi Dieter, glad to see you find the SandMath worth looking at. Which revision are you using? Are there two appendices 9a and 9b, or only one appendix 9?

The version I had was 4.44.44, now updated to 5.55.55 that indeed has additional information on the Normal quantile function.

Actually I do not use Sandmath. I heard about it quite often here, and when I tried a first look at MCODE I came across the manual. At the moment this is all what I got. ;-)

(03-19-2014 06:47 AM)Ángel Martin Wrote:  I say this because in the latest versions (Revision "M", see above link) there are a couple of functions directly related to the quantile, one (ICPF) using the inverse error function on a general Normal distribution (s,m), and another (QNTL) using Halley's iterative method on a standard Normal (0,1).

Halley's method is nice, but it can be done better and faster. ;-)

I now remember that some time ago I set up an algorithm specially taylored for the HP41 and other 10-digit machines. It first calculates an initial estimate with a very simple (1, 2) rational function, requiring just four constants with few digits, and then one single correction step is sufficient for a result with error less than 0,055 ULP, i.e. 5,5 units in the 12th significant digit. All it needs is a sufficiently accurate CDF (cumulative distribution function).

The following method assumes p < 0,5. The rest is handled by CDF(-x) = 1 - CDF(x).

$$u := \sqrt{-2 ln p}$$

$$x := max( 0 , u - \large \frac{2,374 + 0,39 u}{1 + 1,11 u + 0,07157 u^2} )$$

$$t := \large \frac{CDF(x) - p}{PDF(x)}$$

Now apply a third order (!) correction, cf. Abramowitz & Stegun, 10th ed., p. 954:

$$x := x + t + \frac{x}{2}t^2 + \frac{2x^2+1}{6}t^3$$

If evaluated exactly, the result has an error within approx. +0 and -5,5 units in the 12th significant digit throughout the 41's working range (even down to 1E-121).

Here is an implementation in Visual Basic for Excel:
Code:
Function icdf10(p)   signflag = (p > 0.5)   If signflag Then q = 1 - p Else q = p   u = Sqr(-2 * Log(q))   x = u - (2.374 + 0.39 * u) / ((0.07157 * u + 1.11) * u + 1)   If x < 0 Then x = 0   t = (normcdf(x) - q) / normpdf(x)   x = t * t * t * (2 * x * x + 1) / 6 + t * t * x / 2 + t + x   If signflag Then x = -x   icdf10 = x End Function
Here normcdf() and normpdf() refer to the Normal CDF resp. PDF. Please note that the PDF refers to the upper integral, i.e. from x to infinity. Otherwise two small adjustments are required.

The essential magic is in the combination of a dedicated initial guess and a very effective correction afterwards. If the term at t³ is omitted, the result quite exactly matches that of Halley's method and the error is less than 2 units in the 9th digit.

All this should work for p = 0,2 as well as p = 1E-12 or 1E-99 (that's why $$erf^{-1}(2p-1)$$ is not a good idea). Cases close to the distribution center and and the lower working limit (1E-99) require some care at the calculation of CDF(x) - p since this difference may become numerically zero due to underflow or digit cancellation. But there are ways to handle this.

Dieter

Edit: corrected an error and tweaked one of the constants for even better accuracy.
 « Next Oldest | Next Newest »

 Messages In This Thread Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-09-2014, 07:12 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-10-2014, 12:20 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-10-2014, 08:23 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-12-2014, 04:34 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-14-2014, 10:05 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-14-2014, 09:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:00 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:21 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:04 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 09:32 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-15-2014, 12:26 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-15-2014, 06:12 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-15-2014, 07:17 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 12:56 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 08:47 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-18-2014, 09:24 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-19-2014, 06:47 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-19-2014 01:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-19-2014, 07:32 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-19-2014, 08:31 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-20-2014, 06:46 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-20-2014, 12:44 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-20-2014, 02:21 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-20-2014, 09:07 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-21-2014, 05:54 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-21-2014, 01:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-28-2014, 07:56 PM

User(s) browsing this thread: 1 Guest(s)