[frames] | no frames]

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

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:13 2011 http://epydoc.sourceforge.net