1 ###########################################################################
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 #
7 # All Rights Reserved. #
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: #
17 # The above copyright notice and this permission notice shall be #
18 # included in all copies or substantial portions of the Software. #
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 #
28 ###########################################################################
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
35 from math import log, exp, floor
38 cof = [ 76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5]
45 temp -= (x + 0.5) * log(temp)
50 return -temp + log(2.50662827465 * ser)
54 cof = [676.5203681218835, -1259.139216722289, 771.3234287757674, -176.6150291498386,
55 12.50734324009056, -0.1385710331296526, 0.9934937113930748e-5, 0.1659470187408462e-6]
57 ser = 0.9999999999995183
64 result = log(ser) - 5.58106146679532777 - float(Z) + (float(Z) - 0.5) * log(float(Z) + 6.5)
70 return round(exp(gammaln(Z)), 8)
74 return round(exp(gammaln2(Z)), 8)
79 arr = [1.0, 1.0, 2.0, 6.0, 24.0]
82 print "error: factorial only works for non-negative integers"
90 arr.append(arr[j] * ntop)
96 print "error: factorial only works for non-negative integers"
101 return gammaln2(N + 1.0)
105 return factln(n) - factln(k) - factln(n - k)
109 return floor(0.5 + exp(bicoln(n, k)))
113 return exp(gammaln(z) + gammaln(w) - gammaln(z + w))
116 def hyperGeometric(N, r, n, x):
117 """ hyperGeometric(N, r, n, x):
119 r is the number in the population identified as a success
121 x is the number in the sample identified as a success
127 return exp(bicoln(r, x) + bicoln(N - r, n - x) - bicoln(N, n))