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