from math import factorial
import numpy as np
from numpy.polynomial import Polynomial as Po
from random import random, randint
#This should run in Python 3 where big integers allowed, otherwise it should be swapped
#out with a smarter implentation, which I couldn't find in scipy for some reason.
def binomial(n, k): return factorial(n)/(factorial(k)*factorial(n-k))
#the original cumulative binomial function we are going to invert.
def cumulative_binomial(p, n): return sum([binomial(n, k)*p**k*(1-p)**(n-k) for k in range(int(n/2), n+1)])
# rephrasing of the same function, for mathematical convenience to see how the below works
def r_cumulative_binomial(p, n):
return sum([sum([binomial(n, k)*binomial(n-k, l)*(-1)**l*p**(l+k) for l in range(n-k+1)]) for k in range(int(n/2), n+1)])
#this returns the cumulative binomial as an array with which we can create a numpy polynomial, summed from s to n
def polynomial_form(n, s):
degree = 0 #holds the highest exponent, to know how to size numpy array
ll = [] #holds the many terms, as tuple (coefficient, integer power)
for k in range(s, n+1):
for l in range(n-k+1):
ll.append((binomial(n, k)*binomial(n-k, l)*(-1)**l, l+k)) #tuple of current term
if l+k > degree: degree= l+k #if biggest seen so far, set degree
out = np.array([complex(i) for i in np.zeros(degree+1)]) #make the array that will become the polynomial
for i in ll: out[i[1]] += i[0] #set its terms
return out
#this solves it, given epsilon, s, and n
def inverse_cumulative_binomial(eps, s, n):
pp = polynomial_form(n, s) #get the polynomial term
pp[0] -= eps #subtract epsilon from it
pp = Po(pp) #covert to numpy polynomial
return pp.roots() #return all its roots
#this filters the solution into real positive numbers
def get_real_positive(l, err=.00001): return [i.real for i in l if abs(i.imag)<err and i.real > 0]
#test it with random p
test_p = randint(1,10)/1000
test_n = randint(4,16)*2
test_s = int(test_n/2)
sol = cumulative_binomial(test_p, test_n)
print ("cumulative_binomial for p=", test_p, " n=",str(test_n),": \n", sol, "\n for same r_cumubino:\n", r_cumulative_binomial(test_p, test_n))ans = inverse_cumulative_binomial(sol, int(test_n/2), test_n)
print ("real positive solutions\n", get_real_positive(ans))
print ("full solution set\n", ans)
