Normal Distribution Stuff for 33S (LONG) Message #5 Posted by Les Wright on 2 Mar 2007, 2:31 a.m., in response to message #4 by Les Wright
Here are those listings.
Please note that the original routines for the HP41 are due to J-M Baillard and come from here, though I have necessarily had to retool things to work with the limitations of the 33S.
Routine E computes erf and erfc for any real input using the alternating series expansion (http://mathworld.wolfram.com/Erf.html equation (6)) for x < 1.8 and the continued fraction expansion (a version of equation (33) on the same page, computed out to 35 terms) for larger values.
Routine Q computes the upper tail and cumulative normal distribution using routine E and the relationship 1- P(z) = Q(z) = 0.5 erfc(z/sqrt(2)), where Q(z) is the upper tail and P(z) is the cumulative probability associated with z from a standard normal distribution (mean 0, variance 1).
Routine I is basically a programmed equation to compute the inverse distribution using Solve.
Thanks to the speed of the 33S, E and Q are surprisingly fast, returning a result in no more than a couple of seconds or so.
Computing the inverse probabilities takes longer since it uses the solver and requires repeat evaluations of Q, but with decent initial guesses will return a result in under a minute, I find. Note in routine I I take the difference of logarithms. This is especially helpful when the probabilities get really tiny and it helps the Solver avoid converging prematurely since the difference between two tinies is a tiny and may signal convergence before it is desired. Taking logarithms lets the solver work with numbers of more normal magnitude.
E0001 LBL E
E0002 STO A
E0003 ABS
E0004 1.8
E0005 x<=y?
E0006 GTO C
E0007 CLx
E0008 STO i
E0009 Rv
E0010 STO R
E0011 STO S
E0012 x^2
E0013 STO X
E0014 LASTx
L0001 LBL L
L0002 STO T
L0003 ISG i
L0004 SGN
L0005 2
L0006 RCL* i
L0007 1
L0008 +
L0009 1/x
L0010 ENTER
L0011 +
L0012 1
L0013 -
L0014 RCL/ i
L0015 RCL* X
L0016 STO* R
L0017 RCL R
L0018 STO+ S
L0019 RCL T
L0020 RCL S
L0021 x#y?
L0022 GTO L
L0023 ENTER
L0024 +
L0025 PI
L0026 SQRT
L0027 /
L0028 ENTER
L0029 +/-
L0030 1
L0031 +
L0032 x<>y
L0033 RCL A
L0034 x<0?
L0035 GTO N
L0036 Rv
L0037 RTN
C0001 LBL C
C0002 CLx
C0003 35
C0004 STO i
C0005 Rv
C0006 STO X
C0007 CLx
M0001 LBL M
M0002 RCL+ X
M0003 1/x
M0004 0.5
M0005 RCL* i
M0006 *
M0007 DSE i
M0008 GTO M
M0009 RCL+ X
M0010 1/x
M0011 RCL X
M0012 x^2
M0013 +/-
M0014 e^x
M0015 *
M0016 PI
M0017 SQRT
M0018 /
M0019 ENTER
M0020 +/-
M0021 1
M0022 +
M0023 RCL A
M0024 x<0?
M0025 GTO N
M0026 Rv
M0027 RTN
N0001 LBL N
N0002 Rv
N0003 ENTER
N0004 ENTER
N0005 1
N0006 +
N0007 x<>y
N0008 +/-
N0009 RTN
Q0001 LBL Q
Q0002 2
Q0003 SQRT
Q0004 /
Q0005 XEQ E
Q0006 Rv
Q0007 2
Q0008 /
Q0009 ENTER
Q0010 +/-
Q0011 1
Q0012 +
Q0013 x<>y
Q0014 RTN
I0001 LBL I
I0002 INPUT Q
I0003 INPUT Z
I0004 RCL Z
I0005 XEQ Q
I0006 LOG
I0007 RCL Q
I0008 LOG
I0009 -
I0010 RTN
Given any real x for input, XEQ E will return erf(x) in the X register and erfc(x) in the Y register.
Given any real z for input, XEQ Q will return the upper tail normal probability associated with z in register X, and the complement lower tail cumulative probability in Y.
For example 2 XEQ E returns erf(2) = 9.99532223502e-1 in register X and erfc(2) = 4.677734981e-3 in register Y.
Similarly 2 XEQ Q returns Q(2) = 2.2750131948e-2 in register X (the upper tail probability associated with z = 2), and P(2) = 1 - Q(2) = 9.7724986805e-1 in register Y.
To compute the inverse normal distribution, use the routine I with the Solver.
For example, to compute the z associated with an upper tail probability of 0.0001:
1. Set FN = I.
2. Place guesses for Z on the stack--I just use 3 ENTER 3 here.
3. Execute SOLVE Z.
4 At the prompt Q?, enter 1e-4 R/S.
This takes a few seconds but will return Z = 3.71901648545. Checking it with XEQ Q returns 0.0001 in register X, 0.9999 in register Y.
Hope you get some use out of this.
Les
|