Accurate Normal Distribution for the HP67/97

[Dieter]

(12-03-2018 12:44 AM)Albert Chan Wrote:  That is a cool formula ! I was expecting something much worse.

Something much worse would not have fit into the limited HP67/97 memory. ;-)

(12-03-2018 12:44 AM)Albert Chan Wrote:  To prove above correction work:
(...)
Sum it all, and ignore O(h⁴) terms, we get:

Ah, thank you very much. I have never tried a proof for that formula myself. But I wonder what that neglected fourth term with t⁴ will look like. Due to the limited hardware precision it will not make much sense, but anyway. I have a faint idea, but you will know better – Albert?

BTW, the tricky part is calculating t in the first place. Especially evaluating Q(x)–p without too much digit cancellation.

Dieter

[Albert Chan]

(12-03-2018 12:44 AM)Albert Chan Wrote:  Sum it all, and ignore O(h⁴) terms, we get:

x correction = h = t + x/2 · t² + (2x²+1)/6 · t³

I am too lazy to do 4th order correction by hand. Using Mathematica:

z = c Exp[-x²/2];
D[z, {x, 3}] / z ==> -x (-3 + x²)

t = h - (x/2) h^2 + (x²-1)/6 h^3 - x*(x²-3)/24 h^4;

mess = Expand[t + (x/2) t^2 + (2x²+1)/6 t^3 + k t^4];

Simplify[(mess /. h^4 -> c /. h -> 0) / c] ==> k - 1/24 x(7 + 6x²)

4th order x correction = t + x/2 · t² + (2x²+1)/6 · t³ + x*(6x²+7)/24 t⁴

[Dieter]

(12-03-2018 08:25 PM)Albert Chan Wrote:  4th order x correction = t + x/2 · t² + (2x²+1)/6 · t³ + x*(6x²+7)/24 t⁴

Thank you very much. I just found an older VBA function in the Excel sheet where the initial approxmation for the Normal quantile was tested:

Code:

x = x + t * t * t * t * (6 * x ^ 3 + 7 * x) / 24 + t * t * t * (2 * x * x + 1) / 6 + t * t * x / 2 + t
        ----------------------------------------

But I can't say how or when I got to that conclusion. ;-)

With four terms even a not too exciting four-digit approximation like Hastings' yields 15 digit double precision accuracy with just one correction step. Probably even more -- but Excel is limited to these 15 digits.

Dieter

[Albert Chan]

(06-26-2016 09:34 PM)Dieter Wrote:  Evaluating the PDF seems trivial, but accuracy may degrade significantly for large arguments of the exponential function. For example e^-1000/7 = 9,076766360 E-63, but the 67/97 returns 9,076765971 E-63. The error is caused by the fact that the fractional part of the argument carries only seven decimals which leads to an accuracy of merely seven significant digits. That's why the PDF is evaluated in a different way that requires three calls of e^x, but achieves better accuracy.

It might be possible to reduce 3 calls of exp(x) to 2
Let B = INT(x²)/2, I = INT(x), F = FRAC(x):

x²/2 = B + (-B + I²/2 + I F + F²/2)

Example:

10.23456789² / 2
= 52 + (-52 + 10²/2 + 10*0.23456789 + 0.23456789²/2)
~ 52 + (-52 + 50 + 2.3456789 + 0.02751104751)
~ 52 + 0.3731899475

exp(-10.23456789² / 2) ~ exp(-52) exp(-0.3731899475) = 1.797267025e-23

[Albert Chan]

We can use Taylor expansion to calculate Z(x+h) with 1 Exp call

Z(x+h) = Z(x) + Z'(x) h + Z''(x)/2 h² + Z'''(x)/6 h³ + ...
Z(x+h) / Z(x) = 1 - x h + (x² - 1)/2 h² - x (x² - 3)/6 h³ + ...

Example: calculate Z(20.3333 333333)

x = z-score rounded to FIX-4, thus exp(-x²/2) can be calculate accurately.
h = z-score - x, thus |h| ≤ 5e-5

-> above example, x=20.3333, h=3.33333e-5
-> Z(x) = exp(-x²/2) / √(2 Pi) = 6.65095 520534 e-91

Correction to get to Z(x+h):

Code:

c1 = -x h =          -677775.98889 e-9
c2 = (c1+h)(c1-h) / 2 = +229.13459 e-9
c3 = c1 (c2 - h²) / 3 =   -0.05152 e-9
c  = c3 + c2 + c1 =  -677546.90582 e-9

Z(x+h) ~ Z(x) + c Z(x) = 6.64644 887122 e-91, matched true value

[Albert Chan]

Above 1 Exp method, if done on the computer, is almost as fast as direct method.
In other words, the extra precision almost without cost.

PHP Code:

double pdf(double z)            // 1 Exp method
{
    double x = (float) z;       // clear low bits
    double y = x - z;           // y = -h
    double y2 = y * y;
    z = 0.3989422804014327 * exp(-0.5 * x * x);
    x *= y;                     // correction 1
    y = 0.5 * (x+y) * (x-y);    // correction 2
    x += y + (1./3) * x*(y-y2); // correction 3
    return z + x * z;
} 


Updated: fix typo in constant 1/√(2 Pi)

[Dieter]

(12-09-2018 08:17 PM)Albert Chan Wrote:  

PHP Code:

z = 0.398942280401327 * exp(-0.5 * x * x); 


There is a digit missing in the constant.

Dieter

[Albert Chan]

Discovered a simpler and faster 1 Exp method.

Instead of correcting Z(x) to Z(x+h) via Taylor expansion of Z(x+h), do this:

Z(x+h) / Z(x) = exp(-½(x+h)² + ½(x²)) = exp(-x h - h²/2) = exp(y)

exp(y) = 1 + y + y²/2 + y³/6 + ...

|y| ~ |x h|, a very small number, above ... can be ignored.

PHP Code:

double pdf(double z)             // 1 Exp method, revised
{
    double x = (float) z;        // clear low bits
    double y = 0.5*(x-z)*(x+z);  // PDF(z) = exp(y) PDF(x)
    double y2 = y*y;
    z = 0.3989422804014327 * exp(-0.5*x*x);
    y += (y*(1./6) + 0.5) * y2;  // ≈ exp(y) - 1
    return z + y * z;
} 


Note: for calculator program, return z * (1. + y) instead
Calculator might underflow the correction term.

[John Keith]

That is a neat program, although I can only begin to understand it. I tried it on the HP 50 and while it is a smaller and faster program than the 2-Exp program in this thread, it is not as accurate if the mean and variance are other than 0 and 1 respectively. For "standard" values of x it is accurate to +/- 1 ULP.

My implementation is straightforward except for the second line,
double x = (float) z;
which turned into -6 RND i.e. round to 6 digits.

I did not experience any underflow, just a less accurate result.

[Albert Chan]

(12-12-2018 10:49 PM)John Keith Wrote:  My implementation is straightforward except for the second line,
double x = (float) z;
which turned into -6 RND i.e. round to 6 digits.

Round to 6 digits does not work. Example: 4.56789²/2 = 10.4328 095260 5
Try round z to FIX-4, we get:

|y| < 21.4 * 5e-5 = 0.00107, y⁴/24 < 5.5e-14

3rd order correction also enough for C double:
From header file float.h: DBL_EPSILON = 2^-52 ~ 2.22e-16

|y| < 38.6 * 2^-19 ~ 0.000074, y⁴/24 < 1.2e-18

Quote:I did not experience any underflow, just a less accurate result.

An example of underflow, try Z(20.5000 00002):

z = Z(20.5) = 2.21198 438021 e-92
y = 0.5*(20.5 - 20.5000 00002)*(20.5 + 20.5000 00002) ~ -4.1e-8
Other correction terms y²/2, y³/6 are too small, and can be ignored.

-> correction = y z ~ -9.069e-100, underflow to 0

That was why I mentioned do z*(1. + y) for calculator program.

-> Z(20.5000 00002) = Z(20.5) * 0.99999 9959 = 2.21198 428952 e-92, matched true value

[Dieter]

(12-13-2018 12:59 AM)Albert Chan Wrote:  

(12-12-2018 10:49 PM)John Keith Wrote:  My implementation is straightforward except for the second line,
double x = (float) z;
which turned into -6 RND i.e. round to 6 digits.

Round to 6 digits does not work.

That's not what John is doing. Look closer: it's not 6 RND but –6 RND.
Negative arguments round to a certain number of significant digits!
It's like SCI 5 RND on a classic RPN device.

Edit: I just noticed that your example also rounds to 6 significant digits. So you are correct here. If the division by 2 also has to be exact, the argument should be rounded to 5 significant digits instead, i.e. –5 RND. I think this is better than FIX 4 as the latter again rounds to 6 significant digits for all x≥10.

BTW the largest x is not 38,6 but 47,918. The working range of the 12-digit HPs is larger than IEEE 754 double precision, it goes down to 1E–499. You will know what this means for the question we're discussing here. ;-)

Dieter

[pier4r]

(12-03-2018 03:53 AM)Valentin Albillo Wrote:  

(12-02-2018 07:44 PM)Albert Chan Wrote:  I do not have access to the Abramovitz & Stegun book.

You can download it for free using this link (470-page PDF, perfectly legal download):

Handbook ...

V.

As valentin says.

Also wiki https://en.wikipedia.org/wiki/Abramowitz_and_Stegun

"Because the Handbook is the work of U.S. federal government employees acting in their official capacity, it is not protected by copyright in the United States. While it could be ordered from the Government Printing Office, it has also been reprinted by commercial publishers, most notably Dover Publications (ISBN 0-486-61272-4), and can be legally viewed on and downloaded from the web. "

Some torrent with files related to calculators have it as well.

[Dieter]

(12-13-2018 10:02 PM)pier4r Wrote:  

(12-03-2018 03:53 AM)Valentin Albillo Wrote:  You can download it for free using this link (470-page PDF, perfectly legal download):

Handbook ...

V.

As valentin says.

Also wiki https://en.wikipedia.org/wiki/Abramowitz_and_Stegun

For the record: the 470 page version is a really nice and clean scan, but compared to the printed book it is a heavily stripped down version that lacks most of the tables. Yes, nowadays most of these are obsolete, but they add to the somewhat nostalgic feelings you may have while you look at the pages of this book. ;-)

That's why I prefer a different PDF version. With 53 MB It's not much larger than the 470-page version, but it offers no less than 1060 pages (!). No, I can't remember from where I got this one.

But I really like all these tables, e.g. table 23.4 with the sum of positive powers Σkⁿ. This is on page 817 which is missing in the stripped down version. Did you know that the sum of the 10th powers of 1...100 equals 9 59924 14243 42419 24250 ?-)

Dieter

[Albert Chan]

Hi, Dieter

I did not know about the ±499 exponent range. That is huge !

FIX-4 was based on ±99 exponent range, we have Z(x) = 0.0 if |x| ≥ 21.4
So, FIX-4 work: -5 RND for 1 ≤ |x| < 10, -6 RND for |x| ≥ 10

With ±499 exponent range, FIX-4 is too much.

Luckily, with revised 1 Exp method, adding correction terms is easy.

• instead of FIX-4, do -5 RND, this fixed |y| < 48*5e-4 = 0.024
• add 2 correction terms: y⁴/24, y⁵/120. Dropped y⁶/720 < 2.7e-13
  

OTTH, with FIVE correction terms, you might as well just call exp(y)

[Albert Chan]

Another 1 Exp method, but correct from direct method: (Revision 2)

Example: Z(z = 20.3333 333333), using Emu48

B = z²/2 = 206.722 222222
D = exp(-B) / √(2 Pi) = 6.64644 886819 e-91, error = -303 ULP

D is very close to true result, only linear correction is enough !

x = z rounded to 5 digits = 20.333
h = z - x = 0.0003 333333

Z(z) = D (B - x²/2 - x h - h²/2 + 1) = D * 1.00000 000046 = 6.64644 887125 e-91, error = +3 ULP

Edit: if correction underflow is not an issue, D + D (B - x²/2 - x h - h²/2) may be more accurate.
We get D + D * 45.5555 555 e-11 = 6.64644 887122 e-91, matched true result.

[Dieter]

(12-13-2018 11:13 PM)Albert Chan Wrote:  I did not know about the ±499 exponent range. That is huge !

The WP34s (in DP mode) and Free42 Decimal support numbers as low as 1E–6143. ;-)
And even less where with each power of ten below one mantissa digit is lost.

Dieter

[rprosperi]

(12-14-2018 07:42 PM)Dieter Wrote:  The WP34s (in DP mode) and Free42 Decimal support numbers as low as 1E–6143. ;-)
And even less where with each power of ten below one mantissa digit is lost.

Though it's hard to imagine any application where the number of significant digits matter when using numbers of this magnitude... [preparing to be educated how this may be naive...]

[John Keith]

(12-13-2018 08:43 PM)Dieter Wrote:  That's not what John is doing. Look closer: it's not 6 RND but –6 RND.
Negative arguments round to a certain number of significant digits!
It's like SCI 5 RND on a classic RPN device.

Edit: I just noticed that your example also rounds to 6 significant digits. So you are correct here. If the division by 2 also has to be exact, the argument should be rounded to 5 significant digits instead, i.e. –5 RND. I think this is better than FIX 4 as the latter again rounds to 6 significant digits for all x≥10.



Dieter

I did try both -5 RND and -7 RND and both made the error much worse. I did not try other values, however. I just saw Albert's latest posts and haven't tried any of the ideas there.

[pier4r]

(12-13-2018 10:53 PM)Dieter Wrote:  Yes, nowadays most of these are obsolete, but they add to the somewhat nostalgic feelings you may have while you look at the pages of this book. ;-)

I don't think they are useless (and I remember my primary school math books had them, then they disappeared). Why do I think that math tables are somewhat useful today?

It is a bit like paper dictionaries. You search a word and you may stumble in others while you search it. Or you can do a bit of random looking.
Of course digital dictionaries can simulate this, but not as good as paper ones yet.

The same with the mathematical tables. Sometimes curiosity arises from what one experiences. A clear table with certain numbers may induce the reader to explore more. With a calculator, yes, one can reproduce all the tables one wants (actually it could be quite a relaxing task) but one doesn't think about a result until one needs it.
Therefore one doesn't have the table of numbers already compiled in front of him to start to notice patterns or ask questions.

Aside from the above, I believe one could find mathematical tables in PDF or, as I said, one can make a little relaxing project to compile them with a calculator (not a computer!). Especially for large numbers it may be challenging.
Actually it won't be much relaxing if one wants to be sure there are no typos or errors.

I remember Jürgen posted a video where a 19C printed the first 10'000 prime numbers (n1). Well one can go and print (maybe in a digital way?) the results of the mathematical tables.

n1: http://www.hpmuseum.org/forum/thread-10243.html

[Dieter]

(12-14-2018 10:25 PM)John Keith Wrote:  I did try both -5 RND and -7 RND and both made the error much worse. I did not try other values, however. I just saw Albert's latest posts and haven't tried any of the ideas there.

Then you have found your best way to round the input. And indeed rounding to 6 significant digits may be the best choice here. The square of such a number is exact, and the division by 2 will still be exact if the mantissa of x² does not exceed 2. This is true if the mantissa of x does not exceed √20, i.e. up to 4,47213. Which accounts for 44,72% of all possible x. The remaining 55,28% are also exact if the final digit of x² is even. Otherwise the value of x²/2 may be 1/2 ULP high. So an error only occurs in only 27,64% of all cases, and even if this happens the result if only half an ULP off. I have not analyzed what this means for the final Z(z) result, but a slight error in the last digit may always remain in such calculations.

After all there is no guarantee that the e^x function is accurate to half an ULP. Take a look at the 15C Advanced Functions Handbook where the possible errors of different functions are explained in detail. Here an error of slightly more than 1 ULP (but less than 3 ULP) is possible.

Please correct me if I'm wrong here. ;-)

Dieter