Post Reply 
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
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.
Find all posts by this user
Quote this message in a reply
12-16-2015, 06:05 PM (This post was last modified: 12-16-2015 06:06 PM by Dieter.)
Post: #2
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
Find all posts by this user
Quote this message in a reply
12-16-2015, 07:55 PM
Post: #3
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?

   
Find all posts by this user
Quote this message in a reply
12-16-2015, 08:22 PM (This post was last modified: 12-16-2015 08:24 PM by Dieter.)
Post: #4
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
Find all posts by this user
Quote this message in a reply
12-16-2015, 09:03 PM
Post: #5
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
Find all posts by this user
Quote this message in a reply
12-16-2015, 09:08 PM
Post: #6
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 12:52 AM (This post was last modified: 12-19-2015 02:35 PM by Namir.)
Post: #7
RE: F distribution on the 41C
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 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
Find all posts by this user
Quote this message in a reply
12-17-2015, 02:56 PM
Post: #8
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 03:31 PM
Post: #9
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.
Find all posts by this user
Quote this message in a reply
12-17-2015, 03:36 PM
Post: #10
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 07:37 PM
Post: #11
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 09:32 PM (This post was last modified: 12-17-2015 09:42 PM by Dieter.)
Post: #12
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 10:35 PM
Post: #13
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. ;-)
Find all posts by this user
Quote this message in a reply
12-17-2015, 10:53 PM (This post was last modified: 12-17-2015 10:56 PM by Dieter.)
Post: #14
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
Find all posts by this user
Quote this message in a reply
12-17-2015, 11:05 PM
Post: #15
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.
Find all posts by this user
Quote this message in a reply
12-18-2015, 08:38 PM (This post was last modified: 12-25-2015 02:58 PM by Dieter.)
Post: #16
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. ;-)
Find all posts by this user
Quote this message in a reply
12-18-2015, 10:25 PM (This post was last modified: 12-18-2015 10:41 PM by Dieter.)
Post: #17
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
Find all posts by this user
Quote this message in a reply
12-18-2015, 11:26 PM
Post: #18
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. :-)

Thanks for your efforts!

David
Find all posts by this user
Quote this message in a reply
12-19-2015, 01:39 AM
Post: #19
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
Find all posts by this user
Quote this message in a reply
12-19-2015, 05:41 AM
Post: #20
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
Find all posts by this user
Quote this message in a reply
Post Reply 




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