06-26-2016, 09:34 PM
Normal distribution for the HP67/97
This program for the HP67 and 97 evaluates various functions of the Standard Normal distribution:
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 ex, but achieves better accuracy.
The program
Here comes the listing.
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:
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:
Now compute the lower tail CDF:
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.
In which interval around the mean will 98% of the production fall?
So we are here looking for the two-sided quantile.
The tolerance interval is ±2,326 standard deviations.
In absolute milliliters this is...
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:
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
[attachment=3694]
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)
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 ex, 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
[attachment=3694]