[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
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
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
55
56          Returns the incomplete gamma ratio I(x,alpha) where x is the upper
57                 limit of the integration and alpha is the shape parameter.
58          returns (-1) if in error
59          ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
60          (1) series expansion     if alpha>x or x<=1
61          (2) continued fraction   otherwise
62          RATNEST FORTRAN by
63          Bhattacharjee GP (1970) The incomplete gamma integral.  Applied Statistics,
64          19: 285-287 (AS32)
65
66      """
67      p = alpha
68      g = _ln_gamma_function(alpha)
69      accurate = 1e-8
70      overflow = 1e30
71      gin = 0
72      rn = 0
73      a = 0
74      b = 0
75      an = 0
76      dif = 0
77      term = 0
78      if x == 0:
79          return 0
80      if x < 0 or p <= 0:
81          return -1
82      factor = exp(p * log(x) - x - g)
83      if x > 1 and x >= p:
84          a = 1 - p
85          b = a + x + 1
86          term = 0
87          pn = [1, x, x + 1, x * b, None, None]
88          gin = pn[2] / pn[3]
89      else:
90          gin = 1
91          term = 1
92          rn = p
93          while term > accurate:
94              rn += 1
95              term *= x / rn
96              gin += term
97          gin *= factor / p
98          return gin
99      while True:
100          a += 1
101          b += 2
102          term += 1
103          an = a * term
104          for i in range(2):
105              pn[i + 4] = b * pn[i + 2] - an * pn[i]
106          if pn[5] != 0:
107              rn = pn[4] / pn[5]
108              dif = abs(gin - rn)
109              if dif > accurate:
110                  gin = rn
111              elif dif <= accurate * rn:
112                  break
113          for i in range(4):
114              pn[i] = pn[i + 2]
115          if abs(pn[4]) < overflow:
116              continue
117          for i in range(4):
118              pn[i] /= overflow
119      gin = 1 - factor * gin
120      return gin
121
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Mon Jul 10 15:16:48 2017 http://epydoc.sourceforge.net