F distribution on the 41C
12-16-2015, 05:02 PM (This post was last modified: 12-16-2015 05:05 PM by quantalume.)
Post: #1
 quantalume Member Posts: 101 Joined: Dec 2015
F distribution on the 41C
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:

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:

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:

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.
12-16-2015, 06:05 PM (This post was last modified: 12-16-2015 06:06 PM by Dieter.)
Post: #2
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-16-2015, 07:55 PM
Post: #3
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
Dieter, thanks so much for your input. Is this the sort of method you had in mind for Ix?

12-16-2015, 08:22 PM (This post was last modified: 12-16-2015 08:24 PM by Dieter.)
Post: #4
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-16-2015, 09:03 PM
Post: #5
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
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
12-16-2015, 09:08 PM
Post: #6
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-17-2015, 12:52 AM (This post was last modified: 12-19-2015 02:35 PM by Namir.)
Post: #7
 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: F distribution on the 41C

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 Function

I 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
12-17-2015, 02:56 PM
Post: #8
 Ángel Martin Senior Member Posts: 1,091 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-17-2015, 03:31 PM
Post: #9
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
(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.
12-17-2015, 03:36 PM
Post: #10
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
(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
12-17-2015, 07:37 PM
Post: #11
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-17-2015, 09:32 PM (This post was last modified: 12-17-2015 09:42 PM by Dieter.)
Post: #12
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-17-2015, 10:35 PM
Post: #13
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
(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. ;-)
12-17-2015, 10:53 PM (This post was last modified: 12-17-2015 10:56 PM by Dieter.)
Post: #14
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-17-2015, 11:05 PM
Post: #15
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
(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.
12-18-2015, 08:38 PM (This post was last modified: 12-25-2015 02:58 PM by Dieter.)
Post: #16
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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 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. ;-)
12-18-2015, 10:25 PM (This post was last modified: 12-18-2015 10:41 PM by Dieter.)
Post: #17
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: F distribution on the 41C
(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
12-18-2015, 11:26 PM
Post: #18
 quantalume Member Posts: 101 Joined: Dec 2015
RE: F distribution on the 41C
(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. :-)

David
12-19-2015, 01:39 AM
Post: #19
 rprosperi Senior Member Posts: 3,815 Joined: Dec 2013
RE: F distribution on the 41C
(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!

--Bob Prosperi
12-19-2015, 05:41 AM
Post: #20
 Ángel Martin Senior Member Posts: 1,091 Joined: Dec 2013
RE: F distribution on the 41C
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)