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