+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: F distribution on the 41C (/thread-5348.html)
F distribution on the 41C - quantalume - 12-16-2015 05:02 PM
The Stat Pac module (00041-15002) has a couple of ANOVA programs that provide the usual outputs including the F ratio. However, you have to then consult the F tables in order to determine the significance level of your result. Who wants to carry around printed tables when you have a computer in your pocket? I've been looking for a good implementation of the F distribution on the 41C, and the only one I could find is from the HP Business Stat User's Library book (can be found on TOS). The author gives this equation for P value:
[attachment=2962]
but actually doesn't use it at all in the implementation. Instead, he uses 26.6.4 and 26.6.5 from Abramowitz and Stegun:
[attachment=2963][attachment=2964]
The problem is that he doesn't bother to check for the case of v1 and v2 both odd, and the program gives wildly inaccurate results. In this case, equation 26.6.8 from A&S applies:
[attachment=2965]
This looks a bit messy for 41C implementation. Does anyone know of a better F distribution program for the 41C before I dive in on this? Thanks.
RE: F distribution on the 41C - Dieter - 12-16-2015 06:05 PM
(12-16-2015 05:02 PM)quantalume Wrote: I've been looking for a good implementation of the F distribution on the 41C, and the only one I could find is from the HP Business Stat User's Library book (can be found on TOS). The author gives this equation for P value:
Yes, that's the upper tail F integral.
(12-16-2015 05:02 PM)quantalume Wrote: but actually doesn't use it at all in the implementation. Instead, he uses 26.6.4 and 26.6.5 from Abramowitz and Stegun:
Of course the program *does* use the P(x) formula, but since an integral has to be computed and there is no Integrate function on the 41, this integral is evaluated by means of the series expansions you mentioned.
(12-16-2015 05:02 PM)quantalume Wrote: The problem is that he doesn't bother to check for the case of v1 and v2 both odd, and the program gives wildly inaccurate results.
The usual mix of poor choice of method and sloppy implementation. #-\
(12-16-2015 05:02 PM)quantalume Wrote: This looks a bit messy for 41C implementation. Does anyone know of a better F distribution program for the 41C before I dive in on this? Thanks.
Simple. Just write a routine that implements the regularized Incomplete Beta function Ix(a, b) and calculate the F integral with j and k d.o.f. by setting x = k/(j·F+k), a=k/2 and b=j/2 – cf. A&S 26.6.2. Ix can be evaluated with a continued fraction approach, using the (modified) Lentz method. For best accuracy use the reflection formula (let x := 1–x, exchange a and b, result is 1–Ix). This method universally works for odd and even d.o.f.
If you add an implementation of the regularized Incomplete Gamma function you can evaluate all major distributions (Normal, Student's t, Chi², Fisher's F, Binomial, Poisson) by means of these two functions. This is the way my distributions pack for the 35s works. You will find the respective thread in the General Software Library, especially note post #54 ff. and the last post which contains a status file for the 35s emulator.
Dieter
RE: F distribution on the 41C - quantalume - 12-16-2015 07:55 PM
Dieter, thanks so much for your input. Is this the sort of method you had in mind for Ix?
[attachment=2966]
RE: F distribution on the 41C - Dieter - 12-16-2015 08:22 PM
(12-16-2015 07:55 PM)quantalume Wrote: Dieter, thanks so much for your input. Is this the sort of method you had in mind for Ix?
I'm not sure since the formula's png file is displayed as black text on almost-black background... #-)
Take a look at the implementation suggested in "Numerical Recipes in C", section 6.4. The respective pages are freely available online in PDF format. However, you still need a Gamma or lnGamma function... or a another way of calculating B(a,b) where a and b may be half-integers. ;-)
Dieter
RE: F distribution on the 41C - quantalume - 12-16-2015 09:03 PM
Hmmm...the PNG renders satisfactorily on a couple of different devices I tried. Perhaps I should use JPG in the future. Anyway, you are right about needing Gamma. I thought I was home free until I remembered we are dividing the d.o.f.s by 2 before substitution. At least B(a,b) only needs to be computed once.
This looks to have everything I need in one place:
http://www.hpmuseum.org/software/41/41gamdgm.htm
RE: F distribution on the 41C - Dieter - 12-16-2015 09:08 PM
(12-16-2015 09:03 PM)quantalume Wrote: Hmmm...the PNG renders satisfactorily on a couple of different devices I tried. Perhaps I should use JPG in the future.
No, PNG is much better in this case. Just avoid transparent layers. After deactivating the Alpha channel the formula displays fine.
(12-16-2015 09:03 PM)quantalume Wrote: Anyway, you are right about needing Gamma. I thought I was home free until I remembered we are dividing the d.o.f.s by 2 before substitution. At least B(a,b) only needs to be computed once.
This looks to have everything I need in one place:
http://www.hpmuseum.org/software/41/41gamdgm.htm
Yes. Or try a different approach. There are also series representations that work fine for most cases.
Dieter
RE: F distribution on the 41C - Namir - 12-17-2015 12:52 AM
This is your lucky day!
Here is a Visual Basic implementation for the cumulative distribution function for the F statistics. It is based on the Handbook of Mathematical Function. I just coded it today for my upcoming Console Programmable RPN Calculation app for Windows. It handles odd and even degrees of freedom. I checked the results against the HP Prime and the function did well. The function uses an implementation of Gamma function, which I also include.
Code:
Function Fpdf(ByVal F As Double, ByVal Df1 As Integer, ByVal Df2 As Integer) As Double
' Use equations on page 946 (sections 26.6.4 and 26/6/5) of "Handbook of Mathematical Functions", 9th printing, 1970,
' by Milton Abramowitz And Irene Stegun
Dim ddf1, ddf2 As Double
Dim Prod, Sum As Double
Dim Theta, X0 As Double
Dim X, A, B As Double
Dim CosTheta, SinTheta As Double
Dim CosPwr, SinPwr As Double
X = Df2 / (Df2 + Df1 * F)
Sum = 0
If Df1 Mod 2 = 0 Then
If Df1 = 2 Then
Sum = 1
ElseIf Df1 = 4 Then
Sum = 1 + Df2 / 2 * (1 - X)
Else
Prod = Df2 * (Df2 + 2) / 2 / 4
Sum = 1 + Df2 / 2 * (1 - X) + Prod * (1 - X) ^ 2
For ddf1 = 8 To Df1 Step 2
Prod *= (Df2 + ddf1 - 4) / (ddf1 - 2)
Sum += Prod * (1 - X) ^ ((ddf1 - 2) \ 2)
Next
End If
F = X ^ (Df2 / 2) * Sum
ElseIf Df2 Mod 2 = 0 Then
If Df2 = 2 Then
Sum = 1
ElseIf Df2 = 4 Then
Sum = 1 + Df1 / 2 * X
Else
Prod = Df1 * (Df1 + 2) / 2 / 4
Sum = 1 + Df1 / 2 * X + Prod * X ^ 2
For ddf2 = 8 To Df2 Step 2
Prod *= (Df1 + ddf2 - 4) / (ddf2 - 2)
Sum += Prod * X ^ ((ddf2 - 2) \ 2)
Next
End If
F = 1 - (1 - X) ^ (Df1 / 2) * Sum
Else ' df1 and df2 are odd
Theta = Math.Atan(Math.Sqrt(Df1 / Df2 * F))
CosTheta = Math.Cos(Theta)
SinTheta = Math.Sin(Theta)
' Calculate A
If Df2 > 1 Then
Sum = CosTheta
Prod = 1
CosPwr = CosTheta
For ddf2 = 5 To Df2 Step 2
Prod *= (ddf2 - 3) / (ddf2 - 2)
CosPwr *= CosTheta * CosTheta
Sum += Prod * CosPwr
Next
A = 2 / Math.PI * (Theta + SinTheta * Sum)
Else
A = 2 * Theta / Math.PI
End If
' Calculate B
If Df2 > 1 And Df1 > 1 Then
Sum = 1
Prod = 1
SinPwr = 1
For ddf1 = 5 To Df1 Step 2
Prod *= (Df2 + ddf1 - 4) / (ddf1 - 2)
SinPwr *= SinTheta * SinTheta
Sum += Prod * SinPwr ' SinTheta ^ (ddf1 - 3)
Next
B = 2 / Math.Sqrt(Math.PI) * Gamma((Df2 - 1) / 2 + 1) / Gamma((Df2 - 2) / 2 + 1) * SinTheta * (CosTheta ^ Df2) * Sum
ElseIf Df1 = 1 Then
B = 0
Else
B = 0
End If
' calculate result
F = 1 - A + B
End If
Return 1 - F
End Function
Function Gamma(ByVal X As Double)
Const MAX = 26
Dim Sum, XPower, Prod As Double
Static C(MAX) As Double, DataFlag As Integer
Dim I As Integer
If DataFlag <> 123 Then
DataFlag = 123
C(1) = 1
C(2) = 0.577215664901533
C(3) = -0.655878071520254
C(4) = -0.0420026350340952
C(5) = 0.166538611382292
C(6) = -0.0421977345555443
C(7) = -0.009621971527877
C(8) = 0.007218943246663
C(9) = -0.0011651675918591
C(10) = -0.0002152416741149
C(11) = 0.0001280502823882
C(12) = -0.0000201348547807
C(13) = -0.00000125049348311
C(14) = 0.00000113302723202
C(15) = -0.00000020563384173
C(16) = 0.000000006116095
C(17) = 0.00000000500200753
C(18) = -0.0000000011812746
C(19) = 0.00000000010434271
C(20) = 0.00000000000778234
C(21) = -0.0000000000036968
C(22) = 0.00000000000051
C(23) = -0.0000000000000206
C(24) = -0.0000000000000054
C(25) = 0.0000000000000014
C(26) = 0.0000000000000001
End If
Prod = 1
Do While X > 1
X -= 1
Prod *= X
Loop
Sum = 0
XPower = X
For I = 1 To MAX
Sum += C(I) * XPower
XPower *= X
Next I
Return Prod / Sum
End FunctionI know you want HP-41C code, but the above VB code gives you an idea what the HP-41C code must do!!! Not a small task!
Namir
RE: F distribution on the 41C - Ángel Martin - 12-17-2015 02:56 PM
(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?
I's suggest you look at Jean-Marc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at:
http://hp41programs.yolasite.com/statistic.php
RE: F distribution on the 41C - quantalume - 12-17-2015 03:31 PM
(12-17-2015 12:52 AM)Namir Wrote: This is your lucky day!
Any day that I can get free, expert advice from the bright minds that inhabit this forum is a great day, my friend!
Namir Wrote:Here is a Visual Basic implementation for the cumulative distribution function for the F statistics.
Thanks so much for that. I'll print it out now so I can study the algorithm.
RE: F distribution on the 41C - quantalume - 12-17-2015 03:36 PM
(12-17-2015 02:56 PM)Ángel Martin Wrote: I's suggest you look at Jean-Marc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at:
http://hp41programs.yolasite.com/statistic.php
Thank you, Ángel. Unfortunately, my office network has flagged that site as containing malware, and I will have to wait until I arrive home later today to visit it. I have seen Jean-Marc's excellent work in the past.
David
RE: F distribution on the 41C - Dieter - 12-17-2015 07:37 PM
(12-17-2015 03:36 PM)quantalume Wrote:(12-17-2015 02:56 PM)Ángel Martin Wrote: I's suggest you look at Jean-Marc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at:
http://hp41programs.yolasite.com/statistic.php
Thank you, Ángel. Unfortunately, my office network has flagged that site as containing malware, and I will have to wait until I arrive home later today to visit it. I have seen Jean-Marc's excellent work in the past.
JMB also uses the incomplete regularized Beta function for the F-distribution. Which indeed is the obvious choice here.
I now got a short program for the HP41 running. It uses one unified algorithm for all cases (odd and even d.o.f. in any combination) via the regularized Beta function. Here the (complete) Beta function is evaluated directly, which can be done without much effort because the underlying Gamma code only has to handle integers and half-integers so that the 41's FACT function can be used. This makes Beta reasonably fast and accurate, only both d.o.f. must not exceed 69.
At the moment the incomplete Beta function uses code from the HP67/97 library, which of course can be optimized. That's why I do not want to post this more or less experimental code now. Just want to say it's not that complicated (about 160 lines now). ;-)
Dieter
RE: F distribution on the 41C - Dieter - 12-17-2015 09:32 PM
(12-17-2015 07:37 PM)Dieter Wrote: At the moment the incomplete Beta function uses code from the HP67/97 library, which of course can be optimized. That's why I do not want to post this more or less experimental code now. Just want to say it's not that complicated (about 160 lines now). ;-)
OK, I now got something that seems to work, based on the regularized Beta algorithm of the HP35s statistical distributions pack. Which in turn is a variant of the modified Lentz method suggested in the Numerical Recipes book.
The whole thing consists of about 200 lines of code. The regularized Beta function Ix(a,b) and Euler's Beta B(a,b) (...for integers and half-integers) can be called separately. So if you like you can also set up Student's t and the Binomial distribution quite easily. Compared to the timings given on JMB's website this version runs about 30% faster.
Maybe I should do some tests now...
Dieter
RE: F distribution on the 41C - quantalume - 12-17-2015 10:35 PM
(12-17-2015 09:32 PM)Dieter Wrote: Compared to the timings given on JMB's website this version runs about 30% faster.
Great! I was able to get to JMB's site using my wireless connection, and I did enter and run his "FDF" program and got the expected results. There is only one thing I don't like about his implementation: the BETA program expects the arguments to B(a,b) entered in the order (b,a) but the IBETA program expects the arguments to Bx(a,b) entered in the order (a,b,x). I know, consistency is the hobgoblin of small minds. ;-)
RE: F distribution on the 41C - Dieter - 12-17-2015 10:53 PM
(12-17-2015 10:35 PM)quantalume Wrote: Great! I was able to get to JMB's site using my wireless connection, and I did enter and run his "FDF" program and got the expected results. There is only one thing I don't like about his implementation: the BETA program expects the arguments to B(a,b) entered in the order (b,a) but the IBETA program expects the arguments to Bx(a,b) entered in the order (a,b,x). I know, consistency is the hobgoblin of small minds. ;-)
Hmmm... if BETA is Euler's Beta function the order of the arguments doesn't matter.
B(a,b) = B(b,a). ;-)
Dieter
RE: F distribution on the 41C - quantalume - 12-17-2015 11:05 PM
(12-17-2015 10:53 PM)Dieter Wrote: Hmmm... if BETA is Euler's Beta function the order of the arguments doesn't matter.
B(a,b) = B(b,a). ;-)
Dieter
Ahh then, it's only an inconsistency in the documentation.
RE: F distribution on the 41C - Dieter - 12-18-2015 08:38 PM
(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 ENDAnd 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-integersExample: evaluate QF(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. ;-)
RE: F distribution on the 41C - Dieter - 12-18-2015 10:25 PM
(12-17-2015 12:52 AM)Namir Wrote: Here is a Visual Basic implementation for the cumulative distribution function for the F statistics. It is based on the Handbook of Mathematical Function.
Namir, thank you very much for this extensive implementation.
But why do you do it this complicated instead of a simple regularized Beta function? Which could be used for other distributions as well.
Could you try x=3 with 1 and 5 degrees of freedom? I tried your code in Excel VBA (with the usual syntax modifications for my older version) and got a result of 0,535... while the correct result is 0,856... The problem occurs for df1=1 and df2 odd. Maybe I did something wrong here?
BTW I like the Gamma implementation which seems very exact. It seems to use a polynomial approximation for x=0...1. For comparison, here's my humble Gamma(k) function which I use for the distributions – after all only integers and half-integers have to be handled:
Code:
If k = Int(k) Then g = 1 Else g = sqrt(Pi)
Do While k > 1.25
k = k - 1
g = g * k
Loop
Gamma = g;-)
Dieter
RE: F distribution on the 41C - quantalume - 12-18-2015 11:26 PM
(12-18-2015 08:38 PM)Dieter Wrote: 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.
Looks very good! Preliminary testing found no errors, and reduction in execution time compared to JMB's solution varies but is up to 50%. I made a few minor changes; I replaced the LBL 70 - 90 with LBL 02 - 05 (single byte). This reduced program length by 7 bytes and improved speed slightly. Your code is 74 bytes shorter than JMB's, although JMB's includes a standalone continued fractions module.
Dieter Wrote:The major limitation of this program is the Beta function which ... means that odd d.o.f. can be up to 70 and even d.o.f. up to 140.
I don't think anyone has any business doing problems with dof over 70 on a pocket calculator, let alone a vintage one. :-)
Thanks for your efforts!
David
RE: F distribution on the 41C - rprosperi - 12-19-2015 01:39 AM
(12-18-2015 11:26 PM)quantalume Wrote: I don't think anyone has any business doing problems with dof over 70 on a pocket calculator, let alone a vintage one. :-)
Ahh, I see you are new here. I'm afraid common sense about what folks want and expect out of pocket calculators does not apply here. Which of course is what makes it so interesting!
RE: F distribution on the 41C - Ángel Martin - 12-19-2015 05:41 AM
(12-18-2015 08:38 PM)Dieter Wrote: 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.
Definitely worth an article in the Forum section... or even an entry in the Program Library