converted some parts of binomial_conf_interval to logarithmic scale to avoid underflo...
authorForrest Voight <forrest@forre.st>
Thu, 9 Feb 2012 19:50:11 +0000 (14:50 -0500)
committerForrest Voight <forrest@forre.st>
Thu, 9 Feb 2012 19:50:11 +0000 (14:50 -0500)
p2pool/util/math.py

index 614a0f2..0d66c78 100644 (file)
@@ -129,14 +129,14 @@ else:
         if n == 0:
             left = random.random()*(1 - conf)
             return left, left + conf
-        b = special.beta(x+1, n-x+1)
+        bl = float(special.betaln(x+1, n-x+1))
         def f(left_a):
-            left, right = max(1e-8, special.betaincinv(x+1, n-x+1, left_a)), min(1-1e-8, special.betaincinv(x+1, n-x+1, left_a + conf))
-            top = right**(x+1) * (1-right)**(n-x+1) * left*(1-left) - left**(x+1) * (1-left)**(n-x+1) * right * (1-right)
+            left, right = max(1e-8, float(special.betaincinv(x+1, n-x+1, left_a))), min(1-1e-8, float(special.betaincinv(x+1, n-x+1, left_a + conf)))
+            top = math.exp(math.log(right)*(x+1) + math.log(1-right)*(n-x+1) + math.log(left) + math.log(1-left) - bl) - math.exp(math.log(left)*(x+1) + math.log(1-left)*(n-x+1) + math.log(right) + math.log(1-right) - bl)
             bottom = (x - n*right)*left*(1-left) - (x - n*left)*right*(1-right)
-            return top/bottom/b
+            return top/bottom
         left_a = find_root(f, (1-conf)/2, bounds=(0, 1-conf))
-        return special.betaincinv(x+1, n-x+1, left_a), special.betaincinv(x+1, n-x+1, left_a + conf)
+        return float(special.betaincinv(x+1, n-x+1, left_a)), float(special.betaincinv(x+1, n-x+1, left_a + conf))
 
 def binomial_conf_center_radius(x, n, conf=0.95):
     assert 0 <= x <= n and 0 <= conf < 1