[frames] | no frames]

# Source Code for Module Bio.CodonAlign.chisq

```  1  """Python implementation of chisqprob, to avoid SciPy dependency.
2
4  """
5
6  import math
7  import sys
8
9  # Cephes Math Library Release 2.0:  April, 1987
10  # Copyright 1985, 1987 by Stephen L. Moshier
11  # Direct inquiries to 30 Frost Street, Cambridge, MA 02140
12  MACHEP = 0.0000001     # the machine roundoff error / tolerance
13  BIG = 4.503599627370496e15
14  BIGINV =  2.22044604925031308085e-16
15
16
17 -def chisqprob(x, df):
18      """
19      Probability value (1-tail) for the Chi^2 probability distribution.
20
22
23      Parameters
24      ----------
25      x : array_like or float > 0
26
27      df : array_like or float, probably int >= 1
28
29      Returns
30      -------
31      chisqprob : ndarray
32          The area from `chisq` to infinity under the Chi^2 probability
33          distribution with degrees of freedom `df`.
34
35      """
36      if x <= 0:
37          return 1.0
38      if x == 0:
39          return 0.0
40      if df <= 0:
41          raise ValueError("Domain error.")
42      if x < 1.0 or x < df:
43          return 1.0 - _igam(0.5*df, 0.5*x)
44      return _igamc(0.5*df, 0.5*x)
45
46
47 -def _igamc(a, x):
48      """Complemented incomplete Gamma integral.
49
50      SYNOPSIS:
51
52      double a, x, y, igamc();
53
54      y = igamc( a, x );
55
56      DESCRIPTION:
57
58      The function is defined by::
59
60          igamc(a,x)   =   1 - igam(a,x)
61
62                                  inf.
63                                     -
64                            1       | |  -t  a-1
65                      =   -----     |   e   t   dt.
66                           -      | |
67                          | (a)    -
68                                      x
69
70      In this implementation both arguments must be positive.
71      The integral is evaluated by either a power series or
72      continued fraction expansion, depending on the relative
73      values of a and x.
74      """
75      # Compute  x**a * exp(-x) / Gamma(a)
76      ax = math.exp(a * math.log(x) - x - math.lgamma(a))
77
78      # Continued fraction
79      y = 1.0 - a
80      z = x + y + 1.0
81      c = 0.0
82      pkm2 = 1.0
83      qkm2 = x
84      pkm1 = x + 1.0
85      qkm1 = z * x
86      ans = pkm1 / qkm1
87      while True:
88          c += 1.0
89          y += 1.0
90          z += 2.0
91          yc = y * c
92          pk = pkm1 * z - pkm2 * yc
93          qk = qkm1 * z - qkm2 * yc
94          if qk != 0:
95              r = pk/qk
96              t = abs((ans - r) / r)
97              ans = r
98          else:
99              t = 1.0;
100          pkm2 = pkm1
101          pkm1 = pk
102          qkm2 = qkm1
103          qkm1 = qk
104          if abs(pk) > BIG:
105                  pkm2 *= BIGINV;
106                  pkm1 *= BIGINV;
107                  qkm2 *= BIGINV;
108                  qkm1 *= BIGINV;
109          if t <= MACHEP:
110              return ans * ax
111
112
113 -def _igam(a, x):
114      """Left tail of incomplete Gamma function:
115
116              inf.      k
117       a  -x   -       x
118      x  e     >   ----------
119               -     -
120              k=0   | (a+k+1)
121      """
122      # Compute  x**a * exp(-x) / Gamma(a)
123      ax = math.exp(a * math.log(x) - x - math.lgamma(a))
124
125      # Power series
126      r = a
127      c = 1.0
128      ans = 1.0
129
130      while True:
131          r += 1.0
132          c *= x/r
133          ans += c
134          if c / ans <= MACHEP:
135              return ans * ax / a
136
137
138  # --- Speed ---
139
140  #try:
141  #    from scipy.stats import chisqprob
142  #except ImportError:
143  #    pass
144
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu May 29 11:45:52 2014 http://epydoc.sourceforge.net