# Fitting a Gaussian to a histogram when the bin size is significant

I'd like to fit a Gaussian to some experimental data that is binned (the binning is a result of the physical limits of the device). Importantly, the bin size is significant enough that the gaussian cannot be considered flat in the bin window (see pic below). The data is actually 3D but let's just consider the 1D example to start . How does one write a likelihood function for the goodness-of-fit?

My intuition is to simply consider each bin independent and compare the density verses the integrated Gaussian density in the bin window: \begin{align} p(D|\Theta) &= \prod_i^N p(d_i|\Theta) \\ &= \prod_i^N f\left(d_i - \int_{x_i}^{x_{i+1}}\phi(x|\mu,\sigma)dx\right) \end{align} Where N is the number of bins, $d_i$ is the bin height for bin $i$, $\phi(x|\mu,\sigma)$ is the Gaussian PDF, and the integral is over the bin width. My question is: what should I use for $f$? In other words, how is the agreement between $d_i$ and $\phi$ distributed?

• How does this likelihood function change for higher dimensions?
• The integration of a Gaussian over a finite bin size is pretty expensive to compute. Since again my problem is 3D, I'm going to have to do numerical integration MANY times for millions of bins. Is there a faster way to do it?

• Although you discuss goodness of fit, it sounds like what you really want to do is the fitting itself--that is, to estimate the 3D means and covariance matrices. Would this in fact be the case?
– whuber
Commented Nov 14, 2014 at 3:40
• Indeed the goal is to do the fitting. The reason I'm asking for help with the likelihood function is that my eventual goal is to incorporate priors, so everything needs to be probabilistic. If that makes any sense :) Commented Nov 14, 2014 at 3:44
• Sure, it makes sense. I also think things might not be as bleak as you make out: you likely can get an excellent initial fit by sampling each bin at its center (and computing the density there), and then improve on that fit if you like by slightly intensifying that sampling, such as moving to $8$ points per bin, etc. This actually promises to be faster than if you had all the unbinned data.
– whuber
Commented Nov 14, 2014 at 3:47
• @Neil Yes, I'm sure. Remember, I'm talking about using an approximation to achieve an initial fit quickly and inexpensively and then "polishing" that fit based on the full and correct likelihood. Worked examples of how that polishing step would go are presented at stats.stackexchange.com/a/12491 and stats.stackexchange.com/a/68238.
– whuber
Commented Nov 14, 2014 at 16:24
• @Neil That is correct--but in finding the initial ML approximation it is an unnecessary complication, especially in more than one dimension (where the analog of Sheppard's corrections can be expected to be appropriate only for uncorrelated Gaussians). The idea here is to sneak up on a correct ML solution in two steps: the first step essentially ignores the binning, pretending all data are located either in the bin centers, randomly within their bins, or equally spaced within them; after initial estimates of the parameters are obtained, then a better calculation is performed.
– whuber
Commented Nov 14, 2014 at 17:43

If you know that $y_i \in [x_j, x_{j+1})$, where $x_j$'s are cut points from a bin, then you can treat this as interval censored data. In other words, for your case, you can define your likelihood function as

$\displaystyle \prod_{i = 1}^n (\Phi(r_i|\mu, \sigma) - \Phi(l_i|\mu, \sigma) )$

Where $l_i$ and $r_i$ are the upper and lower limits of the bin which the exact value lines in.

A note is that the log likelihood is not strictly concave for many of the models for interval censored data, but in practice this is not of much consequence.

You should treat each bin as if it were generating random points uniformly within its bounds. Therefore calculate a weighted average for each bin $(x_l, x_h]$ of $E(x) = \frac{x_h + x_l}{2}$ and $E(x^2) = \frac{x_h^2 + x_lx_h + x_l^2}{3}$. This weighted average determines a Gaussian.

You can incorporate a prior by treating this Gaussian as a likelihood.

• That's a little confusing, can you elaborate? I understand the idea of each bin generating uniform points, but what do you mean they determine a Gaussian? Commented Nov 14, 2014 at 5:00
• The mean and second moment are one way to parametrize a Gaussian. You could have done mean and variance if you like. Commented Nov 14, 2014 at 5:03
• Ah right of course. Then how do you compare these bin gaussians to the fitting function? Commented Nov 14, 2014 at 5:05
• It is inconsistent to assume the distribution is Gaussian and is uniform within each bin. That inconsistency is going to cause more and more trouble as the dimension increases.
– whuber
Commented Mar 1, 2015 at 21:37
• @whuber: you're right. I guess this could be a starting point of an iterative approach. Commented Mar 1, 2015 at 23:36

Given whuber's comment on my last answer, I suggest you use that answer to find a mean and variance $\mu, \sigma^2$ as a starting point. Then, calculate the log-likelihood of having observed the bin counts you got $\ell$. Finally, optimize the mean and variance by gradient descent. It should be easy to calculate the gradients of the log-likelihood with respect the parameters. This log-likelihood seems to me to be convex.