Obtain correct error bars from observed counts (fraction in sample, detected events, etc.)

poissonerr.py

"""
If you ever made a plot of "fractions" or "rates" with symmetric error 
bars, like in the plot shown, this mini-tutorial is for you. Here is how
 to compute correct error bars, so that uncertainties in the fractions 
do not go below 0% or above 100%.

 If you have a histogram for 
instance, and you detected k objects in a given bin. What is the rate 
underlying at quantiles q=0.1, 0.5(median), 0.9?
"""

import scipy.special
k = 10
q = 0.5 # median
rate = scipy.special.gammaincinv(k + 1, q)

# Gives 10.66 (q=10%-90% -> 7 - 15.4). Note the asymmetric error bars.

"""
If you have n=20 objects and k=17 of them are of a certain class (e.g. AGN). What is the success rate to be in that class?
"""

n = 20
k = 17
q = 0.5 # median
scipy.special.betaincinv(k+1, n+1-k, q)
# Gives a success rate of 83% (71% - 91%).

"""
The incomplete gamma function is underlying the definition of the 
Poisson probability distribution. Inverting the cumulative distribution 
to go from quantiles to a rate is what gammaincinv does.

The incomplete beta function is underlying the definition of the Binomial 
probability distribution. Inverting the cumulative distribution to go 
from quantiles to a success rate is what betaincinv does.

Further reading:
 https://en.wikipedia.org/wiki/Poisson_distribution (CDF in the box)
 https://en.wikipedia.org/wiki/Binomial_distribution
"""