r - 计算大 df2 的 F 值

Question

如果这不是发布此问题的正确位置，我深表歉意，它与 R 中统计计算的数值稳定性有关。

我正在尝试计算非常高的 df2 值的 F 值，但它看起来在数值上不稳定:

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

看起来像这样:

基本上在 df2=400000 左右它不再准确。问题是 - 有人知道我如何解决这个问题吗？例如，F 分布可以近似为两个卡方(例如 http://en.wikipedia.org/wiki/F-distribution#Related_distributions_and_properties)，并且在 的文档中qf 它说了一些关于使用 qchisq 用于大型 d2 的事情。实际上 qchisq 在这些值上看起来确实准确，但对我来说如何实现它并不明显。例如

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

和

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

不要给出相同的结果。

那么，问题是如何获得大 df2 的准确 F 值？任何帮助将不胜感激。谢谢!

Answer 1

开源项目的一个有用之处在于它们是开源的

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)

如果你查看qf的源代码

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

特别是这一点

/* 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;
}

您会看到它们在 4e5 以上的值是捏造的。 (通过假设与 df2 == Inf 时相同的结果完全忽略 df2)