Threaded Mode | Linear Mode

quantalume · (This post was last modified: 12-16-2015 05:05 PM by quantalume.)

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.

Dieter · (This post was last modified: 12-16-2015 06:06 PM by Dieter.)

(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 I_x(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. I_x 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–I_x). 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

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?

Dieter · (This post was last modified: 12-16-2015 08:24 PM by Dieter.)

(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

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

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

Namir · (This post was last modified: 12-19-2015 02:35 PM by Namir.)

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

Á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

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.

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

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

Dieter · (This post was last modified: 12-17-2015 09:42 PM by Dieter.)

(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 I_x(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

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. ;-)

Dieter · (This post was last modified: 12-17-2015 10:56 PM by Dieter.)

(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

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.

Dieter · (This post was last modified: 12-25-2015 02:58 PM by 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. ;-)

Dieter · (This post was last modified: 12-18-2015 10:41 PM by Dieter.)

(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

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

**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!

Á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