Source Code for Module gmisclib.ftest

  1  # Copyright (c) 1999-2000 Gary Strangman; All Rights Reserved. 
  2  # 
  3  # This software is distributable under the terms of the GNU 
  4  # General Public License (GPL) v2, the text of which can be found at 
  5  # http://www.gnu.org/copyleft/gpl.html. Installing, importing or otherwise 
  6  # using this module constitutes acceptance of the terms of this License. 
  7  # 
  8  # Disclaimer 
  9  #  
 10  # This software is provided "as-is".  There are no expressed or implied 
 11  # warranties of any kind, including, but not limited to, the warranties 
 12  # of merchantability and fittness for a given application.  In no event 
 13  # shall Gary Strangman be liable for any direct, indirect, incidental, 
 14  # special, exemplary or consequential damages (including, but not limited 
 15  # to, loss of use, data or profits, or business interruption) however 
 16  # caused and on any theory of liability, whether in contract, strict 
 17  # liability or tort (including negligence or otherwise) arising in any way 
 18  # out of the use of this software, even if advised of the possibility of 
 19  # such damage. 
 20  # 
 21  # Comments and/or additions are welcome (send e-mail to: 
 22  # strang@nmr.mgh.harvard.edu). 
 23  #  
 24  import math 
 25  import numpy as N 
 26   
 27 -def fprob (dfnum, dfden, F): 
 28      """ 
 29  Returns the (1-tailed) significance level (p-value) of an F 
 30  statistic given the degrees of freedom for the numerator (dfR-dfF) and 
 31  the degrees of freedom for the denominator (dfF). 
 32   
 33  Usage:   lfprob(dfnum, dfden, F)   where usually dfnum=dfbn, dfden=dfwn 
 34  """ 
 35      p = betai(0.5*dfden, 0.5*dfnum, dfden/float(dfden+dfnum*F)) 
 36      return p 
 37   
 38   
 39 -def betai(a,b,x): 
 40      """Returns the incomplete beta function:: 
 41   
 42          I-sub-x(a,b) = 1/B(a,b)*(Integral(0,x) of t^(a-1)(1-t)^(b-1) dt) 
 43   
 44          where a,b>0 and B(a,b) = G(a)*G(b)/(G(a+b)) where G(a) is the gamma 
 45          function of a.  The continued fraction formulation is implemented here, 
 46          using the betacf function.  (Adapted from: Numerical Recipies in C.) 
 47   
 48          Usage:   lbetai(a,b,x) 
 49  """ 
 50      if (x<0.0 or x>1.0): 
 51          raise ValueError, 'Bad x in lbetai' 
 52      if (x==0.0 or x==1.0): 
 53          bt = 0.0 
 54      else: 
 55          bt = math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*math.log(x)+b* 
 56                        math.log(1.0-x)) 
 57      if (x<(a+1.0)/(a+b+2.0)): 
 58          return bt*betacf(a,b,x)/float(a) 
 59      else: 
 60          return 1.0-bt*betacf(b,a,1.0-x)/float(b) 
 61   
 62   
 63   
 64   
 65 -def betacf(a,b,x): 
 66      """ 
 67  This function evaluates the continued fraction form of the incomplete 
 68  Beta function, betai.  (Adapted from: Numerical Recipies in C.) 
 69   
 70  Usage:   lbetacf(a,b,x) 
 71  """ 
 72      ITMAX = 200 
 73      EPS = 3.0e-7 
 74   
 75      bm = az = am = 1.0 
 76      qab = a+b 
 77      qap = a+1.0 
 78      qam = a-1.0 
 79      bz = 1.0-qab*x/qap 
 80      for i in range(ITMAX+1): 
 81          em = float(i+1) 
 82          tem = em + em 
 83          d = em*(b-em)*x/((qam+tem)*(a+tem)) 
 84          ap = az + d*am 
 85          bp = bz+d*bm 
 86          d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem)) 
 87          app = ap+d*az 
 88          bpp = bp+d*bz 
 89          aold = az 
 90          am = ap/bpp 
 91          bm = bp/bpp 
 92          az = app/bpp 
 93          bz = 1.0 
 94          if (abs(az-aold)<(EPS*abs(az))): 
 95              return az 
 96      print 'a or b too big, or ITMAX too small in Betacf.' 
 97   
 98   
 99 -def gammln(xx): 
100      """Returns the gamma function of xx. 
101          Gamma(z) = Integral(0,infinity) of t^(z-1)exp(-t) dt. 
102          (Adapted from: Numerical Recipies in C.) 
103          @param xx: float 
104          @rtype: float 
105          """ 
106   
107      coeff = [76.18009173, -86.50532033, 24.01409822, -1.231739516, 
108               0.120858003e-2, -0.536382e-5] 
109      x = xx - 1.0 
110      tmp = x + 5.5 
111      tmp = tmp - (x+0.5)*math.log(tmp) 
112      ser = 1.0 
113      for j in range(len(coeff)): 
114          x = x + 1 
115          ser = ser + coeff[j]/x 
116      return -tmp + math.log(2.50662827465*ser) 
117   
118   
119 -def agammln(xx): 
120      """Returns the gamma function of xx. 
121      C{Gamma(z) = Integral(0,infinity) of t^(z-1)exp(-t) dt}. 
122      Adapted from: Numerical Recipies in C.  Can handle multiple dims ... but 
123      probably doesn't normally have to. 
124   
125      Usage:   agammln(xx) 
126      """ 
127      coeff = [76.18009173, -86.50532033, 24.01409822, -1.231739516, 
128               0.120858003e-2, -0.536382e-5] 
129      x = xx - 1.0 
130      tmp = x + 5.5 
131      tmp = tmp - (x+0.5)*N.log(tmp) 
132      ser = 1.0 
133      for j in range(len(coeff)): 
134          x = x + 1 
135          ser = ser + coeff[j]/x 
136      return -tmp + N.log(2.50662827465*ser) 
137