cistematic/cisstat/helper.py

   1 ###########################################################################
   2 #                                                                         #
   3 # C O P Y R I G H T   N O T I C E                                         #
   4 #  Copyright (c) 2003-10 by:                                              #
   5 #    * California Institute of Technology                                 #
   6 #                                                                         #
   7 #    All Rights Reserved.                                                 #
   8 #                                                                         #
   9 # Permission is hereby granted, free of charge, to any person             #
  10 # obtaining a copy of this software and associated documentation files    #
  11 # (the "Software"), to deal in the Software without restriction,          #
  12 # including without limitation the rights to use, copy, modify, merge,    #
  13 # publish, distribute, sublicense, and/or sell copies of the Software,    #
  14 # and to permit persons to whom the Software is furnished to do so,       #
  15 # subject to the following conditions:                                    #
  16 #                                                                         #
  17 # The above copyright notice and this permission notice shall be          #
  18 # included in all copies or substantial portions of the Software.         #
  19 #                                                                         #
  20 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
  21 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
  22 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
  23 # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
  24 # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
  25 # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
  26 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
  27 # SOFTWARE.                                                               #
  28 ###########################################################################
  29 #
  30 # special functions that are helpful for our statistics
  31 # gammaln, gamma, factorial, factln, bico, and beta are inspired by:
  32 #           Numerical Recipes in C
  33 #           by Press, Flannery, Teukolsky, Vetterling
  34 #
  35 from math import log, exp, floor
  36
  37 def gammaln(Z):
  38     cof = [ 76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5]
  39
  40     x = float(Z) - 1.0
  41     temp = x + 5.5
  42     ser = 1.0
  43     j = 0
  44
  45     temp -= (x + 0.5) * log(temp)
  46     for j in range(6):
  47         x += 1.0
  48         ser += cof[j] / x
  49
  50     return -temp + log(2.50662827465 * ser)
  51
  52
  53 def gammaln2(Z):
  54     cof = [676.5203681218835, -1259.139216722289, 771.3234287757674, -176.6150291498386,
  55            12.50734324009056, -0.1385710331296526, 0.9934937113930748e-5, 0.1659470187408462e-6]
  56
  57     ser = 0.9999999999995183
  58     x = float(Z) - 1.0
  59
  60     for j in range(7):
  61         x += 1.0
  62         ser += cof[j] / x
  63
  64     result = log(ser) - 5.58106146679532777 - float(Z) + (float(Z) - 0.5) * log(float(Z) + 6.5)
  65
  66     return result
  67
  68
  69 def gamma(Z):
  70     return round(exp(gammaln(Z)), 8)
  71
  72
  73 def gamma2(Z):
  74     return round(exp(gammaln2(Z)), 8)
  75
  76
  77 def factorial(N):
  78     ntop = 4
  79     arr = [1.0, 1.0, 2.0, 6.0, 24.0]
  80
  81     if N < 0:
  82         print "error: factorial only works for non-negative integers"
  83         return -1
  84     elif N > 32:
  85         return gamma(N + 1)
  86     else:
  87         while ntop < N:
  88             j = ntop
  89             ntop += 1
  90             arr.append(arr[j] * ntop)
  91         return arr[N]
  92
  93
  94 def factln(N):
  95     if N < 0:
  96         print "error: factorial only works for non-negative integers"
  97         return -1
  98     if N <= 1:
  99         return 0.0
 100     else:
 101         return gammaln2(N + 1.0)
 102
 103
 104 def bicoln(n, k):
 105     return factln(n) - factln(k) - factln(n - k)
 106
 107
 108 def bico(n, k):
 109     return floor(0.5 + exp(bicoln(n, k)))
 110
 111
 112 def beta(z, w):
 113     return exp(gammaln(z) + gammaln(w) - gammaln(z + w))
 114
 115
 116 def hyperGeometric(N, r, n, x):
 117     """ hyperGeometric(N, r, n, x):
 118                 N   is population size
 119                 r   is the number in the population identified as a success
 120                 n   is the sample size
 121                 x   is the number in the sample identified as a success
 122     """
 123     N = float(N)
 124     r = float(r)
 125     n = float(n)
 126     x = float(x)
 127     return exp(bicoln(r, x) + bicoln(N - r, n - x) - bicoln(N, n))