Accurate Normal Distribution for the HP67/97

[Dieter]

Normal distribution for the HP67/97

This program for the HP67 and 97 evaluates various functions of the Standard Normal distribution:

• The lower tail cumulative distribution function P(z), i.e. the Normal integral from –∞ to z
• The upper tail cumulative distribution function Q(z). i.e. the Normal integral from z to +∞
• The two-sided cumulative distribution function A(z), i.e. the Normal integral from –z to +z
• The inverse of the one-sided CDF (quantile) z(P)
• The inverse of the two-sided CDF (quantile) z(A)
• The probability distribution function (PDF) Z(z)
  

Unlike the well-known programs e.g. in the HP Statistic Pacs this one uses different algorithms to achieve much better accuracy. even far out in the distribution tails. There are several methods to do so. Since the HP67 and 97 run rather slow and memory is limited, while on the other hand there are always 26 available data registers that cannot be exchanged for more program steps, the chosen approach uses rational approximations which run reasonably fast. Here the nine coefficients are prestored in the data registers. These can be loaded by a data card – or by your preferred HP67/97 emulator.

Method and accuracy

The Normal CDF is evaluated by two different methods.

For 0 ≤ z ≤ 5 the upper-tail integral Q(z) is approximated by a (Near-)Minimax rational approximation:

\(\large Q(z) \approx e^{-\frac{z^2}{2}} \cdot \frac{1+a_1z+a_2z^2+a_3z^3+a_4z^4}{2+b_1z+b_2z^2+b_3z^3+b_4z^4+b_5z^5}\)

Using 10-digit coefficients the values are

\(\begin{array}{ll}
a_1=0,7981006015 & b_1=3,191970353\\
a_2=0,3111510842 & b_2=2,169125520\\
a_3=0,06328636234 & b_3=0,7932255604\\
a_4=0,005716530175 & b_4=0,1587036976\\
& b_5=0,01432712100
\end{array} \)

If evaluated with sufficient precision, the relative error over the given domain is less than 2,6 E–10.
With a few more digits the error can be reduced to approx. 2,41 E–10.

For z > 5 the well known continued fraction expansion is applied, here with 8 terms:

\(\large Q(z) \approx \frac{1}{\sqrt{2 \pi}} \cdot e^{-\frac{z^2}{2}} \cdot \cfrac{1}{z+\cfrac{1}{z+\cfrac{2}{z+\cfrac{3}{z+\dotsb}}}}\)

The expression is calculated from right to left, starting not with 8/z but 8/(z+1,38) which significantly improves the resulting accuracy for smaller z. This way the relative error for z > 5 stays below 1 E-10 – provided the calculation is performed with sufficient precision.

Due to the HP67/97's limitation to 10 significant digits and some numeric pitfalls the actual results on the HP67/97 will be less accurate. Since it is virtually impossible to verify the results over the complete domain I can only say that according to my results usually 9 significant digits (±1 unit) are achieved. See below for two exceptions.

If you find substantially larger errors, please report here.

The algorithm for the inverse (quantile function) first calculates a rough estimate by means of a simple rational approximation with an error of about ±0,002. The error of this first approximation is taylored for the following correction step that provides the final result. This is a very effective third order extrapolation due to Abramovitz & Stegun (1972) p. 954. With sufficient precision this method is good for about 11 significant digits over the calculator's complete working range down to 1 E–99. Again, the actual accuracy on the 67/97 is less and may drop to about 9 digits. But there is an exception: due to digit cancellation results very close to zero carry less significant digits, e.g. the quantile for a probability of 0,50003 is calculated as 7,5199 E–5. In such cases usually the remaining digits are fine, maybe within ±1 unit tolerance. So the result in FIX DSP 9 (0,000075199) should be OK.

A similar limitation applies to the two-sided CDF for arguments very close to zero. Here you should not expect more than what you see in FIX DSP 9 mode (±1 digit).

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.

The program

Here comes the listing.

Code:

001  LBL B
002  CF 2
003  SF 3
004  GTO 0
005  LBL A
006  CF 2
007  CF 3
008  x<0?
009  SF 2
010  LBL 0
011  ABS
012  GSB 9
013  ENTER
014  F3?
015  +
016  ENTER
017  CHS
018  1
019  +
020  F2?
021  X<>Y
022  RTN
023  LBL 9
024  STO 0
025  8
026  STO I
027  5
028  RCL 0
029  x<=y?
030  GTO 1
031  RCL 0
032  RCL 0
033  1
034  .
035  3
036  8
037  +
038  LBL 2
039  RCL I
040  X<>Y
041  /
042  +
043  DSZ
044  GTO 2
045  RCL 0
046  GSB E
047  X<>Y
048  /
049  RTN
050  LBL 1
051  RCL 0
052  RCL 0
053  RCL 9
054  *
055  RCL 8
056  +
057  *
058  RCL 7
059  +
060  *
061  RCL 6
062  +
063  *
064  RCL 5
065  +
066  *
067  2
068  +
069  STO I
070  R↓
071  RCL 4
072  *
073  RCL 3
074  +
075  *
076  RCL 2
077  +
078  *
079  RCL 1
080  +
081  *
082  1
083  +
084  RCL I
085  /
086  RCL 0
087  GSB E
088  R↓
089  RCL I
090  *
091  RTN
092  LBL D
093  ABS
094  CHS
095  1
096  +
097  2
098  /
099  CF 2
100  GTO 0
101  LBL C
102  CF 2
103  ENTER
104  CHS
105  1
106  +
107  x>y?
108  SF 2
109  x>y?
110  X<>Y
111  LBL 0
112  STO D
113  LN
114  ENTER
115  +
116  CHS
117  √x
118  ENTER
119  ENTER
120  ENTER
121  .
122  3
123  6
124  7
125  *
126  2
127  .
128  3
129  5
130  8
131  +
132  X<>Y
133  .
134  0
135  6
136  6
137  5
138  *
139  1
140  .
141  0
142  8
143  5
144  +
145  R↑
146  *
147  1
148  +
149  /
150  -
151  x<0?
152  CLX
153  GSB 9
154  EEX
155  5
156  *
157  RCL E
158  /
159  RCL D
160  EEX
161  5
162  *
163  RCL E
164  /
165  -
166  EEX
167  5
168  /
169  ENTER
170  ENTER
171  RCL 0
172  x²
173  2
174  *
175  1
176  +
177  6
178  /
179  *
180  RCL 0
181  2
182  /
183  +
184  *
185  *
186  +
187  RCL 0
188  +
189  F2?
190  CHS
191  RTN
192  LBL E
193  STO 0
194  INT
195  x²
196  2
197  /
198  CHS
199  e^x
200  RCL 0
201  INT
202  LSTX
203  FRAC
204  *
205  e^x
206  /
207  RCL 0
208  FRAC
209  x²
210  e^x
211  √x
212  /
213  STO I
214  2
215  PI
216  *
217  1/x
218  √x
219  *
220  STO E
221  RTN

The program expects the coefficients of the rational approximation in R1...R9. If it runs on a real (hardware) 67/97 this can be done by preparing a (single track) data card. The values for the constants have already been mentioned. Be sure to enter all ten digits:

Code:

R1 = 7,981006015 E-01
R2 = 3,111510842 E-01
R3 = 6,328636234 E-02
R4 = 5,716530175 E-03
R5 = 3,191970353
R6 = 2,169125520
R7 = 7,932255604 E-01
R8 = 1,587036976 E-01
R9 = 1,432712100 E-02

The coefficients of the simple rational approximation for the quantile estimate are part of the program code. Of course they can just as well be stored in, say, the secondary registers S0...S3 and recalled from there. This will shorten the program and make the quantile calculation a tiny bit faster, but it requires a double-track data card.

Usage

Calculate the cumulative distribution function:
z   [A]
The lower tail CDF P(z) is returned in X, the upper tail CDF Q(z) in Y.

Calculate the symmetric two-sided cumulative distribution function:
z   [B]
The two-sided CDF A(z) is returned in X, the complement 1–A(z) in Y.

Calculate the quantile for a given lower-tail probability p:
p   [C]

Calculate the quantile for a given two-sided symmetric probability p:
p   [D]

Calculate the probability distribution function Z(z):
z   [E]

Some examples

In a soda water factory a machine fills bottles with an average volume of 503 ml.
The content of the bottles varies slightly with a standard deviation of 5 ml.
Determine the probability that a random soda bottle contains less than 490 ml.

First calculate the Standard Normal variable z:

Code:

 490 [ENTER] 503 [-] 5 [÷]     -2,600000000

Now compute the lower tail CDF:

Code:

     [A]                        0,004661188
     [x<>y]                     0,995338812

So only 0,47% of all bottles will contain less than 490 ml while 99,53% exceed this volume.

How much of the production will fall within ±10 ml around the mean volume?
±10 ml equals ±2 standard deviations.

Code:

  2  [B]                        0,954499736

In which interval around the mean will 98% of the production fall?
So we are here looking for the two-sided quantile.

Code:

0,98 [D]                        2,326347874

The tolerance interval is ±2,326 standard deviations.
In absolute milliliters this is...

Code:

  5  [x]                        11,63173937

So 98% of the production is within 503 ± 11,63 ml.

In the above example all digits displayed in FIX DSP 9 mode are exact.
Here are some other results and their accuracy:

Code:

Q(0,01)     =  0,4960106435     (-2 ULP)
Q(1)        =  0,1586552539     (exact)
P(2)        =  0,9772498681     (exact)
P(-3,14)    =  8,447391737 E-4  (+2 ULP)
P(-6,3)     =  1,488228221 E-10 (-1 ULP)
A(2,3)      =  0,9785517800     (exact)
A(0,0001)   =  0,0000797888     (FIX DSP 9 result ...79789 is within 1 digit of 7,978845595 E-5)

z(0,9)      =  1,281551565      (-1 ULP, truncated)
z(0,95)     =  1,644853627      (exact)
z(0,99)     =  2,326347874      (exact)
z(1E-99)    = -21,16517934      (exact)
z(0,500001) =  0,000002506      (four digits of 2,506628275 E-6)

Cave: this does not mean that this accuracy level can be guaranteed. I have not found a case where the result was not within 1 unit in the 9th place, but please do your own tests. As usual, remarks, corrections and error reports are welcome.

Finally, here are two (zipped) files for use with the Panamatic HP67 emulator. The first version implements the program listed above, the second version has the coefficients of the quantile estimate in registers S0...S3 and thus is a bit shorter.

Dieter

  NormDist_67.zip (Size: 1.21 KB / Downloads: 28)