Back to index

python-biopython  1.60
chi2.py
Go to the documentation of this file.
00001 # Copyright (C) 2011 by Brandon Invergo (b.invergo@gmail.com)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license. Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 #
00006 # This code is adapted (with permission) from the C source code of chi2.c, 
00007 # written by Ziheng Yang and included in the PAML software package:
00008 # http://abacus.gene.ucl.ac.uk/software/paml.html
00009 
00010 from math import sqrt, log, exp
00011 
00012 def cdf_chi2(df, stat):
00013     if df < 1:
00014         raise ValueError, "df must be at least 1"
00015     if stat < 0:
00016         raise ValueError, "The test statistic must be positive"
00017     x = 0.5 * stat
00018     alpha = df / 2.0
00019     prob = 1 - _incomplete_gamma(x, alpha)
00020     return prob
00021                       
00022 def _ln_gamma_function(alpha):
00023     """Compute the log of the gamma function for a given alpha.
00024     
00025     Comments from Z. Yang:
00026     Returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.  
00027     Stirling's formula is used for the central polynomial part of the procedure.
00028     Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
00029     Communications of the Association for Computing Machinery, 9:684
00030     """
00031     if alpha <= 0:
00032         raise ValueError
00033     x = alpha
00034     f = 0
00035     if x < 7:
00036         f = 1
00037         z = x
00038         while z<7:
00039             f *= z 
00040             z += 1
00041         x = z
00042         f = -log(f)
00043     z = 1 / (x * x)
00044     return  f + (x-0.5)*log(x) - x + .918938533204673             \
00045           + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z \
00046                +.083333333333333)/x
00047         
00048 def _incomplete_gamma(x, alpha):
00049     """Compute an incomplete gamma ratio.
00050     
00051     Comments from Z. Yang:
00052     Returns the incomplete gamma ratio I(x,alpha) where x is the upper 
00053            limit of the integration and alpha is the shape parameter.
00054     returns (-1) if in error
00055     ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
00056     (1) series expansion     if alpha>x or x<=1
00057     (2) continued fraction   otherwise
00058     RATNEST FORTRAN by
00059     Bhattacharjee GP (1970) The incomplete gamma integral.  Applied Statistics,
00060     19: 285-287 (AS32)
00061     """
00062     p = alpha
00063     g = _ln_gamma_function(alpha)
00064     accurate = 1e-8
00065     overflow = 1e30
00066     gin = 0
00067     rn = 0
00068     a = 0
00069     b = 0
00070     an = 0
00071     dif = 0
00072     term = 0
00073     if x == 0:
00074        return 0
00075     if x < 0 or p <= 0:
00076         return -1
00077     factor = exp(p*log(x)-x-g)  
00078     if x > 1 and x >= p:
00079         a = 1 - p
00080         b = a + x + 1
00081         term = 0
00082         pn = [1, x, x+1, x*b, None, None]
00083         gin = pn[2] / pn[3]
00084     else:
00085         gin=1
00086         term=1
00087         rn=p
00088         while term > accurate:
00089             rn += 1
00090             term *= x / rn
00091             gin += term
00092         gin *= factor / p
00093         return gin
00094     while True:
00095         a += 1
00096         b += 2
00097         term += 1   
00098         an = a * term
00099         for i in range(2):
00100             pn[i + 4] = b * pn[i + 2] - an * pn[i]
00101         if pn[5] != 0:
00102             rn = pn[4] / pn[5]
00103             dif = abs(gin - rn)
00104             if dif > accurate:
00105                 gin=rn
00106             elif dif <= accurate*rn:
00107                 break
00108         for i in range(4):
00109             pn[i] = pn[i+2]
00110         if abs(pn[4]) < overflow:
00111             continue
00112         for i in range(4):
00113             pn[i] /= overflow
00114     gin = 1 - factor * gin
00115     return gin