X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fcisstat%2Fhelper.py;fp=cistematic%2Fcisstat%2Fhelper.py;h=7abf82712580a173c9c3105772bca556a38cad6d;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181

diff --git a/cistematic/cisstat/helper.py b/cistematic/cisstat/helper.py
new file mode 100644
index 0000000..7abf827
--- /dev/null
+++ b/cistematic/cisstat/helper.py
@@ -0,0 +1,127 @@
+###########################################################################
+#                                                                         #
+# C O P Y R I G H T   N O T I C E                                         #
+#  Copyright (c) 2003-10 by:                                              #
+#    * California Institute of Technology                                 #
+#                                                                         #
+#    All Rights Reserved.                                                 #
+#                                                                         #
+# Permission is hereby granted, free of charge, to any person             #
+# obtaining a copy of this software and associated documentation files    #
+# (the "Software"), to deal in the Software without restriction,          #
+# including without limitation the rights to use, copy, modify, merge,    #
+# publish, distribute, sublicense, and/or sell copies of the Software,    #
+# and to permit persons to whom the Software is furnished to do so,       #
+# subject to the following conditions:                                    #
+#                                                                         #
+# The above copyright notice and this permission notice shall be          #
+# included in all copies or substantial portions of the Software.         #
+#                                                                         #
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
+# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
+# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
+# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
+# SOFTWARE.                                                               #
+###########################################################################
+#
+# special functions that are helpful for our statistics
+# gammaln, gamma, factorial, factln, bico, and beta are inspired by:
+#           Numerical Recipes in C
+#           by Press, Flannery, Teukolsky, Vetterling
+#
+from math import log, exp, floor
+
+def gammaln(Z):
+    cof = [ 76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5]
+
+    x = float(Z) - 1.0
+    temp = x + 5.5
+    ser = 1.0
+    j = 0
+
+    temp -= (x + 0.5) * log(temp)
+    for j in range(6):
+        x += 1.0
+        ser += cof[j] / x
+
+    return -temp + log(2.50662827465 * ser)
+
+
+def gammaln2(Z):
+    cof = [676.5203681218835, -1259.139216722289, 771.3234287757674, -176.6150291498386, 
+           12.50734324009056, -0.1385710331296526, 0.9934937113930748e-5, 0.1659470187408462e-6]
+
+    ser = 0.9999999999995183
+    x = float(Z) - 1.0
+
+    for j in range(7):
+        x += 1.0
+        ser += cof[j] / x
+
+    result = log(ser) - 5.58106146679532777 - float(Z) + (float(Z) - 0.5) * log(float(Z) + 6.5)
+
+    return result
+
+
+def gamma(Z):
+    return round(exp(gammaln(Z)), 8)
+
+
+def gamma2(Z):
+    return round(exp(gammaln2(Z)), 8)
+
+
+def factorial(N):
+    ntop = 4
+    arr = [1.0, 1.0, 2.0, 6.0, 24.0]
+
+    if N < 0:
+        print "error: factorial only works for non-negative integers"
+        return -1
+    elif N > 32:
+        return gamma(N + 1)
+    else:
+        while ntop < N:
+            j = ntop
+            ntop += 1
+            arr.append(arr[j] * ntop)
+        return arr[N]
+
+
+def factln(N):
+    if N < 0:
+        print "error: factorial only works for non-negative integers"
+        return -1
+    if N <= 1:
+        return 0.0
+    else:
+        return gammaln2(N + 1.0)
+
+
+def bicoln(n, k):
+    return factln(n) - factln(k) - factln(n - k)
+
+
+def bico(n, k):
+    return floor(0.5 + exp(bicoln(n, k)))
+
+
+def beta(z, w):
+    return exp(gammaln(z) + gammaln(w) - gammaln(z + w))
+
+
+def hyperGeometric(N, r, n, x):
+    """ hyperGeometric(N, r, n, x):
+                N   is population size
+                r   is the number in the population identified as a success
+                n   is the sample size
+                x   is the number in the sample identified as a success
+    """
+    N = float(N)
+    r = float(r)
+    n = float(n)
+    x = float(x)
+    return exp(bicoln(r, x) + bicoln(N - r, n - x) - bicoln(N, n))
\ No newline at end of file