Back to index

python-biopython  1.60
Functions
Bio.Phylo.PAML.chi2 Namespace Reference

Functions

def cdf_chi2
def _ln_gamma_function
def _incomplete_gamma

Function Documentation

def Bio.Phylo.PAML.chi2._incomplete_gamma (   x,
  alpha 
) [private]
Compute an incomplete gamma ratio.

Comments from Z. Yang:
Returns the incomplete gamma ratio I(x,alpha) where x is the upper 
       limit of the integration and alpha is the shape parameter.
returns (-1) if in error
ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
(1) series expansion     if alpha>x or x<=1
(2) continued fraction   otherwise
RATNEST FORTRAN by
Bhattacharjee GP (1970) The incomplete gamma integral.  Applied Statistics,
19: 285-287 (AS32)

Definition at line 48 of file chi2.py.

00048 
00049 def _incomplete_gamma(x, alpha):
00050     """Compute an incomplete gamma ratio.
00051     
00052     Comments from Z. Yang:
00053     Returns the incomplete gamma ratio I(x,alpha) where x is the upper 
00054            limit of the integration and alpha is the shape parameter.
00055     returns (-1) if in error
00056     ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
00057     (1) series expansion     if alpha>x or x<=1
00058     (2) continued fraction   otherwise
00059     RATNEST FORTRAN by
00060     Bhattacharjee GP (1970) The incomplete gamma integral.  Applied Statistics,
00061     19: 285-287 (AS32)
00062     """
00063     p = alpha
00064     g = _ln_gamma_function(alpha)
00065     accurate = 1e-8
00066     overflow = 1e30
00067     gin = 0
00068     rn = 0
00069     a = 0
00070     b = 0
00071     an = 0
00072     dif = 0
00073     term = 0
00074     if x == 0:
00075        return 0
00076     if x < 0 or p <= 0:
00077         return -1
00078     factor = exp(p*log(x)-x-g)  
00079     if x > 1 and x >= p:
00080         a = 1 - p
00081         b = a + x + 1
00082         term = 0
00083         pn = [1, x, x+1, x*b, None, None]
00084         gin = pn[2] / pn[3]
00085     else:
00086         gin=1
00087         term=1
00088         rn=p
00089         while term > accurate:
00090             rn += 1
00091             term *= x / rn
00092             gin += term
00093         gin *= factor / p
00094         return gin
00095     while True:
00096         a += 1
00097         b += 2
00098         term += 1   
00099         an = a * term
00100         for i in range(2):
00101             pn[i + 4] = b * pn[i + 2] - an * pn[i]
00102         if pn[5] != 0:
00103             rn = pn[4] / pn[5]
00104             dif = abs(gin - rn)
00105             if dif > accurate:
00106                 gin=rn
00107             elif dif <= accurate*rn:
00108                 break
00109         for i in range(4):
00110             pn[i] = pn[i+2]
00111         if abs(pn[4]) < overflow:
00112             continue
00113         for i in range(4):
00114             pn[i] /= overflow
00115     gin = 1 - factor * gin
00116     return gin

Here is the call graph for this function:

Here is the caller graph for this function:

def Bio.Phylo.PAML.chi2._ln_gamma_function (   alpha) [private]
Compute the log of the gamma function for a given alpha.

Comments from Z. Yang:
Returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.  
Stirling's formula is used for the central polynomial part of the procedure.
Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
Communications of the Association for Computing Machinery, 9:684

Definition at line 22 of file chi2.py.

00022 
00023 def _ln_gamma_function(alpha):
00024     """Compute the log of the gamma function for a given alpha.
00025     
00026     Comments from Z. Yang:
00027     Returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.  
00028     Stirling's formula is used for the central polynomial part of the procedure.
00029     Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
00030     Communications of the Association for Computing Machinery, 9:684
00031     """
00032     if alpha <= 0:
00033         raise ValueError
00034     x = alpha
00035     f = 0
00036     if x < 7:
00037         f = 1
00038         z = x
00039         while z<7:
00040             f *= z 
00041             z += 1
00042         x = z
00043         f = -log(f)
00044     z = 1 / (x * x)
00045     return  f + (x-0.5)*log(x) - x + .918938533204673             \
00046           + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z \
00047                +.083333333333333)/x
        

Here is the caller graph for this function:

def Bio.Phylo.PAML.chi2.cdf_chi2 (   df,
  stat 
)

Definition at line 12 of file chi2.py.

00012 
00013 def cdf_chi2(df, stat):
00014     if df < 1:
00015         raise ValueError, "df must be at least 1"
00016     if stat < 0:
00017         raise ValueError, "The test statistic must be positive"
00018     x = 0.5 * stat
00019     alpha = df / 2.0
00020     prob = 1 - _incomplete_gamma(x, alpha)
00021     return prob
                      

Here is the call graph for this function: