PIG 16 - Model Priors API

docs/development/pigs/pig-026.rst

@@ -0,0 +1,287 @@
+.. include:: ../../references.txt
+
+.. _pig-026:
+
+********************************
+PIG 26 - Model Priors API 
+********************************
+
+* Author: Noah Biederbeck, Katrin Streil
+* Created: ``2023-03-13``
+* Accepted:
+* Status: draft
+* Discussion: `#4381`_
+
+Abstract
+========
+
+This PIG is intended to introduce priors on parameters that are evaluated during fitting.
+
+Motivation
+==========
+
+Using priors on models or on parameters is the next step in full-fledged statistical analyses. 
+A prior can incorporate the analysers' knowledge of the topic at hand or information about
+the estimated IRFs systematics provided by the corresponding experiment, and yield more realistic and trustworthy results. 
+The proposed formalism also includes the application of nuisance parameters and regularization.
+
+Use cases
+=========
+
+In the past priors were a regularly requested feature or the solution to a problem.
+See the following issues and PRs of Gammapy:
+
+Case 1: Background systematics as a nuisance parameter `#3955`_
+--------------
+The goal is to fit energy-dependent systematics in the FoVBackground with nuisance parameters and set priors on the two model parameters ``norm`` and ``tilt``:
+
+
+.. code:: python 
+
+  from gammapy.modeling.model import GaussianPrior, PowerLawNormSpectralModel
+
+  bkg_model = FoVBackgroundModel(spectral_model = PowerLawNormSpectralModel(),
+                          dataset_name = dataset.name)
+  tilt_prior =  GaussianPrior(mu = "0", sigma = "0.05")
+  norm_prior =  GaussianPrior(mu = "1", sigma = "0.1")
+
+  bkg_model.set_prior([bkg_model.parameters['tilt'],bkg_model.parameters['norm']], [tilt_prior, norm_prior])
+
+
+A preliminary version was set up like this and tests on simulated datasets were successful. Note that this setup is quite simple and it can be developed more advanced based on the same principle. 
+
+Case 2: Favoring postivive values for flux amplitudes
+--------------
+
+A step-like prior function can be used to favour positive values for physical properties like the flux amplitude. By setting a prior one avoids defining hard boundary conditions with the ``min`` attribute of the to-be-fitted parameter. 
+The prior is set to ``value`` if the parameter value is between ``xmin`` and ``xmax`` and 0 if not.
+
+.. code:: python 
+
+   from gammapy.modeling.models import PowerLawSpectraModel, StepPrior
+
+   pwl = PowerLawSpectraModel()
+   prior = StepPrior(xmin = "-inf", xmax = "0",  value = "1")
+   pwl.set_prior([pwl.parameters['amplitude']], [prior])
+
+
+
+
+Case 3: Support unfolding methods for spectral flux points `#4122`_
+--------------
+
+The proposed prior class will allow the probability to unfold spectral flux points with Tikonov regularisation. The Tikonov matrix can be defined and set as the covariance matrix in the class ``CovarianceGaussianPrior``. The weight of the ``PriorFitStatistic`` can be interpreted as the regularisation strength tau. However, this requires a more advanced setup since this multi-dimensional prior is set on multiple parameters simultaneously. The goal is to use the proposed prior class as a starting point and develop it into multidimensional use cases in the near future. 
+A suggested implementation example is shown below. Here the prior is set on the model instead of the single parameters.
+
+.. code:: python 
+
+  import numpy as np
+  from astropy import units as u
+  from gammapy.modeling.model import PiecewiseNormSpectralModel, MultivariateGaussianPrior, PowerLawSpectralModel
+
+  n_points = 10
+
+  energy = np.geomspace(1, 10, n_points) * u.TeV
+  norm =  PiecewiseNormSpectralModel(energy=energy)
+
+  prior = MultivariateGaussianPrior.from_covariance_matrix(means=np.ones(10), covariance_type="diagonal")
+  prior.weight = 1.1 
+  norm.prior = prior
+
+  spectral_model = PowerLawSpectraModel() * norm
+
+Additional use case:  Add sigma v estimator functionality for dark matter use case `#2075`_
+
+Implementation
+==============
+
+The following implementation draft is compatible with the current API,
+where the models are set via Dataset.model and the Fit is run via Fit.run(datasets).
+
+.. code:: python
+
+  class PriorModel(ModelBase):
+
+        _weight = 1
+        _type = "prior"
+
+        @property
+        def parameters(self):
+            """PriorParameters (`~gammapy.modeling.PriorParameters`)"""
+            return PriorParameters(
+                [getattr(self, name) for name in self.default_parameters.names]
+            )
+
+
+        @property
+        def weight(self):
+            return self._weight
+
+        @weight.setter
+        def weight(self, value):
+            self._weight = value
+
+        def __call__(self, value):
+            """Call evaluate method"""
+            kwargs = {par.name: par.quantity for par in self.parameters}
+            if isinstance(value, Parameter):
+                return self.evaluate(value.quantity, **kwargs)     
+             else:
+                raise TypeError(f"Invalid type: {value}, {type(value)}")
+
+
+The base ``PriorModel`` class inherits from the ``ModelBase`` class and allows to set (for now one) unit to which the to-be-evaluated parameter the prior is set on is converged. In addition, there can be a weight set. The parameters of the  ``PriorModel`` are ``PriorParameters`` and ``PriorParameter``. They inherit from the ``Parameters`` and ``Parameter`` classes, respectively, however, have a limited amount of attributes. 
+It is assumed that the ``PriorParameters`` are set in the same unit as the ``Parameter`` they get evaluated on.  
+
+Future additional properties of the ``PriorModel`` classes: 
+
+- Write/read from/to a yaml file, ideally also when the corresponding model is written/read (see serialisation example below)
+- Prior registry system
+- Different prior subclasses depending on the use cases 
+
+Exemplary additional prior subclasses:
+---------------------------------------
+.. code:: python
+
+  class GaussianPrior(PriorModel):
+
+        """Gaussian Prior with mu and sigma.
+        """
+        tag = ["GaussianPrior"]
+        _type = "prior"
+        mu = PriorParameter(name="mu", value = 0, unit = '')
+        sigma = PriorParameter(name="sigma", value = 1, unit = '')
+
+        @staticmethod
+        def evaluate(value, mu, sigma):
+            return ((value - mu) / sigma) ** 2
+
+
+.. code:: python
+
+  class UniformPrior(PriorModel):
+
+        """Uniform Prior
+        """
+        tag = ["UniformPrior"]
+        uni = PriorParameter(name="uni", value = 0, min = 0, max = 10, unit = '' )
+
+        @staticmethod
+        def evaluate(value, uni):
+            return uni
+
+
+The priors can be set on a ```Parameter``` as ``.prior`` and get evaluated within:
+
+.. code:: python 
+
+  class Parameter():
+
+        _prior = None
+
+        @property
+        def prior(self):
+            return self._prior
+
+        @prior.setter
+        def prior(self, value):
+            self._prior = value
+
+        def prior_stat_sum(self):
+            if self.prior is not None:
+                return self.prior.weight * self.prior(self)
+
+The ``Parameters`` inherit the priors and evaluate the ``prior_stat_sum`` of all the ``Parameters.parameters``:
+
+.. code:: python 
+
+  class Parameters():
+
+        @property
+        def prior(self):
+              return [par.prior for par in self]
+
+          def prior_stat_sum(self):
+              parameters_stat_sum = 0
+              for par in self:
+                  if par.prior is not None:
+                      parameters_stat_sum += par.prior_stat_sum()
+              return parameters_stat_sum
+
+
+During the Fit the ``stat_sum`` of the Datasets is evaluated where now the statistics due to the priors set on the parameters is added. 
+
+.. code:: python 
+
+  class Datasets():
+
+        def stat_sum(self):
+            """Total statistic given the current model parameters."""
+            stat = self.stat_array()
+
+            if self.mask is not None:
+                stat = stat[self.mask.data]
+            prior_stat_sum = self.models.parameters.prior_stat_sum()
+            return np.sum(stat, dtype=np.float64) + prior_stat_sum
+
+
+A simple method in the ``Models`` class will allow setting priors on multiple parameters simultaneously. 
+
+.. code:: python 
+
+  class Models(Models):
+
+        def set_prior(self, parameters, priors):
+            for parameter, prior in zip(parameters, priors):
+                parameter.prior = prior
+
+
+
+
+Serialisation 
+------------
+
+The priors are serialised and read and writeable to a yaml file. They are also saved if the associated ``Parameter`` or ``Parameters`` is transformed to a dictionary. For this, the prior subclasses are tagged. The following two example shows the resulting yaml file for a ``Parameter``.
+
+.. code-block:: yaml
+
+
+  {'name': 'testpar',
+   'value': 0.1,
+   'unit': '',
+   'error': 0,
+   'min': nan,
+   'max': nan,
+   'frozen': False,
+   'interp': 'lin',
+   'scale_method': 'scale10',
+   'is_norm': False,
+   'prior': {'tag': 'GaussianPrior',
+            'parameters': [{'name': 'mu', 'value': 0},
+                           {'name': 'sigma', 'value': 0.1}]}}
+
+
+Implementation Outline
+--------------
+
+The PIG implementation will be piece-wise within individual issues in the following order:
+
+1. ``Prior`` base class and subclasses (See `#4620`_)
+2. ``PriorParameter`` and ``PriorParameters`` (See `#4620`_)
+3. ``.prior`` attribute in the ``Parameter``, ``Parameters`` and ``Model`` classes
+4.  Evaluation of the prior in the ``Datasets`` class
+
+
+
+Decision
+========
+
+
+
+
+.. _#2075: https://github.com/gammapy/gammapy/issues/2075
+.. _#3955: https://github.com/gammapy/gammapy/issues/3955
+.. _#4122: https://github.com/gammapy/gammapy/issues/4122
+.. _#4208: https://github.com/gammapy/gammapy/issues/4208
+.. _#4620: https://github.com/gammapy/gammapy/pull/4620
+.. _#4381: https://github.com/gammapy/gammapy/issues/4381