Poisson confidence interval with numpy

Question

I am trying to put Poisson error strings on a histogram that I create using matplotlib, but I cannot find a numpy function that will give me a 95% confidence interval suggesting Poisson data. Ideally, the solution does not depend on scipy, but everything will work. Is there such a function? I learned a lot about downloading, but in my case it seems a little too much.

Answer 1

Using scipy.stats.poisson and the interval method:

 >>> scipy.stats.poisson.interval(0.95, [10, 20, 30]) (array([ 4., 12., 20.]), array([ 17., 29., 41.])) 

Although calculating the Poisson distribution for non-integer values ​​has limited meaning, it is possible to calculate the exact confidence intervals requested by the OP, which can be done as follows:

 >>> data = np.array([10, 20, 30]) >>> scipy.stats.poisson.interval(0.95, data) (array([ 4., 12., 20.]), array([ 17., 29., 41.])) >>> np.array(scipy.stats.chi2.interval(.95, 2 * data)) / 2 - 1 array([[ 3.7953887 , 11.21651959, 19.24087402], [ 16.08480345, 28.67085357, 40.64883744]]) 

---

You can also use the ppf method:

 >>> data = np.array([10, 20, 30]) >>> scipy.stats.poisson.ppf([0.025, 0.975], data[:, None]) array([[ 4., 17.], [ 12., 29.], [ 20., 41.]]) 

But since the distribution is discrete, the returned values ​​will be integers, and the confidence interval will not be 95%:

 >>> scipy.stats.poisson.ppf([0.025, 0.975], 10) array([ 4., 17.]) >>> scipy.stats.poisson.cdf([4, 17], 10) array([ 0.02925269, 0.98572239]) 

Answer 2

In the end, I wrote my own function based on some properties that I found on Wikipedia .

 def poisson_interval(k, alpha=0.05): """ uses chisquared info to get the poisson interval. Uses scipy.stats (imports in function). """ from scipy.stats import chi2 a = alpha low, high = (chi2.ppf(a/2, 2*k) / 2, chi2.ppf(1-a/2, 2*k + 2) / 2) if k == 0: low = 0.0 return low, high 

This returns continuous (rather than discrete) boundaries, which is more standard in my area.

Answer 3

This problem is much in astronomy (my field!), And this article is a reference to these confidence intervals: Gehrels 1980

It has a lot of mathematics for an arbitrary confidence interval with Poisson statistics, but for a 95% two-sided confidence interval (corresponding to a 2-sigma Gaussian confidence interval or S = 2 in the context of this paper) some simple analytical formulas for the upper and lower confidence limits when measuring N events:

 upper = N + 2. * np.sqrt(N + 1) + 4. / 3. lower = N * (1. - 1. / (9. * N) - 2. / (3. * np.sqrt(N))) ** 3. 

where I already posted them in Python format. All you need is numpy or your other favorite square root module. Keep in mind that they will give you upper and lower limits for events, not +/- values. You simply subtract N from both to get them.

Consult paper for the accuracy of these formulas for the required confidence interval, but they should be more than accurate enough for most practical applications.