PIG 26 - Model Priors API#
Author: Noah Biederbeck, Katrin Streil
Created:
2023-03-13Accepted:
2023-11-03Status: Accepted
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:
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 positive 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.
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.
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).
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:#
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
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:
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:
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 Datasets.stat_sum() function is evaluated. Now the function has to compute the statistics due to the priors set on the parameters and add to the stat_sum.
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.
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.
{'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:
Decision#
The PIG was discussed and reviewed among Gammapy developers. It has been reviewed by the CC and is now considered accepted.