Everyone_Else - 3 months ago 19x
Python Question

# Calculate F-distribution p values in python?

Suppose that I have an F value and the associated degrees of freedom, df1 and df2. How can I use python to programmatically calculate the p value associated with these numbers?

Note: I would not accept a solution using scipy or statsmodels.

The CDF for the F-distribution (and hence the p-value) can be calculated with the regularized (incomplete) beta function `I(x; a, b)`, see, e.g., MathWorld. Using the code for `I(x; a, b)` from this blog, which uses only `math`, the p-value is

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

Here the result for some sample values, matching `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
``````

Just in case the blog disappears, here the code:

``````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.'
``````
Source (Stackoverflow)