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

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 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

[attachment=3694]