Package Bio :: Package Phylo :: Package PAML :: Module chi2
[hide private]
[frames] | no frames]

Source Code for Module Bio.Phylo.PAML.chi2

  1  # Copyright (C) 2011 by Brandon Invergo (b.invergo@gmail.com) 
  2  # This code is part of the Biopython distribution and governed by its 
  3  # license. Please see the LICENSE file that should have been included 
  4  # as part of this package. 
  5  # 
  6  # This code is adapted (with permission) from the C source code of chi2.c, 
  7  # written by Ziheng Yang and included in the PAML software package: 
  8  # http://abacus.gene.ucl.ac.uk/software/paml.html 
  9   
 10  from math import log, exp 
 11   
 12  __docformat__ = "restructuredtext en" 
 13   
14 -def cdf_chi2(df, stat):
15 if df < 1: 16 raise ValueError("df must be at least 1") 17 if stat < 0: 18 raise ValueError("The test statistic must be positive") 19 x = 0.5 * stat 20 alpha = df / 2.0 21 prob = 1 - _incomplete_gamma(x, alpha) 22 return prob
23 24
25 -def _ln_gamma_function(alpha):
26 """Compute the log of the gamma function for a given alpha. 27 28 Comments from Z. Yang: 29 Returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places. 30 Stirling's formula is used for the central polynomial part of the procedure. 31 Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function. 32 Communications of the Association for Computing Machinery, 9:684 33 """ 34 if alpha <= 0: 35 raise ValueError 36 x = alpha 37 f = 0 38 if x < 7: 39 f = 1 40 z = x 41 while z<7: 42 f *= z 43 z += 1 44 x = z 45 f = -log(f) 46 z = 1 / (x * x) 47 return f + (x-0.5)*log(x) - x + .918938533204673 \ 48 + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z 49 +.083333333333333)/x
50 51
52 -def _incomplete_gamma(x, alpha):
53 """Compute an incomplete gamma ratio. 54 55 Comments from Z. Yang:: 56 57 Returns the incomplete gamma ratio I(x,alpha) where x is the upper 58 limit of the integration and alpha is the shape parameter. 59 returns (-1) if in error 60 ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant. 61 (1) series expansion if alpha>x or x<=1 62 (2) continued fraction otherwise 63 RATNEST FORTRAN by 64 Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics, 65 19: 285-287 (AS32) 66 67 """ 68 p = alpha 69 g = _ln_gamma_function(alpha) 70 accurate = 1e-8 71 overflow = 1e30 72 gin = 0 73 rn = 0 74 a = 0 75 b = 0 76 an = 0 77 dif = 0 78 term = 0 79 if x == 0: 80 return 0 81 if x < 0 or p <= 0: 82 return -1 83 factor = exp(p*log(x)-x-g) 84 if x > 1 and x >= p: 85 a = 1 - p 86 b = a + x + 1 87 term = 0 88 pn = [1, x, x+1, x*b, None, None] 89 gin = pn[2] / pn[3] 90 else: 91 gin=1 92 term=1 93 rn=p 94 while term > accurate: 95 rn += 1 96 term *= x / rn 97 gin += term 98 gin *= factor / p 99 return gin 100 while True: 101 a += 1 102 b += 2 103 term += 1 104 an = a * term 105 for i in range(2): 106 pn[i + 4] = b * pn[i + 2] - an * pn[i] 107 if pn[5] != 0: 108 rn = pn[4] / pn[5] 109 dif = abs(gin - rn) 110 if dif > accurate: 111 gin=rn 112 elif dif <= accurate*rn: 113 break 114 for i in range(4): 115 pn[i] = pn[i+2] 116 if abs(pn[4]) < overflow: 117 continue 118 for i in range(4): 119 pn[i] /= overflow 120 gin = 1 - factor * gin 121 return gin
122