In [64]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import math
import heapq
In [65]:
MAX_ITER = 64
eps = 1.19209e-07
cp_n = 0
cnx = []
nx = {}
#range_array = []
def get_possibilities(n, p):
global cp_n, cnx, nx
if cp_n != n:
cp_n = n
cnx = np.zeros(n+1)
cnx[0] = 1.0
for i in range(1, n+1):
cnx[i] = cnx[i-1] * (n - i + 1) / i
nx = {}
if p not in nx:
ps = np.zeros(n+1)
for i in range(n+1):
ps[i] = cnx[i] * (p**i) * ((1-p)**(n-i))
nx[p] = ps
return nx[p]
def compute_possibility(n, p, x):
nx = get_possibilities(n, p)
return nx[x]
def compute_confidence(n, L, t, p):
nx = get_possibilities(n, p)
s = np.sum(nx[(L <= p) & ((L + t) >= p)])
return s
def bs_next_helper(n, c, t, L, m, p):
nx = get_possibilities(n, p)
nx = nx[:m+1]
Lm = L[:m+1]
s = np.sum(nx[(Lm <= p) & ((Lm + t) >= p)])
#s = np.sum(nx[(L <= p) & ((L + t) >= p)])
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+eps))
L[i+1] = beg
if not bs_next_helper(n, c, t, L, i+1, end):
heapq.heappush(hq, end)
return beg
return 1.0
def bci_helper(n, c, t, L):
#L -= 3.0
L[0] = 0.0
# wrong
#hq = [0.0, 1.0]
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):
#global range_array
#range_array = np.array(range(n))
L = np.zeros(n+1)
beg = 0.0
end = 1.0
for i in range(MAX_ITER):
mid = (beg + end) / 2.0
if mid == beg:
break
ok = bci_helper(n, c, mid, L)
print(i, mid, ok)
if ok:
end = mid
else:
beg = mid
bci_helper(n, c, end, L)
return L, L + end, end
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)
# %time L, R, t = compute_and_test(200, 0.99, 0.01)
In [66]:
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[66]:
In [67]:
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)