F distribution on the 41C

[Dieter]

(12-16-2015 05:02 PM)quantalume Wrote:  Does anyone know of a better F distribution program for the 41C before I dive in on this? Thanks.

For the record, here is my attempt at a F distribution program for the HP41. I hope there are not too many errors – try yourself and see what you get, I did not do any thorough tests.

Code:

 01 *LBL"FCDF"
 02  RDN
 03  STO 01
 04  X<>Y
 05  STO 02
 06  R↑
 07  *
 08  +
 09  RCL 01
 10  X<>Y
 11  /
 12  STO 03
 13  2
 14  ST/ 01
 15  ST/ 02
 16  GTO 70
 17 *LBL"IX"
 18  STO 03
 19  RDN
 20  STO 02
 21  RDN
 22  STO 01
 23 *LBL 70
 24  RCL 01
 25  1
 26  +
 27  RCL 01
 28  RCL 02
 29  +
 30  2
 31  +
 32  /
 33  RCL 03
 34  X<Y?
 35  GTO 80
 36  1
 37  RCL 03
 38  -
 39  STO 03
 40  RCL 01
 41  X<> 02
 42  STO 01
 43  XEQ 80
 44  X<>Y
 45  RTN
 46 *LBL 80
 47  1
 48  STO 04
 49  STO 06
 50  RCL 01
 51  RCL 02
 52  +
 53  RCL 03
 54  *
 55  RCL 01
 56  1
 57  +
 58  /
 59  -
 60  X=0?
 61  XEQ 90
 62  1/x
 63  STO 05
 64 *LBL 10
 65  STO 00
 66  ISG 06
 67  ENTER
 68  RCL 06
 69  2
 70  /
 71  INT
 72  LASTx
 73  X≠Y?
 74  GTO 00
 75  RDN
 76  RCL X
 77  RCL 02
 78  -
 79  *
 80  GTO 01
 81 *LBL 00
 82  RDN
 83  RCL 01
 84  +
 85  RCL X
 86  RCL 02
 87  +
 88  *
 89 *LBL 01
 90  RCL 03
 91  *
 92  RCL 01
 93  RCL 06
 94  +
 95  RCL X
 96  1
 97  -
 98  *
 99  /
100  CHS
101  RCL X
102  RCL 05
103  *
104  1
105  +
106  X=0?
107  XEQ 90
108  1/x
109  STO 05
110  X<>Y
111  RCL 04
112  /
113  1
114  +
115  X=0?
116  XEQ 90
117  STO 04
118  *
119  RCL 00
120  ST* Y
121  X<>Y
122  X≠Y?
123  GTO 10
124  RCL 03
125  RCL 01
126  y^x
127  *
128  1
129  RCL 03
130  -
131  RCL 02
132  y^x
133  *
134  RCL 01
135  /
136  STO 00
137  XEQ 70
138  ST/ 00
139  1
140  RCL 00
141  -
142  RCL 00
143  RTN
144 *LBL 90
145  RDN
146  1 E-90
147  RTN
148 *LBL"BETA"
149  STO 02
150  X<>Y
151  STO 01
152 *LBL 70
153  RCL 01
154  XEQ 88
155  STO 06
156  RCL 02
157  XEQ 88
158  ST* 06
159  RCL 01
160  RCL 02
161  +
162  XEQ 88
163  ST/ 06
164  RCL 06
165  RTN
166 *LBL 88
167  ENTER
168  INT
169  X=Y?
170  GTO 11
171  RCL X
172  ST+ X
173  FACT
174  RCL Y
175  FACT
176  /
177  4
178  RCL Z
179  y^x
180  /
181  Pi
182  SQRT
183  *
184  RTN
185 *LBL 11
186  1
187  -
188  FACT
189  END

And here are the instructions:

Code:

df1 [ENTER]                           ' nominator degrees of freedom
df2 [ENTER]                           ' denominator degrees of freedom
 F            XEQ"FCDF"               ' random variable
                             Q(F)     ' upper tail F integral
    [x<>y]                   P(F)     ' lower tail F integral

 a  [ENTER]
 b  [ENTER]
 x            XEQ"IX"
                             Ix(a,b)  ' regularized Beta function
    [x<>y]                 1-Ix(a,b)  ' and its 1-complement

 a  [ENTER]                           ' a and b must be integers
 b            XEQ"BETA"      B(a,b)   ' or half-integers

Example: evaluate Q_F(2,7) for 10 and 15 d.o.f.

Code:

 10  [ENTER]
 15  [ENTER]
 2,7           XEQ"FCDF"     0,040332   ' Q(2,7 | 10, 15)
     [x<>y]                  0,959668   ' P(2,7 | 10, 15)

The major limitation of this program is the Beta function which on the one hand runs fast and does not require any iterations because the Gamma code (at LBL 88) calculates Gamma for half-integers directly, so the result is exact (± roundoff errors). On the other hand the value of (2x)! has to be calculated, which means that BETA can only handle input for half-integers up to 34,5. For integer arguments the usual limit of 70 applies. This means that odd d.o.f. can be up to 69 and even d.o.f. up to 140.

If this is serious limitation the program may be modified to work with ln(Gamma), coded in a new routine that would have to be added.

As usual, all comments and corrections are welcome.

Dieter

Edit: of course the largest odd integer ≤ 70 is not 70 but 69. ;-)