import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammaincc

Compare limits Bayesian (two diff. prior) and frequentist
on counting experiments with background
No uncertainty on the background.

# Here are counts and backgrouns
N = [0,     2,    0,   2,   0, 2]
B = [0.5, 0.5,  2.0, 2.0, 3.5, 3.5]

To calculate the frequentist limit use result $\sum_{n=0}^N \frac{\mu^n e^{-\mu}}{N!} = \frac{1}{N!}\Gamma(N+1, \mu)$

where $\Gamma(a,b) = \int_b^{\infty} t^{a-1}e^{-t} dt$

for n,b in zip(N,B):

    print("N =",n, "and B =",b)
    # Frequentist first.
    s   = -b
    ds  = 0.01
    while gammaincc(n+1, s+b) > 0.05:
        s = s + ds
    print('Frequentist                  limit = ', round(s,1))
        
    # Now Bayesian. A bit brute force 
    # (25 is about infty for our purposes)
    sarray        = np.linspace(0, 25, 2500)
    integrand     = np.power((sarray+b),n) * np.exp(-(sarray+b))
    totalIntegral = integrand.sum()  # up to factor ds
    integral      = np.cumsum(integrand) / totalIntegral
    index         = (integral>0.95).argmax()  # index of limit
    s             = sarray[index]
    print('Bayesian flat prior          limit = ', round(s,1))
    
    integrand     = np.power((sarray+b),n) * np.exp(-(sarray+b)) / np.sqrt(sarray+b)
    totalIntegral = integrand.sum()  # up to factor ds
    integral      = np.cumsum(integrand) / totalIntegral
    index         = (integral>0.95).argmax()  # index of limit
    s             = sarray[index]
    print('Bayesian 1/sqrt(S+B) prior   limit = ', round(s,1))
    
    print(" ")
    
        

N = 0 and B = 0.5
Frequentist                  limit =  2.5
Bayesian flat prior          limit =  3.0
Bayesian 1/sqrt(S+B) prior   limit =  2.4
 
N = 2 and B = 0.5
Frequentist                  limit =  5.8
Bayesian flat prior          limit =  5.8
Bayesian 1/sqrt(S+B) prior   limit =  5.1
 
N = 0 and B = 2.0
Frequentist                  limit =  1.0
Bayesian flat prior          limit =  3.0
Bayesian 1/sqrt(S+B) prior   limit =  2.7
 
N = 2 and B = 2.0
Frequentist                  limit =  4.3
Bayesian flat prior          limit =  4.8
Bayesian 1/sqrt(S+B) prior   limit =  4.3
 
N = 0 and B = 3.5
Frequentist                  limit =  -0.5
Bayesian flat prior          limit =  3.0
Bayesian 1/sqrt(S+B) prior   limit =  2.7
 
N = 2 and B = 3.5
Frequentist                  limit =  2.8
Bayesian flat prior          limit =  4.3
Bayesian 1/sqrt(S+B) prior   limit =  3.9