In [337]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import math
import heapq
In [334]:
MAX_ITER = 64
eps = 1.19209e-07
cp_n = 0
cp_p = 0
cnx = []
nx = []
def init_compute_possibility(n, p):
global cp_n, cp_p, cnx, nx
cp_n = n
cp_p = p
cnx = [0 for _ in range(n+1)]
nx = [0 for _ in range(n+1)]
cnx[0] = 1.0
for i in range(1, n+1):
cnx[i] = cnx[i-1] * (n - i + 1) / i
for i in range(n+1):
nx[i] = cnx[i] * p**i * (1-p)**(n-i)
def compute_possibility(n, p, x):
global cp_n, cp_p, nx
if cp_n != n or cp_p != p:
init_compute_possibility(n, p)
return nx[x]
#def compute_possibility(n, p, x):
# s = 1.0
# for i in range(x):
# s *= (n - i) / (i + 1)
# s *= (p**x) * ((1-p)**(n-x))
# return s
def bs_next_helper(n, c, t, L, m, p):
s = 0.0
for i in range(m+1):
#print(p, L[i], L[i] + t)
if p < L[i] or p > L[i] + t:
continue
#print("s", s, compute_possibility(n, p, m))
s += compute_possibility(n, p, i)
if s >= c:
return True
return False
#print(bs_next_helper(10, 0.9, 1.0, [0.0,0.0], 0, 0.0001))
def bs_next(n, c, t, L, i, hq):
if len(hq) == 0:
return 1.0
beg = heapq.heappop(hq)
while len(hq) > 0:
end = heapq.heappop(hq)
#print(end, bs_next_helper(n, c, t, L, i, end))
if bs_next_helper(n, c, t, L, i, end):
beg = end
continue
if bs_next_helper(n, c, t, L, i, end - eps):
beg = end - eps
else:
for _ in range(MAX_ITER):
mid = (beg + end) / 2.0
if mid == beg:
break
#print("mid:", mid, L[0], t, c, bs_next_helper(n, c, t, L, i+1, mid))
if bs_next_helper(n, c, t, L, i, mid):
beg = mid
else:
end = mid
heapq.heappush(hq, beg)
#heapq.heappush(hq, min(1.0, beg+t))
heapq.heappush(hq, min(1.0, beg+t+eps))
heapq.heappush(hq, end)
return beg
return 1.0
def bci_helper(n, c, t, L):
for i in range(n+1):
L[i] = 0.0
#L[0] = 0.0
#L[1] = 1.0 - math.pow(c, 1/n)
hq = [i / n for i in range(n+1)]
heapq.heappush(hq, t)
for i in range(n):
L[i+1] = bs_next(n, c, t, L, i, hq)
#print(hq)
#tmp = test_bci(n, L, t, 0.01)
#if tmp >= c:
# return True
#return False
while len(hq) > 0:
p = heapq.heappop(hq)
if compute_confidence(n, L, t, p) < c:
return False
if p == 1.0:
break
return True
def compute_binomail_confidence_interval(n, c):
L = [0.0 for _ in range(n+1)]
beg = 0.0
end = 1.0
for _ in range(MAX_ITER):
mid = (beg + end) / 2.0
if mid == beg:
break
ok = bci_helper(n, c, mid, L)
#print(mid, ok)
if ok:
end = mid
else:
beg = mid
bci_helper(n, c, end, L)
return L, [min(1.0, l + end) for l in L], end
def compute_confidence(n, L, t, p):
s = 0.0
for i in range(n+1):
if L[i] <= p and p <= L[i] + t:
s += compute_possibility(n, p, i)
return s
def test_bci(n, L, t, step):
ret = compute_confidence(n, L, t, 1.0)
ind = 1.0
p = 0
while p <= 1.0:
tmp = compute_confidence(n, L, t, p)
if tmp < ret:
ret = tmp
ind = p
p += step
return ret, ind
def compute_and_test(n, c, step):
L, R, t = compute_binomail_confidence_interval(n, c)
ret, ret_p = test_bci(n, L, t, step)
print("result:", c, ret, ret_p)
return L, R, t
#compute_binomail_confidence_interval(10, 0.9)
#compute_and_test(10, 0.9, 0.01)
In [335]:
n = 20
c = 0.7
steps = 100
L, R, t = compute_and_test(n, c, (1.0 / steps))
print("interval:", t)
x = [i * (1.0 / steps) for i in range(steps+1)]
y = [compute_confidence(n, L, t, i) for i in x]
#print(x, y)
plt.plot(x,y)
#for i,j in zip(x, y):
# print(i,j)
Out[335]:
In [336]:
print("n:", n, ", confidence:", c, ", interval:", t)
for i, l in enumerate(L):
print(i, [l, min(1.0, l + t)], i / n)
#compute_confidence(n, L, t, 0.713201540802663)