Calculate exact Binomial Limits

limits.py

import scipy.special
import scipy.stats

def binomial_limit(n, k, sigma=1):
    """
    Calculates the upper and lower limit for the probability p of a binomial distribution
    if an experiment yielded k successes for n trials.
    The confidence level for the limits is given in sigmas of the gaussian distribution.
    In contrast to other methods (see https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval),
    this method is exact in the sense that there is no approximation involved for the binomial distribution.
    @param n number of trials
    @param k number of successes
    @param sigma confidence level expressed in sigmas of a gaussian distribution (e.g. 2 sigma ~ 95% CL)
    """
    alpha = 2*(1-scipy.stats.norm.cdf(sigma))
    upper_limit = 1-scipy.special.betaincinv(n-k+1, k+1, (1-k/n)*alpha)
    lower_limit = 1-scipy.special.betaincinv(n-k+1, k+1, 1-k/n*alpha)
    return (upper_limit, lower_limit)

# n=10 measurements and no successes k=0 (only upper limit, lower limit is always 0)
binomial_limit(10, 0, 1)

# n=100 measurements and 50 successes k=50
binomial_limit(100, 50, 1)

# Rule of three Upper limit is 3/n for n trial with no successes @ 95 CL = 2 sigma
binomial_limit(100, 0, 2)

import numpy as np
# Coverage
N =10000
n = np.random.randint(1,100, size=N)
p =  np.random.random(size=N)
k = scipy.stats.binom.rvs(n, p)
u, l = binomial_limit(n, k)
coverage = ((u > p) & (l < p)).mean()
print(coverage)