F distribution on the 41C

12162015, 05:02 PM
(This post was last modified: 12162015 05:05 PM by quantalume.)
Post: #1




F distribution on the 41C
The Stat Pac module (0004115002) 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. 

12162015, 06:05 PM
(This post was last modified: 12162015 06:06 PM by Dieter.)
Post: #2




RE: F distribution on the 41C
(12162015 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. (12162015 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. (12162015 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. #\ (12162015 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 

12162015, 07:55 PM
Post: #3




RE: F distribution on the 41C  
12162015, 08:22 PM
(This post was last modified: 12162015 08:24 PM by Dieter.)
Post: #4




RE: F distribution on the 41C
(12162015 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 almostblack 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 halfintegers. ;) Dieter 

12162015, 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 

12162015, 09:08 PM
Post: #6




RE: F distribution on the 41C
(12162015 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. (12162015 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. Yes. Or try a different approach. There are also series representations that work fine for most cases. Dieter 

12172015, 12:52 AM
(This post was last modified: 12192015 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:
I know you want HP41C code, but the above VB code gives you an idea what the HP41C code must do!!! Not a small task! Namir 

12172015, 02:56 PM
Post: #8




RE: F distribution on the 41C
(12162015 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 JeanMarc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at: http://hp41programs.yolasite.com/statistic.php 

12172015, 03:31 PM
Post: #9




RE: F distribution on the 41C
(12172015 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. 

12172015, 03:36 PM
Post: #10




RE: F distribution on the 41C
(12172015 02:56 PM)Ángel Martin Wrote: I's suggest you look at JeanMarc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at: 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 JeanMarc's excellent work in the past. David 

12172015, 07:37 PM
Post: #11




RE: F distribution on the 41C
(12172015 03:36 PM)quantalume Wrote:(12172015 02:56 PM)Ángel Martin Wrote: I's suggest you look at JeanMarc Baillard's program collection, which has a very comprehensive section on Statistical Distributions posted at: JMB also uses the incomplete regularized Beta function for the Fdistribution. 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 halfintegers 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 

12172015, 09:32 PM
(This post was last modified: 12172015 09:42 PM by Dieter.)
Post: #12




RE: F distribution on the 41C
(12172015 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 halfintegers) 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 

12172015, 10:35 PM
Post: #13




RE: F distribution on the 41C
(12172015 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. ;) 

12172015, 10:53 PM
(This post was last modified: 12172015 10:56 PM by Dieter.)
Post: #14




RE: F distribution on the 41C
(12172015 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 

12172015, 11:05 PM
Post: #15




RE: F distribution on the 41C  
12182015, 08:38 PM
(This post was last modified: 12252015 02:58 PM by Dieter.)
Post: #16




RE: F distribution on the 41C
(12162015 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" And here are the instructions: Code: df1 [ENTER] ' nominator degrees of freedom Example: evaluate Q_{F}(2,7) for 10 and 15 d.o.f. Code: 10 [ENTER] 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 halfintegers 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 halfintegers 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. ;) 

12182015, 10:25 PM
(This post was last modified: 12182015 10:41 PM by Dieter.)
Post: #17




RE: F distribution on the 41C
(12172015 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 halfintegers have to be handled: Code: If k = Int(k) Then g = 1 Else g = sqrt(Pi) ;) Dieter 

12182015, 11:26 PM
Post: #18




RE: F distribution on the 41C
(12182015 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 

12192015, 01:39 AM
Post: #19




RE: F distribution on the 41C
(12182015 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 

12192015, 05:41 AM
Post: #20




RE: F distribution on the 41C
(12182015 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)