Error Bars for Histogram with Uncertain Data

Question

Context

I have a set of data points $\{x_1, \dots, x_N \}$ along with the respective measurement uncertainties $\{\epsilon_1, \dots, \epsilon_N\}$ in them ($N \approx 100$). These data are basically the measured distances to the occurrences of some astrophysical process, and I am trying to estimate the spatial distribution of these events without assuming any model (because I really don't have a reasonable model). So to do that, I built a histogram out of my data with bins of equal size $\{B_0, \dots, B_M\}$, and now I want to also put some error bars on my histogram, with my measurement uncertainties taken into account. But after I have looked around for how to do this, I got even more confused.

(I don't have much experience with statistics, so the real problems may just be my lack of understanding in statistics.)

Histogram with no measurement uncertainty

First of all, I found that I don't seems to even understand what these error bars suppose to mean. Let's first ignore the $\epsilon_i$'s and compute the error of a histogram of "perfect data". I have come across the following calculation in several different places:

Denote the number of data points fall in bin $B_k$ correspondingly as $N_k$. We estimate the probability of fall in this bin as $p_k = \frac{N_k}{N}$. Then since we can thought of $N_k$ as a sum of Bernoulli variable $Ber(p_k)$, the variance of $N_k$ is just $\sigma^2[N_k] = Np_k(1-p_k) = N_k(1-\frac{N_k}{N})$. For large enough $N$, we can ignore the second term, and we have the error bar $\sigma_k = \sqrt{N_k}$.

But I don't understand:

1. I saw people often refer this as a "Poisson noise", but I am not sure if I see where is that underlying Poisson process generating this Poisson noise.

2. This also suggest that bins with zero count has no error, which doesn't sound right to me. Indeed, I have come across this article discussing exactly what's wrong with assigning a Poisson error bar $\sigma_k = \sqrt{N_k}$. In particular, the author says

If we observe N, that measurement has NO uncertainty: that is what we saw, with 100% probability. Instead, we should apply a paradigm shift, and insist that the uncertainty should be drawn around the model curve we want to compare our data points to, and not around the data points!

But that doesn't sound right neither. While my measurements are deterministic numbers (ignoring measurement uncertainty), I am trying to estimate a distribution using a finite sample, so there still got to be uncertainty associated with my estimation. So what should be the correct way to understand these issues?

3. I have also been suggested to use bootstrapping to estimate these error bars, but again I don't quite understand why should it work. If $N_k=0$ for my original data set, no matter how I resample my data, I will always have zero count in $B_k$, so I am again forced to conclude that $p_k = 0$ with zero uncertainty. So intuitively I don't see how bootstrapping my data can give me any new insight about my distribution estimate. Well, it may just be that I don't understand how resampling methods work in general.

Histogram with measurement uncertainty

Coming back to my original problem. I did find some answers about how to put in measurement uncertainties such as in this answer. The method basically is to find the probability $q_i(B_k)$ of the $i$-th data point falling in bin $B_k$ assuming the $i$-th measurement is normal distributed with $\mathcal{N}(x_i, \epsilon_i^2)$:

$$ q_i(B_k) = \int_{B_k} \frac{1}{\sqrt{2\pi}\epsilon_i} e^{-\frac{(x-x_i)^2}{2\epsilon_i^2}} \ dx$$

And then use these $q_i(B_k)$ to construct the Bernoulli variance in $B_k$ as

$$ \sum_{i=1}^{N} q_i(B_k)(1 - q_i(B_k)) $$

But my question is, where does that "Poisson noise" go in this method? The bin count $N_k$ doesn't even show up anymore, and this make me feels like something is missing. Or maybe I have overlooked something.

So I guess what I really want, is to see a complete treatment of error estimation for histogram, which I couldn't find anywhere.

Answer 1

I thought about it some more, and I have a couple of ideas.

(1) About measurement uncertainty: from what you said, it's big enough to take into account. I agree with the formula for qi -- this is just the mass of the distribution for x[i] which falls into B[k]. From that, it looks to me that the mean of the proportion of x which falls into B[k] (let's call that q(B[k])) is the sum of those bits over all the data, i.e., q(B[k]) = sum(qi, i, 1, N). Then the height of the histogram bar k is q(B[k]). and its variance is q(B[k])*(1 - q(B[k])).

So I disagree about the variance -- I think the summation over i should be inside q in variance = q*(1 - q), not outside.

It occurs to me that you'll want to ensure that the q(B[k]) sum to 1 -- maybe that's guaranteed by construction. In any event you'll want to verify that. EDIT: Also, as the measurement error becomes smaller and smaller, you should find that the q(B[k]) converges to the simple n[k]/sum(n[k]) estimate.

(2) About prior information about nonempty bins, I recall that adding a fixed number to the numerator and denominator in n[k]/n, i.e., (n[k] + m[k])/(n + sum(m[k])), is equivalent to assuming a prior over the bin proportion, with the prior mean being m[k]/sum(m[k]). As you can see, the larger m[k], the stronger the influence of the prior. (This business about the prior count is equivalent to assuming a conjugate prior for the bin proportion -- "conjugate prior beta binomial" is a topic you can look up.)

Since q(B[k]) is not just a proportion of counts, it's not immediately clear to me how to incorporate the prior count. Maybe you need (q(B[k]) + m[k])/Z where Z is whatever makes the adjusted proportions sum to 1.

However, I don't know how hard you should try to fix up the bin proportions. You were saying you don't have enough prior information to pick a parametric distribution -- if so, maybe you also don't have enough to make assumptions about bin proportions. That's a kind of higher-level question you can consider.

Good luck and have fun, it seems like an interesting problem.

Answer 2

I have a similar problem and I have a solution in mind, though it's more complicated than I'd like, which is how I stumbled upon this answer: seeking an easier answer.

With that preamble out of the way, I'll share the solution I had in mind. It is more complicated than I'd like. And I am sorry that this is coming over a year after your original post. This is almost certainly no longer relevant.

To begin with, you do not have a model for the underlying probability distribution. There are several ways that you infer one, and the way we'll be discussing is using a technique called kernel density estimation, which is essentially that you put a gaussian (or other kernel function) at each observation point, then sum all of the gaussians up. The variance of each gaussian is related to a parameter called the bandwidth, and there are some algorithms that you may use to come up with a "good" bandwidth, but it's a hyperparameter, and it's mostly guesswork. There's also a method of using one of these guesses as input to another round of the kernel density estimation this time using the frequency of the points in the neighborhood to automatically adjust the bandwidth of each input point. This is called "variable bandwidth" kernel density estimation.

However, your observations also have a significant associated uncertainty. You'll need to convolve each kernel with the gaussian from the measurement uncertainty, and then convolve again during the variable bandwidth stage. (Because convolving a normalized gaussian with another normalized gaussian is as simple as adding variances I recommend using a gaussian kernel. Though the specific kernel you use is irrelevant with enough data points.)

(That being said, I am not certain and I haven't fully thought through the consequences of convolving the measurement uncertainty at both stages vs only one or the other. However, it seems to me that both stages is the correct call, n.b. if we consider each observation to be several hundred observations (without associated uncertainty) distributed according to the original observation's uncertainty, then this is roughly equivalent to the convolution (and exactly equivalent in the limit as the number of samples goes to infinity), and in this concrete case is equivalent to convolving at both stages.)

This should give you a PDF of the resulting distribution. You may then rebin this density, using the total weight of the pdf contained in that bin along with the total count of all observations, to give error bars for each bin. (In this case, you will also have zero error on any bin with zero value; however, because of the use of gaussians, no bin will be exactly zero, and any bin with near zero value for after using this technique is almost certain to actually have near zero weight.)

Variable Bandwidth Kernel Density Estimation

You also mentioned that you only have a few hundred points, I have a few 10s of thousand on the other hand, so I want a fast method for summing gaussians. This is given in the paper "Improved fast Gauss transform with variable source scales".

Answer 3

Check this link about the error bars for cases similar to the one you mention.

https://www.science20.com/quantum_diaries_survivor/those_deceiving_error_bars-85735

We are measuring the number counts within bins. (Similar to building a histogram).

What is usually done? R- The errors in each bin are considered symmetric, with value = \sqrt(N) , where N is the number of events (or counts) in that particular bin. So, the data point N, gets positive and negative errors stretching from (N + \sqrt(N)) down to (N - \sqrt(N) ).

" In other words, the default is to use the fact that the event counts, being a random variable drawn from a Poisson distribution, has a variance equal to the mean. "

When we see this kind of plot, it is said to consider "Poisson errors"... And that is the "typical" history.

But, there is much more to it! From the link https://www.science20.com/quantum_diaries_survivor/those_deceiving_error_bars-85735 . I quote here (for easy of access and legacy):

" Any statistics textbook explains that the Poisson is a discrete distribution describing the probability to observe N counts when an average of m is expected. Its formula, P(N|m)=[exp(-m)* m^N]/N! (where ! is the symbol for the factorial, such that N!=N*(N-1)(N-2)...*1, and P(N|m) should be read as "the probability that I observe N given an expectation value of m").

So what's the problem with Poisson error bars ? The problem is that those error bars are not representing exactly what we would want them to. A "plus-or-minus-one-sigma" error bar is one which should "cover", on average, 68% of the time the true value of the unknown quantity we have measured: 68% is the fraction of area of a Gaussian distribution contained within -1 and +1 sigma. For Gaussian-distributed random variables a 1-sigma error bar is always a sound idea, but for Poisson-distributed data it is not typically so. What's worse, we do not know what the true variance is, because we only have an estimate of it (N), while the variance is equal to the true value (m). In some cases this makes a huge difference.

Take a bin where you observe 9 event counts and the true value was 16: the variance is sqrt(16)=4, so you should assign an error bar of +-4 to your data point at 9. But you do not know the true value, so you correctly estimate it as N=9, whose square root is 3. You thus proceed to plot a point at 9 and draw an error bar of +-3. Upon visually comparing your data point (9+-3) with the expectation from a true model, drawn as a continuous histogram and having value 16 in that bin, you are led to believe you are observing a significant negative departure of event counts from the model, since 9+-3 is over two "sigma" away from 16; 9 is instead less than two sigma away from 16+-4. So the practice of plotting +-sqrt(N) error bars deceives the eye of the user.

Worse still is the fact that the Poisson distribution, for small (m<=50 or so) expected counts, is not really symmetric. This causes the +-sqrt(N) bar to misrepresent the situation very badly for small N. Let us see this with an example.

Imagine you observe N=1 event count in a bin, and you want to draw two models on top of that observation: one model predicts m=0.01 events there, the other predicts m=1.99. Now, regardless of whether m=0.01 or m=1.99 is the expectation of the event counts, if you see 1 event you are going to draw a error bar extending from 0 (i.e., 1-sqrt(1)) to 2 (1+sqrt(1)), thus apparently "covering" the expectation value in both cases; but while for m=1.99 the probability to observe 1 event is very high (and thus the error bar around your N=1 data point should indeed cover 1.99), for m=0.01 the probability to observe 1 event is very small: P(1|0.01)=exp(-0.01)0.01^1/1!=0.01exp(-0.01)=0.0099. N=1 should definitely not belong to a 1-sigma interval if the expectation is 0.01, since almost all the probability is concentrated at N=0 in that case (P(0|0.01)=0.99)!

The solution, of course, is to try and draw error bars that correspond more precisely to the 68% coverage they should be taken to mean. But how to do that ? We simply cannot: as I explained above, we observe N, but we do not know m, so we do not know the variance. Rather, we should realize that the problem is ill-posed. If we observe N, that measurement has NO uncertainty: that is what we saw, with 100% probability. Instead, we should apply a paradigm shift, and insist that the uncertainty should be drawn around the model curve we want to compare our data points to, and not around the data points!

If our model predicts m=16, should we then draw a uncertainty bar, or some kind of shading, around that histogram value, extending from 16-sqrt(16) to 16+sqrt(16), i.e. from 12 to 20 ? That would be almost okay, were it not for the asymmetric nature of the Poisson. Instead, we need to work out some prescription to count the probability of different event counts for any given m (where m, the expectation value of the event counts, is not an integer!), finding an interval around m which contains 68% of it.

Sound prescriptions do exist. One is the "central interval": we start from the value of N which is the nearest integer to m smaller than m, and proceed to move right and left summing the probability of N+1 and N-1 given m, taking in turn the largest of these. We continue to sum until we exceed 68%: this gives us a continuous range of integer values which includes m and "covers" precisely as it should, given the Poisson nature of the distribution.

Another prescription is that of finding the "smallest interval" which contains 68% or more of the total area of the Poisson distribution for a given m. But I do not need to go into that kind of detail. Suffices here to point the interested reader to a preprint recently produced by R.Aggarwal and A.Caldwell, titled "Error Bars for Distributions of Numbers of Events". The paper also deals with the more complicated issue of how to include in the display of model uncertainty the systematics on the model prediction, finding a Bayesian solution to the problem which I consider overkill for the problem at hand. I would be very happy, however, if particle physics experiments turned away from the sqrt(N) error bars and adopted the method of plotting box uncertainties with different colours, as advocated in the cited paper. You would get something like what is shown in the figure on the right.

Note how the data points can now be immediately classified here, and more soundly, according to how much they depart to the model, which is now not anymore a line, but a band giving the extra dimensionality of the problem (the model's probability density function, as colour-coded by green for 68% coverage, yellow for 95% coverage, and red for 99% coverage). I would be willing to get away without the red shading -68% and 95% coverages would suffice, and are more in line with current styles of showing expected values adopted in many recent search results (the so-called "Brazil-bands").

Despite the soundness of the approach advocated in the paper, though, I am willing to bet that it will be very hard to impossible to convince the HEP community to stop plotting Poisson error bars as sqrt(N) and start using central intervals around model expectations. An attempt to do so was made in BaBar, but resulted in a failure -everybody continued to use the standard sqrt(N) prescription. There is a definite amount of serendipity in the behaviour of the average HEP experimentalist, I gather! "

=====

From the paper "Error bars for distributions of number events" https://arxiv.org/pdf/1112.2593.pdf

Firstofall,there is no uncertaint on the number of observed events. We certainly do not mean that there is a high probability that we had 2.3 rather than 2 events in the 7th bin in the plot. Actually, the error bar is intended to represent the uncertainty on a different quantity -the uncertainty on the mean of an assumed underlying Poisson distribution-.

=========

So, the best solution seems to be: i) consider an underlying model

ii) from the model, get the number of expected events 'N_expected' for a given bin, and add errors to it as +- sqrt(N_expected) .

iii) check if [N_observed] is compatible to [N_expected +- sqrt(N_expected)]

iv) Now you could argue that the observation is/isn't compatible with expectations, within x sigma . (where the 1 sigma interval ~ sqrt(N_expected) )

v) To improve: could try to incorporate the symmetry for the positive/negative error bars (as discussed in arXiv:1112.2593).

vi) Take Away: We should get used to adding error bars to the N_expected from modeling, not to N_observed. There should be no Poisson-errors associated to N_observed, since those could be misleading when trying to decide each model is the best one to explain observations (as discussed in arXiv:1112.2593 and other resources)