--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/server/stats.py	Wed Jul 16 13:49:50 2014 -0400
@@ -0,0 +1,98 @@
+#make sure default is floating point division
+from __future__ import division
+from math import log, exp
+
+def ln_gamma(a):
+    c = 2.506628274631
+    s = [0] * 6
+    total = 0
+    temp = 0
+    s[0] =  76.180091729406 / a
+    s[1] = -86.505320327112 / (a + 1)
+    s[2] =  24.014098222230 / (a + 2)
+    s[3] =  -1.231739516140 / (a + 3)
+    s[4] =   0.001208580030 / (a + 4)
+    s[5] =  -0.000005363820 / (a + 5)
+    total =  1.000000000178
+    for i in range(6):
+        total += s[i]
+        temp = (a - 0.5) * log(a + 4.5) - (a + 4.5) + log(c * total)
+    return temp
+
+def ln_beta(a, b):
+    return ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b)
+    
+
+def incomplete_beta(a, b, x):
+    tiny = 1.0e-10                      # convergence parameter 
+    p = [0] * 3
+    q = [0] * 3
+
+    if x > (a + 1) / (a + b + 1):   # to accelerate convergence   
+        swap = True                     # complement x and swap a & b 
+        x    = 1 - x
+        t    = a
+        a    = b
+        b    = t
+    else:                                # do nothing 
+        swap = False
+    if x > 0:
+        factor = exp(a * log(x) + b * log(1 - x) - ln_beta(a,b)) / a
+    else:
+        factor = 0
+    p[0] = 0
+    q[0] = 1
+    p[1] = 1
+    q[1] = 1
+    f = p[1] / q[1]
+    n = 0
+
+    while True:                               # recursively generate the continued 
+        g = f                           # fraction 'f' until two consecutive 
+        n += 1                             # values are small                   
+        if ((n%2) > 0)  :
+            t = (n - 1) / 2
+            c = -(a + t) * (a + b + t) * x / ((a + n - 1) * (a + n))
+
+        else :
+            t = n / 2
+            c = t * (b - t) * x / ((a + n - 1) * (a + n))
+        p[2] = p[1] + c * p[0]
+        q[2] = q[1] + c * q[0]
+        if q[2] != 0:                # rescale to avoid overflow 
+            p[0] = p[1] / q[2]
+            q[0] = q[1] / q[2]
+            p[1] = p[2] / q[2]
+            q[1] = 1.0
+            f = p[1]
+
+        if abs(f - g) < tiny and q[1] == 1:
+            break
+  
+    if swap:
+        return 1.0 - factor * f
+
+    return factor * f
+
+def student_pdf(n, x):
+    t = -0.5 * (n + 1) * log(1 + ((x**2) / n)) - ln_beta(0.5, n / 2.0)
+    return exp(t) / n**0.5
+
+def student_cdf(n, x):
+    t = (x ** 2) / (n + x ** 2)
+    s = incomplete_beta(0.5, n * 0.5, t)
+    if x >= 0:
+        return 0.5 * (1.0 + s)
+    return 0.5 * (1.0 - s)
+
+def student_idf(n, u):
+    tiny = 1.0e-10
+    x = 0
+    while True:
+        t = x
+        x = t + (u - student_cdf(n, t)) / student_pdf(n, t)
+        if abs(x - t) < tiny:
+            break
+    return x
+
+