Post Reply 
Accurate Normal Distribution for the HP67/97
12-03-2018, 07:02 PM
Post: #21
RE: Accurate Normal Distribution for the HP67/97
(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(h4) 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 t4 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
Find all posts by this user
Quote this message in a reply
12-03-2018, 08:25 PM
Post: #22
RE: Accurate Normal Distribution for the HP67/97
(12-03-2018 12:44 AM)Albert Chan Wrote:  Sum it all, and ignore O(h4) 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 t4
Find all posts by this user
Quote this message in a reply
12-03-2018, 08:45 PM (This post was last modified: 12-03-2018 08:50 PM by Dieter.)
Post: #23
RE: Accurate Normal Distribution for the HP67/97
(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 t4

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
Find all posts by this user
Quote this message in a reply
12-04-2018, 01:38 PM (This post was last modified: 12-04-2018 09:02 PM by Albert Chan.)
Post: #24
RE: Accurate Normal Distribution for the HP67/97
(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 ex, 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
Find all posts by this user
Quote this message in a reply
12-09-2018, 03:28 PM (This post was last modified: Yesterday 04:43 PM by Albert Chan.)
Post: #25
RE: Accurate Normal Distribution for the HP67/97
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
Find all posts by this user
Quote this message in a reply
12-09-2018, 08:17 PM (This post was last modified: 12-10-2018 10:25 AM by Albert Chan.)
Post: #26
RE: Accurate Normal Distribution for the HP67/97
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. Smile

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


Updated: fix typo in constant 1/√(2 Pi)
Find all posts by this user
Quote this message in a reply
12-09-2018, 08:47 PM
Post: #27
RE: Accurate Normal Distribution for the HP67/97
(12-09-2018 08:17 PM)Albert Chan Wrote:  
PHP Code:
0.398942280401327 exp(-0.5 x); 

There is a digit missing in the constant.

Dieter
Find all posts by this user
Quote this message in a reply
12-11-2018, 02:33 PM (This post was last modified: Yesterday 02:52 PM by Albert Chan.)
Post: #28
RE: Accurate Normal Distribution for the HP67/97
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;
    
0.3989422804014327 exp(-0.5*x*x);
    
+= (y*(1./6) + 0.5) * y2;  // ≈ exp(y) - 1
    
return z;


Note: for calculator program, return z * (1. + y) instead
Calculator might underflow the correction term.
Find all posts by this user
Quote this message in a reply
Yesterday, 10:49 PM (This post was last modified: Yesterday 10:50 PM by John Keith.)
Post: #29
RE: Accurate Normal Distribution for the HP67/97
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.
Find all posts by this user
Quote this message in a reply
Today, 12:59 AM (This post was last modified: Today 02:11 AM by Albert Chan.)
Post: #30
RE: Accurate Normal Distribution for the HP67/97
(Yesterday 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.43280 95260 5
Try round z to FIX-4, we get:

|y| < 21.4 * 5e-5 = 0.00107, y4/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-18 ~ 0.000147, y4/24 < 2.0e-17

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, match true value
Find all posts by this user
Quote this message in a reply
Post Reply 




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