PIG 16 - Model Priors API

[nbiederbeck]

This pull request adds a rough PIG for Prior (and FitStatistic). I anticipate I will add details with your assistance.

[codecov]

Codecov Report

Merging #4381 (f6d08d2) into main (e6e5dd8) will decrease coverage by 0.28%.
Report is 165 commits behind head on main.
The diff coverage is n/a.

@@            Coverage Diff             @@
##             main    #4381      +/-   ##
==========================================
- Coverage   75.90%   75.62%   -0.28%     
==========================================
  Files         223      226       +3     
  Lines       32312    33026     +714     
==========================================
+ Hits        24527    24977     +450     
- Misses       7785     8049     +264     

see 61 files with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[bkhelifi]

Suggested change

	Using priors on models and on parameters is the next step in full-fledged statistical analyses.
	Using priors on models or on parameters is the next step in full-fledged statistical analyses.

[bkhelifi]

Maybe not only the analyser... Indeed, one idea to deal with the IRFs systematics was to use the priors mechanics. So, if it is implemented like that, the knowledge comes from the people making the IRFs...

By the use, do we add as use case the inclusion of IRFs systematics into priors?

[adonath]

Thanks a lot @nbiederbeck for this minimal draft. The overall structure is good and the content should be extended for details now:

• We need the definition of an abstract Prior base class. Define the evaluation and parameter API.
• Clarify whether priors can have fittable parameters (yes!)
• Some priors do not have fitable parameters...
• Do we need a special PriorParameter class?
• Clarify serialisation of priors with the model e.g. via an example YAML file
• Mention prior registry system
• Add a list of proposed standard priors: GaussianPrior, UniformPrior, MultivariateGaussianPrior
• Strictly speaking one would not need the L2Regularization as it presents a multivariate Gaussian prior.
• Clarify which Gammapy classes will support a .prior attribute: Models, SpectralModel, SkyModel etc.
• Clarify how (hierarchical) priors are combined in the fit-statistic function
• Do we need e.g. a class like PriorFitStatistic(beta=1.)?

These are just some thoughts for now, I'm sure I'll have some more...

[adonath]

Note that this case actually corresponds to the GaussianPrior on each of the individual parameters, where mu=0 and sigma=1. The most general case I could imagine would be MultiVariateGaussianPrior on a list of parameters. The case of L2 regularization would be included there...but probably we just keep it for convenience.

[adonath]

@katrinstreil kindly agreed to take over the writing and further development of the PIG. Please let us of know if you have any questions.

[registerrier]

Thanks @katrinstreil @nbiederbeck . This already quite nice.

I would like to clarify a number of potential issues:

1. it seems a bit complicated to me that priors can live and be defined and evaluated at different levels in the modeling package structure. I don't think that the SkyModel is the correct place anyway. We could have something to set priors on the Models instance with e.g. models.set_prior(parameters, prior).
2. Parameters have units. How do we make sure that units are correctly handled?
3. Is the PriorFitStatistic really needed? As far as I can see, it is mostly a single function.
4. Connected to this, I wonder if the mixing of two use cases in a single prior framework is not a bit confusing. We have:

• priors on individual parameters as usually done in e.g. sherpa, bxa, 3ML.
• penalization terms for the ML which are used for e.g. regularization in unfolding
  Maybe this should be distinct features?

[registerrier]

Does this mean that Prior has to contain Parameter objects? Doesn't this cause problems?

[adonath]

Yes, I think in general a prior should be treated much like a model. In general prior parameter can be fitted in a likelihood fit (maximum a posteriori) as well...so I think they should handle Parameter objects.

[registerrier]

Should this follow the general scheme applied in models where evaluate is a static method that takes all parameters on input?

[adonath]

Yes, I agree that the .evaluate() method should be static if possible. This will allows us to easier move to other tensor libraries later...The stateful call can happen in Prior.__call__

[registerrier]

Does it really make sense to define the prior on the SkyModel? I would assume the more logical place is Parameters.
I suppose there are use cases where the prior applies to parameters of distinct SkyModel.
You could still do model.parameters.prior = some_prior

[adonath]

The reason why it should be attached to the Model is bookkeeping. Priors should be serialized with the model. Parameters do not know, to which model they belong. The alternative is to handle it completely independent as we do for the covariance matrix and we would rely on the order of parameters. But in this case it would still be attached to the Model.

[registerrier]

How does it deal with the Parameter unit?

[adonath]

I guess priors should be defined in such a way, that they are unit less (i.e. a likelihood after evaluation). But maybe we should add a check for this...

[registerrier]

The need of an additional class seems unclear at the moment. It does only stat_sum, you could equivalently do something like this in Datasets.stat_sum():

    def stat_sum(self):
        """Compute joint statistic function value."""
        stat_sum = 0
  
        for dataset in self:
            stat_sum += dataset.stat_sum()
        prior_stat_sum = self.models.parameters.prior_stat_sum()
        return stat_sum + prior_stat_sum

[adonath]

The idea of the additional class is just to make it explicit, that one changes to a MAP estimate. It's like a global switch to account for the prior or not. The alternative is indeed to always evaluate priors, if present...

[adonath]

Connected to this, I wonder if the mixing of two use cases in a single prior framework is not a bit confusing. We have:
priors on individual parameters as usually done in e.g. sherpa, bxa, 3ML.
penalization terms for the ML which are used for e.g. regularization in unfolding
Maybe this should be distinct features?

I don't think they are not distinct features. If you do the math you will find the L2 regularization corresponds to a multivariate Gaussian prior in a Maximum A-Posteriori estimation. The information is all in the covariance matrix. I think the right approach here is to just introduce a MultivariateGaussianPrior and some convenience methods for the creation of covariance matrices. See my sample above.

[registerrier]

Thanks @katrinstreil .

Let's now send it to the CC.

[registerrier]

@katrinstreil , you will also need to update the index.rst file to add the PIG to the list.

[bkhelifi]

Thank you, @katrinstreil and @nbiederbeck

[registerrier]

Review period is over. The PIG is accepted. Merging now.