I work on data from a mass spectrometer that produces billions upon billions of count histograms, and I need a good way to test whether these histograms are consistent with one or several model distributions (Gaussian, heavy-tailed, multimodal, etc). Outliers may be present in a good fraction of the histograms, if not all. The histograms may have anywhere from 0 to 10^6 counts in them, and they come to us already discretized, so the histogram is not losing any info with respect to the original observations.
As a naive jack-of-all-trades data analyst trained by physicists, my instinct is to do something like the following:
For each model distribution,
- estimate its parameters using moments or nonlinear fitting using the Poisson likelihood (since this is count data, each bin is a Poisson random variate)
- calculate the $\chi^2$ of the data vs. the fitted distribution
Then, with the chi squared values from the several models in hand...
- pick the model with the best $\chi^2$ value
- if $\chi^2$ is too big (as referenced against the theoretical $\chi^2$ distribution with the appropriate degrees of freedom), flag the distribution as deviating significantly from the model.
I was curious if more experienced statisticians could advise me on whether this procedure makes sense, limitations I might encounter, better alternatives, etc. Here are a couple things I'm wondering about:
- For histograms with few counts, I feel like it more sense to use the Poisson likelihood / Kullback-Leibler divergence in the goodness-of-fit metric rather than the sum of squares used in the $\chi^2$ test statistic. It's most appropriate to use it in the fitting, why not also in the test? But I don't know any commonly-used test that works this way. I googled around for Poisson histogram goodness-of-fit tests and found nothing.
- I have the vague sense that I should use some AIC type thing to account for the number of parameters in the distribution, but maybe that's already rolled into the $\chi^2$ degrees of freedom.