Package Bio :: Package CodonAlign :: Module chisq
[hide private]
[frames] | no frames]

Source Code for Module Bio.CodonAlign.chisq

  1  """Python implementation of chisqprob, to avoid SciPy dependency. 
  2   
  3  Adapted from SciPy: scipy/special/cephes/{chdtr,igam}.c 
  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 21 Broadcasting rules apply. 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