[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  __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
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
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Wed Dec 17 16:12:37 2014 http://epydoc.sourceforge.net