Obtain complete correlation matrix

[peroju]

Hi gammapy team! I'm Judit (the DM girl) and we are finally finishing the DM analysis pipeline we started some years ago. For this, I wanted also to include some features fo which I encountered some problems.

To set up the discussion, I am using gammapy v-0.18.2 and I am performing a 3D simulation for CTA, including 5 different emission models (DM, CR, 2 AGNs + IRF-BKG), including emission templates for the DM and CR, and modelling the AGNs as point-like sources. This part seems to work perfectly, but I have some problems in the fitting. For the fit, I'm using a template fitting approach, where I leave free up to 8 parameters of these models (mainly amplitudes and spectral indices), and one of the things I would like to get as an output is the correlation matrix between these fitted parameters. The problem that I found is that, even I can obtain the correlation matrix, it only correlates parameters that belong to the same emission, like for example AGN amplitude & index, when I want also to obtain the correlation with the rest, like for example DM amplitude & CR amplitude. I'm obtaining this matrix after the fit though the models attached to the dataset (dataset.models.covariance.correlation). I'm also pretty convinced that I settled correctly the parameters I want to be free, and indeed I obtain them as a result of the fit, but not for the correlation matrix.

Just for completeness, I think it maybe can be useful to show how I create the simulation and perform de fit.

models = Models([model_simu, 
                 cr_model_simu,
                 ps1_model_simu,
                 ps2_model_simu,
                 bkg_model])

dataset = maker.run(empty, obs)
dataset.models = models
# Poission fluctuate in memory
dataset.fake(random_state=int(time.time()))

# Fix to which DM model we want to fit the data
flux_model_fit = DarkMatterAnnihilationSpectralModel(
    mass=m*u.Unit("GeV"),
    channel=ch,
    jfactor=JFAC,
    z=redshift)
# Set max and min for fit to search
flux_model_fit.scale.max = 1e7
flux_model_fit.scale.min = -1e7
# Attach EBL
flux_fit = flux_model_fit * absorption
model_fit = SkyModel(spatial_model=spatial_model, spectral_model=flux_fit)

# Impose the model in the dataset
models_fit = Models([model_fit,
                     cr_model_simu,
                     ps1_model_simu,
                     ps2_model_simu,
                     bkg_model])
#  Untab the ones you want be frozen
models_fit["dataset-simu-bkg"].spectral_model.norm.frozen = False
models_fit["dataset-simu-bkg"].spectral_model.tilt.frozen = False
models_fit["cr-model"].spectral_model.parameters['norm'].frozen = False
models_fit["ps-1"].spectral_model.parameters['amplitude'].frozen = False
models_fit["ps-1"].spectral_model.parameters['index'].frozen = False
models_fit["ps-1"].spatial_model.lon_0.frozen = True
models_fit["ps-1"].spatial_model.lat_0.frozen = True
models_fit["ps-2"].spectral_model.parameters['amplitude'].frozen = False
models_fit["ps-2"].spectral_model.parameters['index'].frozen = False
models_fit["ps-2"].spatial_model.lon_0.frozen = True
models_fit["ps-2"].spatial_model.lat_0.frozen = True

dataset.models = models_fit

#  Instantiate the fitter
fit = Fit([dataset])
fit_result = fit.run()

# Save correlation matrix
np.savetxt(output+'correlation_trial_1.txt'.format(ch), dataset.models.covariance.correlation)

I can also provide the definitions of the models. Thanks a lot for your help in advance! :)

[adonath]

Thanks @peroju! Can you maybe attach a plot of the final correlation matrix, e.g. obtained with dataset.models.covariance.plot_correlation()?

[peroju]

Hi @adonath , thanks a lot for the quick response! Here I provide an example of a correlation matrix obtained with the method you proposed

As there are a lot of parameters, let me guide you through them:

• 3 first correspond to the DM model, only scale (which is an amplitude) is free in the fit
• next 3 correspond to CR model, which only norm is free
• next 5 correspond to the first AGN, where index & amplitude are the only ones free
• next 5 correspond to the second AGN, where index & amplitude are the only ones free
• last 3 correspond to the IRFs background, where norm & tilt are the only ones free

Please let me know any other thing you may need! Thanks a lot! :)

[registerrier]

The AGN amplitude and index seem perfectly anti-correlated. Could there be some parameter degeneracy? Does the fit actually converge? What does fit_result contain?

[peroju]

Hi @registerrier , indeed I've just check for this same example and fit_result.success = False and fit_result.message = 'Optimization failed. Estimated distance to minimum too large. Usually I do not only use a value of one fit, I average the results of O(100) fits to different "simulations" (different dataset.fake(random_state=int(time.time())) and change the DM model to different masses and channels. Indeed, we expect some degeneracy between these parameters, but I thought that the correlation matrix should show this. Still, when I compute the 95% upper limits on the DM scale parameter through the projected loglikelihood, the values that we obtain make perfectly sense, and the loglikelihood behaves nicely as well (attach image). The residuals also seem ok (attach image too) and in the results there are some errors attached to the fitted values (attach as a .txt file).

results_0_b_50000.txt

[AtreyeeS]

Hi @peroju ! Sorry for the late reply! We realised that this was an issue in older versions of gammapy (#4149). Please try with gammapy 1.0 and using the models from FitResult as explained here https://docs.gammapy.org/1.0/tutorials/api/fitting.html#covariance-and-parameters-errors

Thanks!