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?

1. 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.

• Great question. About Poisson noise, my guess is that's about approximating a binomial distribution for large N as a Poisson distribution. About measurement uncertainty, how large is it in relation to the bin size? If measurement uncertainty is an appreciable fraction of the bin size, it might be necessary to bring it into the analysis. Finally, the business about p_k = N_k/N is just the simplest way to estimate the proportion falling into bin k. If you have prior information that p_k must be greater than zero, then you can formalize that and bring it into play. I'll try to write more later. Aug 27 '20 at 18:01
• @RobertDodier You can assume the uncertainty is large enough that each data point is likely to fit in more than one bin within measurement uncertainty. Of course, I can try to resize / reposition the bin or even use instead a KDE to avoid this mess, but before I move to more sophisticated method, I want to make sure I understand these basic problems with histogram or with density estimation in general. Aug 28 '20 at 6:31

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.

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".

• Hi, I'm new here. I know in other Stack Exchange sites, linking is frowned upon, and it is encouraged to copy out the relevant portions of the resource and place them here. However, the amount of relevant information in the particular link is very large, and I feel it's best presented in its current form, at least partially because I am not enough of an expert to be certain I wouldn't be significantly misrepresenting it. Here is an arxiv link of the same paper, just in case the original link goes away. arxiv.org/abs/2101.04783 Nov 11 '21 at 19:49
• Similarly, here is an alternative link for the IFGT paper: academia.edu/2791237/… Nov 11 '21 at 20:04