From full to simplified statistical models

[jackaraz]

Hi all, I'm trying to find a "proper" way to use only a subset of regions used in a channel(s). As @alexander-held graciously explained in here it is not very kosher to break a full statistical model apart since it can mess up the original fit. So I was wondering if it would be possible to extract a covariance matrix from pyhf to be used in simplified likelihood construction. I know that there is no such functionality at the moment, but if you can give me some pointers to write such a module that extracts a covariance matrix, possibly for a skewed gaussian distribution, to preserve as much information as possible, that would be great. Thanks a lot!

[alexander-held]

Hi @jackaraz, you can get the covariance matrix when using MINUIT as optimizer. Here is an example:

import pyhf

pyhf.set_backend("numpy", "minuit")

model = pyhf.simplemodels.correlated_background(
    signal=[12.0, 11.0],
    bkg=[50.0, 52.0],
    bkg_up=[51.0, 53.0],
    bkg_down=[49.0, 52.0],
)

res, res_obj = pyhf.infer.mle.fit([63, 64] + model.config.auxdata, model, return_result_obj=True)
print(res_obj.minuit.covariance)
# ┌────────────────────────────┬───────────────────────────────────────────────────────┐
# │                            │ correlated_bkg_uncertainty                         mu │
# ├────────────────────────────┼───────────────────────────────────────────────────────┤
# │ correlated_bkg_uncertainty │                          1                    -0.0662 │
# │                         mu │                    -0.0662                      0.244 │
# └────────────────────────────┴───────────────────────────────────────────────────────┘

The object you get out is https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.covariance.

[matthewfeickert]

@jackaraz yeah, this is what I was recommending yesterday in our informal coffee chat. Let us know if you have any additional questions here, or if you wanted something else.

[jackaraz]

Thank you both! I’ll try this out and get back to you!

[jackaraz]

Hi @alexander-held and @matthewfeickert, so as far as I understand this cov matrix is basically $-E[\partial^2 L / \partial\theta_i\partial\theta_j]$ correct? So the terms refer to the nuisance parameters. Is there a way to extract covariance matrix for regions only? i.e. I want to capture the correlation between regions for bkg only sample.

[alexander-held]

Yes, this is the covariance for the nuisance parameters. It's not obvious to me what you mean exactly by "covariance matrix for regions only". Perhaps it is similar to what we need for calculating post-fit yield uncertainties (per bin, but there is also a way to do this per region that is less meaningful as the fit does not consider "regions"). This is implemented in cabinetry (cabinetry.model_utils.yield_stdev, though the internals are not exposed) so that could perhaps be a starting point.