在python中計算F分佈p值？

Question

假設我有一個F值和相關的自由度df1和df2。我如何使用python以編程方式計算與這些數字相關的p值？

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

Answer 1

可以用正則化（不完整）β函數I(x; a, b)來計算F分佈的CDF（以及因此的p值），參見例如MathWorld。從這個blog，只使用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.'