Poisson Distribution Confidence Intervals in Scipy

2012 June 7

As happens often when doing transient searches, lately I’ve found myself wanting to compute confidence intervals in rate measurements from a small number of events. Clearly I want to be using the properties of the Poisson distribution, but how so, exactly?

The “Parameter Estimation” section of the above Wikipedia page gives some
hints but involves a few notation switches and also special functions that
don’t seem to be implemented in the handy
`scipy.special`

or
`scipy.stats`

`
modules. Once you track down the notation, though, it turns out to be really
easy. We want the inverse incomplete regularized Γ function.

Let’s say you observe *n* events in a period and want to compute the *k*
confidence interval on the true rate — that is, 0 < *k* ≤ 1, and *k* = 0.95
would be the equivalent of 2σ. Let *a* = 1 − *k*, *i.e.* 0.05. The lower bound
of the confidence interval, expressed as a potential number of events, is

scipy.special.gammaincinv(n, 0.5 * a)

and the upper bound is

scipy.special.gammaincinv(n + 1, 1 - 0.5 * a)

Conversion to rates is then just a matter of dividing by the relevant
duration, area, whatever. These apparently give the “exact” interval which may
be minutely more strict than is necessary. The halving of *a* is just because
the 95% confidence interval is made up of two tails of 2.5% each, so the
`gammaincinv`

function is really, once you chop through the obscurity, exactly
what you want.

So in
this example,
*n* is 14, *a* is 0.05, and we divide by 400 to get the interval bounds as
reported (0.019135 to 0.058724, in case that webpage disappears).

It took me a while to understand the numbers in the first table
on this page, but
they’re very helpful once you get over that: they give the 95% confidence
limits for a given number of observed events, when that number is small.
*I.e.*, if I observe one event over some period, my 95% confidence interval is
0.0253 to 5.5616 events. These numbers don’t quite make sense on their own,
since you can only get integral numbers of events in practice, but you can use
them to compute rate limits that *are* perfectly meaningful.

(Once you think about it, the relationship between the expressions for the two
ends of the interval is kind of cute. Because the Poisson distribution can
only take on integral values, the 2.5% upper limit to the value *n* must be
the 2.5% lower limit to the value *n* + 1, because there are no intervening
values that can be chosen.)