# Finding quality of fit for two discrete variables with low statistics

I have data from an experiment which I am trying to explain using a model. I do not have an analytic formula for the prediction of the model but instead I got its prediction through a simulation. The results look something like this (this here is fictional data but hopefully you get the idea):

$\begin{array}{c|c|c} & Experiment & Model \\ \hline \mbox{Bin 1} & 5 \pm 2& 3 \\ \mbox{Bin 2}&6 \pm 2&7 \\ \mbox{Bin 3 }&2 \pm 1 & 3 \\ \vdots & \vdots & \vdots \end{array}$

where the different bins are independent. I'd like to quantify how good the model fits to the data. Chi-squared appears to be the perfect measure if it weren't for the low statistics.

I'm also aware of the likelihood method but the resources I found online explain the method in the context of a continuous distribution with unknown parameters. I've also read about the Fischer test but I'm not sure where it applicable here.

How can I quantify the quality of fit?

• Could you explain what you mean by "low statistics"? Suggesting the use of a $\chi^2$ distribution indicates your data are counts--is that really the case? Exactly how were you able to make a prediction via simulation? How did that work?
– whuber
Mar 27, 2015 at 20:00
• The data is indeed counts. The data arises from an experiment at the Large Hadron Collider (LHC) where they have count the number of events satisfies a set of criteria. However, the criteria are quite stringent so they have very few events in each bin (typically <5 per bin). Mar 28, 2015 at 11:15
• For my simulation I ran a Monte Carlo simulator that reproduces the collisions a the LHC given a certain model. It runs thousands of events and then I apply the same criteria used in the experiment and see how many events I get an each bin, essentially repeating the counting experiment only assuming a particular model. Mar 28, 2015 at 11:18

If you're already simulating from a model, you can use still use the chi-squared test. All you have to do is, instead of analytically calculating the null distribution using the (invalid) large-sample approximation, you figure out its distribution by simulating from your model. That is:

1. Calculate the expected count in each row using your simulations.

2. Calculate the $\chi^2$ statistic of your observed data using these expected counts.

3. Simulate a bunch more data sets from your model and calculate the $\chi^2$ statistic for each one. This is known as (one kind of) bootstrapping.

4. Your p-value is the percentile of the observed $\chi^2$ value among the bootstrapped $\chi^2$ statistics. So if the observed $\chi^2$ is larger than 95% of the simulated ones, then you might reject the hypothesis that the observed data was drawn from your model.

Note that, if your model has high variance and you have low counts, this test may not have very much power, that is, it may not reject even extreme data.

• And this is very easy to do in R. With the example data: chisq.test(x=c(5,6,2), p=c(3,7,3), rescale.p=TRUE, simulate.p.value=TRUE, B=10^4). You can change the number of simulations from the null distribution by changing the value of B (here set to 10,000 simulations). Mar 29, 2015 at 18:21
• Thanks for your response. What I find confusing about this method is it doesn't seem to take into account the error bars of the experiment (i.e., it seems that you would get the same value regardless of these errors)? Is there an assumption made about the size of the systematic errors or am I missing something? Mar 30, 2015 at 10:57
• @JeffDror: sorry--by "model" in the above, I meant "model of the data-generating process", not "model of the physical process". Uncertainty about measurements should be included in your model of the data-generating process. In other words, when you simulate from the model, you should simulate a set of true values from the model, and then simulate a set of measured values from the true values plus your model of the measurement uncertainty. Mar 30, 2015 at 22:13