Documenting/improving the manipulation of covariance with CompoundSpectralModel

[registerrier]

A relatively common use case is plotting flux points and butterfly plot corrected for the EBL absorption (i.e. plotting the de-absorbed model). To properly work, this requires that the covariance is properly handled. We probably need to cover this use case somewhere in the documentation and there might be some issues in the code.

First, currently covariance handling is documented here:

• https://docs.gammapy.org/1.2/tutorials/api/fitting.html?highlight=covariance#covariance-and-parameters-errors
• the subtleties of covariance propagation are not so clear see e.g. this issue Covariance matrix is not properly set on Datasets #4149
  The documentation might be improved here.

Second, that contrary to what the documentation seems to imply, the covariance matrix is not fully propagated from the CompoundSpectralModel to its submodels model1 and model2. Only the diagonal terms are but not the correlations. This might be a bug (see below an example: opinion @adonath , @QRemy ?). A minima this requires improved documentation on how to handle this specific situation.

Third, we might provide a specific tutorial/recipe/howto on how to plot de-absorbed spectra. The case of a fixed redshift/alpha parameter is not too complex (modulo the covariance issue above). The case of an EBL model with a free parameter is likely more complex to handle and we might not want to cover it for now.

To reproduce:

import numpy as np
from gammapy.modeling.models import PowerLawSpectralModel, EBLAbsorptionNormSpectralModel

pwl = PowerLawSpectralModel()
ebl = EBLAbsorptionNormSpectralModel.read_builtin()
model = pwl * ebl

matrix = np.array(
    [[1.e-02, 4.e-13, 0.0, 0.0, 0.0],
     [4.e-13, 2.e-23, 0.0, 0.0, 0.0],
     [0.0,    0.0,    0.0, 0.0, 0.0],
     [0.0,    0.0,    0.0, 0.0, 0.0],
     [0.0,    0.0,    0.0, 0.0, 0.0]
    ])
model.covariance = matrix     

model.covariance.plot_correlation()

model.model1.covariance.plot_correlation()

[adonath]

Yes, I think this is a bug. For models with sub-model components we need specific handling of setting the covariance See for example: https://github.com/gammapy/gammapy/blob/main/gammapy/modeling/models/cube.py#L139

[registerrier]

@chaimain Would you be interested in making an howto/tutorial on extracting de-absorbed spectra/flux points with gammapy?

[chaimain]

@chaimain Would you be interested in making an howto/tutorial on extracting de-absorbed spectra/flux points with gammapy?

I can help with this, but currently, I am occupied with my PhD defense preparations.