Weighted Poisson Uncertainties¶

Poisson distributions turn up a lot in high-energy astronomy. If you take the photons from a source and sort them into phase bins, the number in each phase bin will be Poisson distributed, with an expected value corresponding to the source's emission. This is something we often do, for example to look for evidence that the source's emission varies: in addition to various statistical tests, we plot a histogram of photon arrival phases. Part of plotting such a histogram is plotting error bars on each bin, to give the reader an indication of the likelihood that such variations are real.

The error bars on each bin are supposed to indicate the uncertainty on the value in that bin. Now, a Poisson random variable with expected value $\lambda$ has standard deviation $\sqrt{\lambda}$, but of course we don't actually know $\lambda$. We can estimate $\lambda$ as the number $n$ of observed photons (in fact this is an unbiased estimate and about as good as we can hope for from a single measurement). It turns out that if $n$ is decently large, $\sqrt{n}$ is a decent estimate for $\sqrt{\lambda}$. It's not unbiased, and it's definitely wrong if $n$ is small, but in that case the distribution isn't much like a normal distribution anyway, and so you enter realms of messy statistics. For histograms we try to keep the number $n$ in each bin pretty high anyway. So $\sqrt{n}$ is a decent value for the error bar in each bin.

We can even check the reasonableness of the error estimate. If we pick a $\lambda$ and generate Poisson samples $n_i$, then for each $n_i$ we can compute $z_i = (\lambda-n_i)/\sqrt{n_i}$, which should be distributed as a normal distribution with mean zero and unit standard deviation. So let's test that.

Not too bad. It's going to flunk a Kolmogorov-Smirnov test because it's discrete, and we can't really do much about that without going to higher $n$ where we expect the approximation to be better anyway. Well, I guess we could choose the $\lambda$ independently for each trial, but that's getting really very messy. And it's not what I'm really interested in anyway: the results so far are straight from the textbook.

What interests me right now is a more tricky variation on the problem. With the Fermi space telescope, the tool gtsrcprob gives us a probability for each photon that it's actually from the source (rather than from some nearby contaminating source). So what we want to plot in a histogram is not the number of photons, but the total probability, the expected number of photons. What's the uncertainty on that?

It's not entirely clear how best to pose the problem mathematically. But let's try this: there is some unknown distribution of photon weights. To find the expected number of photons in each bin, we draw a number of photons from a Poisson distribution, then we draw the weight from our (unknown) distribution.

Taking this as our mathematical model, what's the standard deviation of the result? Well, first let's get the mean: if the Poisson distribution has mean $\lambda$, and the unknown distribution has mean $\mu$, the expected value will be

$$E(T) = \sum_{n=0}^{\infty} P(n|\lambda)n\mu = \lambda\mu.$$ Now what about the standard deviation? Let's use $U$ for variables drawn from the unknown distribution, and compute the variance of $T$: $$E((T-E(T))^2) = E(T^2)-E(T)^2 = \sum_{n=0}^\infty P(n|\lambda)nE(U^2) - \lambda^2\mu^2 = \lambda E(U^2) - (\lambda\mu)^2,$$ since the different weights are drawn independently.

This doesn't quite answer our question, since what we're looking for is a way to estimate the standard deviation given a list of $n$ weights $w_i$ drawn from a single sample. But we can take a guess based on the Poisson approach. We can estimate $\lambda$ as $n$ and $E(U^2)$ as the root mean squared of the $w_i$. How does this stand up?

The second of these values is the false positive probability, that is, the chance that you would obtain such a deviation from uniformity drawing from a genuinely uniform distribution.

In short, if you observe $n$ weights $w_i$, the ($1\sigma$) uncertainty you should quote is

$$\sqrt{\sum_{i=1}^n w_i^2 - (\bar w_i)^2}.$$ Note that in spite of the similarity, this is not the formula for the standard deviation of the $w_i$. On the other hand, the second term is generally considerably smaller than the first, by a factor of roughly $n$; omitting it does very little harm, so that it is probably as good to simply add all the weights in quadrature.

A little experimentation shows that for small numbers of photons, the approximation isn't very good. Worse, what matters is the number of photons with weight near the maximum weight. All those small-weight photons don't help improve the approximation substantially.