(12C Platinum) Normal Distribution

[Gamo]

This Normal Distribution program is adapted from TI SR-56 Applications Library.

Instruction is attached here and procedure is mostly the same.

   

Above Example: FIX 6

2.02 [R/S] display 0.051863574 // Z(x)
[R/S] display 0.021691583 // Q(x)

Remark: Answer is a bit less accurate due to the Pi precisions

Program: Normal Distribution (ALG Mode) 94 Program Lines

Code:


STO 1
X^2
÷
2
=
CHS
e^x
÷
(355 ÷ 113 x 2)
√x
=
STO 2
R/S
-----------------------------
.231642
x
RCL 1
X^2
√x
+
1
=
1/x
STO 3
x
1.33027
-
1.82126
=
x
RCL 3
+
1.78148
=
x
RCL 3
-
.356564
=
x
RCL 3
+
.319382
=
x
RCL 2
x
RCL 3
=

Check answer with your calculator from this Normal Distribution Calculator website.

https://keisan.casio.com/exec/system/1180573188

Percentile X = your desire input
Mean u = 0
Standard Deviation q = 1




Gamo

[Dieter]

(12-01-2018 07:53 AM)Gamo Wrote:  Above Example: FIX 9

See below. ;-)

(12-01-2018 07:53 AM)Gamo Wrote:  Remark: Answer is a bit less accurate due to the Pi precisions

Yes, but why don't you simply use the exact value? Gamo, you are wasting no less than 11 steps to get an approximate √(2pi) value with only seven (!) valid digits:

Code:

÷
(355 ÷ 113 x 2)
√x
=

Why don't you simply use the exact value directly (even with less steps):

Code:

x
,3989422804
=

BTW, on a ten digit calculator multiplying by 1/√2pi is slightly more accurate than dividing by √(2pi) because the value rounds much better to ten digits.

Edit: But let's not bee too picky here. For large x the accuracy for Z(x) will lose up to three digits anyway. I also wonder if the 12C (Platinum) will accept the final digit 4 in the constant. Considering its internal 12-digit precision this might be the best way:

Code:

÷
6,283185307
√x
=

This should give 11 digits of √(2pi): 2,50662827460 instead of ...63.
Requires the same 11 steps as before, but has much better accuracy.

(12-01-2018 07:53 AM)Gamo Wrote:  Program: Normal Distribution (ALG Mode) 94 Program Lines

First of all: while Z(x) works for positive and negative x, the program calculates Q(x) only for x≥0. It does not handle negative x, for instance the program will return the same value for x=1 and –1. Here the user has to apply the formula Q(–x) = 1–Q(x) manually. Of course you can also make the program do so.

The cumulative distribution function Q(x) – the upper tail Normal integral – is approximated by a polynomial, the well-known Hastings approximation which for instance is listed in Abramovitz & Stegun, formula 26.2.17. But... the constants are rounded to six digits (probably due to the SR56's limited memory), so the accuracy of the original formula – a largest absolute error of 7,5 E–8 – is not met. That's why FIX 9 makes no sense here. For x≤2 the relative error is within 5 E–6, so you get about 5 valid significant digits. For larger x, similar to the original approximation, the number of valid digits decreases. For x=5 only three significant digits are left, and finally only two or one are correct.

The accuracy can be improved if the original untruncated coefficients are used:

  p = 0,2316419
b1 = 0,31938153
b2 = -0,356563782
b3 = 1,781477937
b4 = -1,821255978
b5 = 1,330274429

This way you can set FIX 7, and what you see is correct within ±1 ULP. If the 12C Platinum does it like the original 12C, it will eventually switch to scientific notation. At this point you know you shouldn't trust more than two and finally only one significant digit: Try x=8, the exact Q-value here is 6,221 E–16 while the Hastings approximation with the above coefficients returns 6,285 E–16.

That's why I'd prefer a different way of calculating the Normal integral. For instance with a customized rational function. Same effort, better accuracy. ;-)

Here is such an approximation (for x≥0) that I set up some time ago:

Q(x) ≈ e^–x²/2 · (1 + a1·x + a2·x² + a3·x³) / (2 + b1·x + b2·x² + b3·x³)

where
a1 = 0,433563
a2 = 0,084144052
a3 = -9,6367 E-5
b1 = 2,46295834
b2 = 1,13288719
b3 = 0,20590615

Again there are six constants, but with less digits. For x≤5 the largest relative error is less than 1,2 E–6, and the absolute error is < 1 unit in the 6th significant digit. For x>5 this approach is still better than Hastings'. For instance it returns Q(8) as 6,220 E–16. Near the end of the working range, at Q(21) = 3,28 E–98 the above method returns 3,26 E–98 while the Hastings approximation yields 3,56 E–98. All in all even at the underflow limit the relative error is less than 0,007.

But for large x you can apply a continued fraction method that keeps the error on the 1 E–6 level also for x>5. Actually the accuracy increases as x gets larger. I have also designed such a CF method that complements the above rational approximation and requires just three terms.

But I think I'm getting lost here... ;-)

Dieter

[Albert Chan]

(12-01-2018 10:59 AM)Dieter Wrote:  Here is such an approximation (for x≥0) that I set up some time ago:

Q(x) ≈ e^–x²/2 · (1 + a1·x + a2·x² + a3·x³) / (2 + b1·x + b2·x² + b3·x³)

where
a1 = 0,433563
a2 = 0,084144052
a3 = -9,6367 E-5
b1 = 2,46295834
b2 = 1,13288719
b3 = 0,20590615

Using x = 0 to 5 step 0.01, confirmed your claim of 1.2e-6 max rel. error ... Nice!

Hasting approximation (even with best constants) is 1000X worse, max rel. error = 1.6e-3

[Dieter]

(12-01-2018 05:37 PM)Albert Chan Wrote:  Using x = 0 to 5 step 0.01, confirmed your claim of 1.2e-6 max rel. error ...

With exact coefficents it's even a tad less, about 1,10 E–6.
I like that it's "exact by design" for x=0. ;-)

(12-01-2018 05:37 PM)Albert Chan Wrote:  Nice!

It gets even nicer: take a look at this thread in the HP67 software forum. This implements a (4;5) rational approximation.

Finally, here is the CF expansion that complements the above rational approximation for x>5:

Q(x) ≈ Z(x) / (x+1/(x+2/(x+3/(x+0,66))))

The two largest errors are at x=5 (–1,04 E–6) and x=6,13 (+8,9 E–7). For larger x the error decreases continuously towards zero. At the underflow limit (x=21,16517934) it is merely 1,4 E–9. If exactly evaluated, that is. Which is not possible with the standard 12C implementation of Z(x). The linked HP67/97 program evaluates the PDF differently to achieve better accuracy.

Dieter

[Albert Chan]

(12-01-2018 07:44 PM)Dieter Wrote:  

(12-01-2018 05:37 PM)Albert Chan Wrote:  Using x = 0 to 5 step 0.01, confirmed your claim of 1.2e-6 max rel. error ...

With exact coefficents it's even a tad less, about 1,10 E–6.
I like that it's "exact by design" for x=0. ;-)

I see it now.
Q(x = 0) = 1/2 * exp(0) = 1/2 (exactly)

Can you post the exact coefficients too ?

Edit: playing with the coefficients, I managed to lower max rel. errors to 1.12e-6

a1 = 0.433563
a2 = 0.084144052
a3 = -9.63659e-5
b1 = 2.4629583
b2 = 1.1328872
b3 = 0.2059062

[Dieter]

(12-01-2018 11:16 PM)Albert Chan Wrote:  Can you post the exact coefficients too ?

According to my results the following set has a largest error of 1,0064 E-6:

a1 = 0,4335629429
a2 = 0,0841440309
a3 = -9,636444 E-5
b1 = 2,462958275
b2 = 1,132887108
b4 = 0,20590614476

(12-01-2018 11:16 PM)Albert Chan Wrote:  Edit: playing with the coefficients, I managed to lower max rel. errors to 1.12e-6

Fine, even with less digits. But it can still be tweaked a bit:

a1 = 0,433563
a2 = 0,084144052
a3 = -9,6366 E-5
b1 = 2,4629583
b2 = 1,1328872
b3 = 0,2059062

This way the error seems to stay a tiny bit below 1,12E-6. And for calculators every digit counts. For such applications I would now recommend this set of coefficients.

Maybe someone with access to Mathematica, Maple or similar software can post the exact largest errors.

Dieter

[Dieter]

(12-02-2018 09:11 AM)Dieter Wrote:  Fine, even with less digits. But it can still be tweaked a bit:

a1 = 0,433563
a2 = 0,084144052
a3 = -9,6366 E-5
b1 = 2,4629583
b2 = 1,1328872
b3 = 0,2059062

This way the error seems to stay a tiny bit below 1,12E-6.

There now is a new approximation with focus on the absolute error. But in a somewhat different way: With the following coefficients the error is less than 6,5 units in the 7th significant digit.

a1 = 0,4367117
a2 = 0,0851264
a3 = -9,34575 E-5
b1 = 2,46927
b2 = 1,1397766
b3 = 0,2085059

The valid domain is even a bit wider than before: the above error threshold is maintained up to x=5,1993..., which is the point where the CDF reaches 1 E–7. Which is easy to memorize: any result down to 1 E–7 has an error less than 7 units in its 7th significant digit. ;-)

This rational approximation can be complemented with the above CF expansion. Replace the constant 0,66 with 0,6444 and switch between both methods at x=5. This way the error stays below 7 units in the 7th significant digit for all x≥0.

Dieter

[Albert Chan]

(12-02-2018 06:59 PM)Dieter Wrote:  There now is a new approximation with focus on the absolute error. But in a somewhat different way: With the following coefficients the error is less than 6,5 units in the 7th significant digit ...

For x = 0 to 5, confirmed on the 6.5e-7 maximum absolute error ...

However, the old constants did even better, max abs. error = 5.1e-7

[Dieter]

(12-02-2018 08:14 PM)Albert Chan Wrote:  

(12-02-2018 06:59 PM)Dieter Wrote:  There now is a new approximation with focus on the absolute error. But in a somewhat different way: With the following coefficients the error is less than 6,5 units in the 7th significant digit ...

For x = 0 to 5, confirmed on the 6.5e-7 maximum absolute error ...

However, the old constants did even better, max abs. error = 5.1e-7

It's 6,5 units in the 7th significant (!) digit. Not an absolute error of 6,5 E-7. The previous approximation had a largest error of about 1 unit in the 6th significant digit, now this is down to 65% of that.

Compare the results for x=3,1:

Result for the coefficients in post #6: 9,676022 E-04
Compare this to the true result 9,676032 E-04
The error is about -1 unit in the 6th significant digit.

Result for the coefficients in post #7: 9,676026 E-04
Compare this to the true result 9,676032 E-04
The error is about -6 units in the 7th significant digit.

The coefficients in post #7 limit this kind of error to less than 6,5 units in the 7th digit. These are simply two different targets both approximations are designed for.

Dieter