3

My apologies if this is not the right place to be posting this question, it is related to numerical stability for stats calculations in R.

I am trying to calculate the F-value for very high values of df2 but it looks like it is numerically unstable:

nrange <- 350000:450000
f <- qf(1e-8, 8, nrange, lower.tail=FALSE)
plot(f ~ nrange)

looks like this: graph of F-values

Essentially at around df2=400000 it is no longer accurate. The question is - does anybody know how I might get round this problem? For example, the F distribution can be approximated to two chi-squares (e.g. http://en.wikipedia.org/wiki/F-distribution#Related_distributions_and_properties), and in the documentation for qf it says something about using qchisq for large d2. Actually qchisq does seem accurate at those values but it is not obvious to me how to implement this. For example

qf(0.05, 8, 100, lower.tail=FALSE)

and

(qchisq(0.05, 8, lower.tail=FALSE)/8) / (qchisq(0.05, 100, lower.tail=FALSE)/100)

do not give the same results.

So, the question is how do I get accurate F-values for large df2? Any help would be greatly appreciated. Thanks!

1 Answer 1

2

One useful thing about open source projects is that they are open source

fortune(250)

As Obi-Wan Kenobi may have said in Star Wars: "Use the source, Luke!"
   -- Barry Rowlingson (answering a question on the documentation of some implementation details)
      R-devel (January 2010)

If you look at the source code for qf

https://svn.r-project.org/R/trunk/src/nmath/qf.c

and specifically this bit

/* fudge the extreme DF cases -- qbeta doesn't do this well.
   But we still need to fudge the infinite ones.
 */

if (df1 <= df2 && df2 > 4e5) {
if(!R_FINITE(df1)) /* df1 == df2 == Inf : */
    return 1.;
/* else */
return qchisq(p, df1, lower_tail, log_p) / df1;
}

You will see that they are fudging values above 4e5. (ignoring df2 entirely by assuming the same outcome as when df2 == Inf)

1
  • Haha thanks. I looked at the source and for the time being I will try using the method they use for df < 400000 (calculating from the qbeta), but I am unsure what the consequences will be! Feb 27, 2013 at 8:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.