--- /dev/null
+###########################################################################
+# #
+# 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