Parameter Inference when Model is a bad fit to the data.

Question

I am working with gamma-ray data from the Fermi Satellite. The data has been binned into energy dependent maps of the sky -- e.g. three dimensions (energy, latitude, longitude) and is extremely high resolution. I then generate very complex models of the gamma-ray emission in the Galaxy, and compare them to the data. There are a fairly small number of important model parameters (tens), while we have something like 10^7 total pixels with counts following a Poisson distribution.

Given this number of pixels, and the relative simplicity of the model, the model does not even come close to describing the sky at the level of Poisson noise. Thus when I fit the data and make tiny adjustments to the parameters, the change in the likelihood can be very dramatic and is a highly noisy surface.

I would now like to run parameter inference, but using the delta chi^2 is not a very sensible thing to do given that there are other much larger uncertainties in the problem, which are not easily quantified... Is there anything I can do to make sense of this without when I do not know the level of systematic uncertainty?

Thank you!

Answer 1

If you have 10^7 pixels and I imagine the 3 dimensions are also per pixel, you might consider collapsing your data. You mentioned pixel counts which are Poisson distributed, so this essentially means your non-zero pixel counts are very rare, so there is sparseness in any matrices used. In this case, either use an optimization method that is based on the Poisson model (with e.g. "rare counts"), or possibly collapse the dimensions. If you can't use a Poisson-based optimization method, maybe try dropping cells in the 3-dimensional grid that have zero counts and then model on what is non-zero. (After all, if you don't go the "Poisson route" which means a tremendous amount of zeroes in the data between the occurrence of ones, twos, threes, etc, then you almost need to drop the zeroes since at this point you wouldn't meet the Poisson assumption anyhow).

Initially, my gut feeling would be to say find a Poisson-based method and then go with it.

For $\chi^2$ and binning for a specific model run, however, you might try binning all of your data and calculating an "observed" chi-squared value called $\chi^2_{obs}=\sum_i(O_i - E_i)^2/E_i$. Then, for this model permute the data (randomly reshuffle) $B=10000$ times calculating observed and expected bin counts after each reshuffling. Each $b$th reshuffling will give you a $\chi^2_{(b)}$, and the number of times the $\chi^2_{(b)} \geq \chi^2_{obs}$ divided by 10,000 is the tail probability for a statistical hypothesis test of how significant (rare) the observed data configuration is when compared to random permuted data. The equation for this is:

$P=\frac{\# \{\chi^2_{(b)} \geq \chi^2_{obs}\}}{B}$

where the numerator is the number of times the $\chi^2_{(b)}$ from permuted data (reshuffled, randomly) exceeds the observed $\chi^2_{obs}$ when the raw data weren't shuffled. A value of $P \leq 0.05$ implies that your data configuration is not random -- that is, you found something. [the statistical interpretation of this would be that, if $P \leq 0.05$, the chances of seeing the pixel counts that you observed is less than 5% (noise is 100%)]. So essentially, what you're doing is making random data from your own data, and trying to determine if bin counts for the observed are significantly different from 10,000 realizations of reshuffled data (not new random data, but your same data which is just reshuffled to change bin counts). These kind of tests are called "randomization tests" or "empirical p-value testing."