在Python中计算F分布的p值？

Question

在Python中计算F分布的p值？

pythonmathstatisticsprobability

5

假设我有一个F值和相关的自由度，df1和df2。如何使用Python编程计算与这些数字相关联的p值？

注意：我不接受使用scipy或statsmodels的解决方案。

- Everyone_Else

2

你不想用scipy，那numpy和statsmodels呢？也许你可以展示一下你尝试过的东西。 - Jeff

你确定numpy有这个能力吗？我已经搜索过他们的文档，但是numpy的统计能力非常有限。使用statsmodels也不可行；考虑到这一点，我将编辑问题。 - Everyone_Else

虽然很多人愿意帮忙，但是已经说过很多次了，SO并不是一个代码编写服务。但是，当你明确地排除了应该用于此的包时，你至少应该向我们展示你尝试过什么。 - Jeff

Jeff，这不是我可以编写代码的公式。这个任务应该通过某个包（除了上面提到的）基本上只需要一行代码就能完成 - 我只是不知道哪一个包。我已经在其他地方搜索解决方案，但它们都使用scipy。我也阅读了numpy文档并dir'd了包，但没有提到这个功能。 - Everyone_Else

2

如果您已经有了检验统计量和自由度，那么您可以直接使用 scipy.stats.f.sf。您可以查看 scipy.stats 或 statsmodels 的源代码以获取示例。scipy.stats.f.sf 是 scipy.special 中相应函数的包装器。 - Josef

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ulrich Stern · Accepted Answer

对于F分布的累积分布函数(CDF)（即p值），可以使用规则化（不完全）beta函数I(x; a, b)进行计算，例如参见MathWorld。使用此博客中仅使用math的I(x; a, b)代码，得到的p值为：

1 - incompbeta(.5*df1, .5*df2, float(df1)*F/(df1*F+df2))

以下是一些样本值的结果，与scipy.stats.f.sf匹配：

In [57]: F, df1, df2 = 5, 20, 18

In [58]: 1 - incompbeta(.5*df1, .5*df2, float(df1)*F/(df1*F+df2))
Out[58]: 0.0005812207389501722

In [59]: st.f.sf(F, df1, df2)
Out[59]: 0.00058122073922042188

以防博客消失，这里是代码：

import math

def incompbeta(a, b, x):

    ''' incompbeta(a,b,x) evaluates incomplete beta function, here a, b > 0 and 0 <= x <= 1. This function requires contfractbeta(a,b,x, ITMAX = 200) 
    (Code translated from: Numerical Recipes in C.)'''

    if (x == 0):
        return 0;
    elif (x == 1):
        return 1;
    else:
        lbeta = math.lgamma(a+b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1-x)
        if (x < (a+1) / (a+b+2)):
            return math.exp(lbeta) * contfractbeta(a, b, x) / a;
        else:
            return 1 - math.exp(lbeta) * contfractbeta(b, a, 1-x) / b;

def contfractbeta(a,b,x, ITMAX = 200):

    """ contfractbeta() evaluates the continued fraction form of the incomplete Beta function; incompbeta().  
    (Code translated from: Numerical Recipes in C.)"""

    EPS = 3.0e-7
    bm = az = am = 1.0
    qab = a+b
    qap = a+1.0
    qam = a-1.0
    bz = 1.0-qab*x/qap

    for i in range(ITMAX+1):
        em = float(i+1)
        tem = em + em
        d = em*(b-em)*x/((qam+tem)*(a+tem))
        ap = az + d*am
        bp = bz+d*bm
        d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
        app = ap+d*az
        bpp = bp+d*bz
        aold = az
        am = ap/bpp
        bm = bp/bpp
        az = app/bpp
        bz = 1.0
        if (abs(az-aold)<(EPS*abs(az))):
            return az

    print 'a or b too large or given ITMAX too small for computing incomplete beta function.'