Prior statistics and set-up for multidimensional priors

gammapy/datasets/core.py

@@ -9,6 +9,7 @@
 from astropy.table import Table, vstack
 from gammapy.data import GTI
 from gammapy.modeling.models import DatasetModels, Models
+from gammapy.stats import prior_fit_statistic
 from gammapy.utils.scripts import make_name, make_path, read_yaml, write_yaml
 
 log = logging.getLogger(__name__)
@@ -75,7 +76,7 @@ def stat_sum(self):
             stat = stat[self.mask.data]
         prior_stat_sum = 0.0
         if self.models is not None:
-            prior_stat_sum = self.models.parameters.prior_stat_sum()
+            prior_stat_sum = prior_fit_statistic(self.models.priors)
         return np.sum(stat, dtype=np.float64) + prior_stat_sum
 
     @abc.abstractmethod
@@ -229,9 +230,16 @@ def stat_sum(self):
         """Compute joint statistic function value."""
         stat_sum = 0
         # TODO: add parallel evaluation of likelihoods
+        prior_stat_sum = 0
         for dataset in self:
             stat_sum += dataset.stat_sum()
-        return stat_sum
+            # remove prior_fit_statistic from individual dataset to avoid double counting
+            if dataset.models is not None:
+                prior_stat_sum -= prior_fit_statistic(dataset.models.priors)
+        # compute prior_fit_statistics from all datasets
+        if self.models is not None:
+            prior_stat_sum += prior_fit_statistic(self.models.priors)
+        return stat_sum + prior_stat_sum
 
     def select_time(self, time_min, time_max, atol="1e-6 s"):
         """Select datasets in a given time interval.

gammapy/datasets/map.py

@@ -16,6 +16,7 @@
     cash,
     cash_sum_cython,
     get_wstat_mu_bkg,
+    prior_fit_statistic,
     wstat,
 )
 from gammapy.utils.deprecation import deprecated_renamed_argument
@@ -220,7 +221,6 @@ def __init__(
     ):
         self._name = make_name(name)
         self._evaluators = {}
-
         self.counts = counts
         self.exposure = exposure
         self.background = background
@@ -354,7 +354,6 @@ def models(self, models):
                         use_cache=USE_NPRED_CACHE,
                     )
                     self._evaluators[model.name] = evaluator
-
         self._models = models
 
     @property
@@ -1114,12 +1113,11 @@ def plot_residuals(
         return ax_spatial, ax_spectral
 
     def stat_sum(self):
-        """Total statistic function value given the current model parameters and priors."""
+        """Total statistic function value given the current model parameters and priors set on the models."""
+        counts, npred = self.counts.data.astype(float), self.npred().data
         prior_stat_sum = 0.0
         if self.models is not None:
-            prior_stat_sum = self.models.parameters.prior_stat_sum()
-
-        counts, npred = self.counts.data.astype(float), self.npred().data
+            prior_stat_sum = prior_fit_statistic(self.models.priors)
 
         if self.mask is not None:
             return (

gammapy/datasets/tests/test_datasets.py

@@ -29,6 +29,9 @@ def test_datasets_types(datasets):
 
 
 def test_datasets_likelihood(datasets):
+
+    # rename the model of the second dataset for unique model names
+    datasets[1].models[0].name = datasets[1].name
     likelihood = datasets.stat_sum()
     assert_allclose(likelihood, 14472200.0002)
 

gammapy/datasets/tests/test_map.py

@@ -720,8 +720,12 @@ def test_prior_stat_sum(sky_model, geom, geom_etrue):
     datasets.models = models
     dataset.counts = dataset.npred()
 
-    uniformprior = UniformPrior(min=-np.inf, max=0, weight=1)
-    datasets.models.parameters["amplitude"].prior = uniformprior
+    UniformPrior(
+        min=-np.inf,
+        max=0,
+        weight=1,
+        modelparameters=datasets.models.parameters["amplitude"],
+    )
     assert_allclose(datasets.stat_sum(), 12825.9370, rtol=1e-3)
 
     datasets.models.parameters["amplitude"].value = -1e-12

gammapy/modeling/models/core.py

@@ -37,7 +37,6 @@
     """get a model class from a dict"""
     from . import (
         MODEL_REGISTRY,
-        PRIOR_REGISTRY,
         SPATIAL_MODEL_REGISTRY,
         SPECTRAL_MODEL_REGISTRY,
         TEMPORAL_MODEL_REGISTRY,
@@ -51,8 +50,6 @@
         cls = SPECTRAL_MODEL_REGISTRY.get_cls(data["spectral"]["type"])
     elif "temporal" in data:
         cls = TEMPORAL_MODEL_REGISTRY.get_cls(data["temporal"]["type"])
-    elif "prior" in data:
-        cls = PRIOR_REGISTRY.get_cls(data["prior"]["type"])
     return cls
 
 
@@ -403,6 +400,19 @@
         """Parameters (`~gammapy.modeling.Parameters`)"""
         return Parameters.from_stack([_.parameters for _ in self._models])
 
+    @property
+    def priors(self):
+        """Priors (list).
+
+        Duplicate prior objects have been removed.
+        """
+        priors = {}
+
+        for parameter in self.parameters:
+            if parameter.prior is not None:
+                priors[parameter.prior] = parameter.prior
+        return list(priors)
+
     @property
     def parameters_unique_names(self):
         """List of unique parameter names as model_name.par_type.par_name"""
@@ -1124,10 +1134,6 @@
 
         self._models.insert(idx, model)
 
-    def set_prior(self, parameters, priors):
-        for parameter, prior in zip(parameters, priors):
-            parameter.prior = prior
-
 
 class restore_models_status:
     def __init__(self, models, restore_values=True):

gammapy/modeling/models/prior.py

@@ -2,7 +2,7 @@
 import logging
 import numpy as np
 import astropy.units as u
-from gammapy.modeling import PriorParameter, PriorParameters
+from gammapy.modeling import Parameter, Parameters, PriorParameter, PriorParameters
 from .core import ModelBase
 
 log = logging.getLogger(__name__)
@@ -32,7 +32,15 @@
 class Prior(ModelBase):
     _unit = ""
 
-    def __init__(self, **kwargs):
+    def __init__(self, modelparameters, **kwargs):
+
+        if isinstance(modelparameters, Parameter):
+            self._modelparameters = Parameters([modelparameters])
+        elif isinstance(modelparameters, Parameters):
+            self._modelparameters = modelparameters
+        else:
+            raise ValueError(f"Invalid model type {modelparameters}")
+
         # Copy default parameters from the class to the instance
         default_parameters = self.default_parameters.copy()
 
@@ -52,6 +60,13 @@
         else:
             self._weight = 1
 
+        for par in self._modelparameters:
+            par.prior = self
+
+    @property
+    def modelparameters(self):
+        return self._modelparameters
+
     @property
     def parameters(self):
         """PriorParameters (`~gammapy.modeling.PriorParameters`)"""
@@ -73,11 +88,11 @@
     def weight(self, value):
         self._weight = value
 
-    def __call__(self, value):
+    def __call__(self):
         """Call evaluate method"""
-        # assuming the same unit as the PriorParamter here
+        # assuming the same unit as the PriorParameter here
         kwargs = {par.name: par.value for par in self.parameters}
-        return self.weight * self.evaluate(value.value, **kwargs)
+        return self.weight * self.evaluate(self._modelparameters.value, **kwargs)
 
     def to_dict(self, full_output=False):
         """Create dict for YAML serialisation"""
@@ -99,30 +114,33 @@
                     ):
                         del par[item]
 
-        data = {"type": tag, "parameters": params, "weight": self.weight}
+        data = {
+            "type": tag,
+            "parameters": params,
+            "weight": self.weight,
+            "modelparameters": self._modelparameters,
+        }
 
-        if self.type is None:
-            return data
-        else:
-            return {self.type: data}
+        return data
 
     @classmethod
     def from_dict(cls, data):
+        from . import PRIOR_REGISTRY
+
+        prior_cls = PRIOR_REGISTRY.get_cls(data["type"])
         kwargs = {}
 
-        key0 = next(iter(data))
-        if key0 in ["prior"]:
-            data = data[key0]
-        if data["type"] not in cls.tag:
+        if data["type"] not in prior_cls.tag:
             raise ValueError(
                 f"Invalid model type {data['type']} for class {cls.__name__}"
             )
-
         priorparameters = _build_priorparameters_from_dict(
-            data["parameters"], cls.default_parameters
+            data["parameters"], prior_cls.default_parameters
         )
         kwargs["weight"] = data["weight"]
-        return cls.from_parameters(priorparameters, **kwargs)
+        kwargs["modelparameters"] = data["modelparameters"]
+
+        return prior_cls.from_parameters(priorparameters, **kwargs)
 
 
 class GaussianPrior(Prior):
@@ -144,8 +162,7 @@
     mu = PriorParameter(name="mu", value=0)
     sigma = PriorParameter(name="sigma", value=1)
 
-    @staticmethod
-    def evaluate(value, mu, sigma):
+    def evaluate(self, value, mu, sigma):
         return ((value - mu) / sigma) ** 2
 
 
@@ -172,8 +189,7 @@
     min = PriorParameter(name="min", value=-np.inf, unit="")
     max = PriorParameter(name="max", value=np.inf, unit="")
 
-    @staticmethod
-    def evaluate(value, min, max):
+    def evaluate(self, value, min, max):
         if min < value < max:
             return 1.0
         else:

gammapy/modeling/models/tests/test_prior.py

@@ -1,23 +1,26 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import pytest
 from numpy.testing import assert_allclose
-import astropy.units as u
-from gammapy.modeling.models import PRIOR_REGISTRY, GaussianPrior, Model, UniformPrior
+from gammapy.modeling import Parameter
+from gammapy.modeling.models import PRIOR_REGISTRY, GaussianPrior, Prior, UniformPrior
 from gammapy.utils.testing import assert_quantity_allclose
 
+TEST_PARAMETER = [dict(testparameter=Parameter(name="test", value=1))]
+
+
 TEST_PRIORS = [
     dict(
         name="gaussian",
-        model=GaussianPrior(mu=4.0, sigma=1.0),
-        val_at_0=16.0,
-        val_at_1=9.0,
-        val_with_weight_2=32.0,
+        model=GaussianPrior(mu=2.3, sigma=1.0, modelparameters=Parameter("index", 2.3)),
+        val_at_default=0.0,
+        val_at_0=5.29,
+        val_with_weight_2=10.58,
     ),
     dict(
         name="uni",
-        model=UniformPrior(min=0.0),
+        model=UniformPrior(min=0.0, modelparameters=Parameter("amplitude", 1e-12)),
+        val_at_default=1.0,
         val_at_0=0.0,
-        val_at_1=1.0,
         val_with_weight_2=0.0,
     ),
 ]
@@ -29,14 +32,15 @@ def test_priors(prior):
     for p in model.parameters:
         assert p.type == "prior"
 
-    value_0 = model(0.0 * u.Unit(""))
-    value_1 = model(1.0 * u.Unit(""))
+    value_default = model()
+    model.modelparameters[0].value = 0.0
+    value_0 = model()
+    assert_allclose(value_default, prior["val_at_default"], rtol=1e-7)
     assert_allclose(value_0, prior["val_at_0"], rtol=1e-7)
-    assert_allclose(value_1, prior["val_at_1"], rtol=1e-7)
 
     model.weight = 2.0
-    value_0_weight = model(0.0 * u.Unit(""))
-    assert_allclose(value_0_weight, prior["val_with_weight_2"], rtol=1e-7)
+    value_2_weight = model()
+    assert_allclose(value_2_weight, prior["val_with_weight_2"], rtol=1e-7)
 
 
 def test_to_from_dict():
@@ -45,9 +49,9 @@ def test_to_from_dict():
     model.weight = 2.0
     model_dict = model.to_dict()
     # Here we reverse the order of parameters list to ensure assignment is correct
-    model_dict["prior"]["parameters"].reverse()
+    model_dict["parameters"].reverse()
 
-    model_class = PRIOR_REGISTRY.get_cls(model_dict["prior"]["type"])
+    model_class = PRIOR_REGISTRY.get_cls(model_dict["type"])
     new_model = model_class.from_dict(model_dict)
 
     assert isinstance(new_model, UniformPrior)
@@ -57,7 +61,7 @@ def test_to_from_dict():
     assert_quantity_allclose(actual, desired)
     assert_allclose(model.weight, new_model.weight, rtol=1e-7)
 
-    new_model = Model.from_dict(model_dict)
+    new_model = Prior.from_dict(model_dict)
 
     assert isinstance(new_model, UniformPrior)
 

gammapy/modeling/parameter.py

@@ -148,7 +148,7 @@
         self.scan_n_sigma = scan_n_sigma
         self.interp = interp
         self.scale_method = scale_method
-        self.prior = prior
+        self._prior = prior
 
     def __get__(self, instance, owner):
         if instance is None:
@@ -179,20 +179,23 @@
         if value is not None:
             from .models import Prior
 
+            # creating Prior out of dict
             if isinstance(value, dict):
-                from .models import Model
 
-                self._prior = Model.from_dict(value)
-            elif isinstance(value, Prior):
-                self._prior = value
+                value = Prior.from_dict(value)
+                if self in value.modelparameters:
+                    self._prior = value
+
+            # testing of self in prior.modelparameter
+            if isinstance(value, Prior):
+                if self in value.modelparameters:
+                    self._prior = value
+                else:
+                    raise TypeError(
+                        f"Invalid modelparameter: {value.modelparameters!r}"
+                    )
             else:
                 raise TypeError(f"Invalid type: {value!r}")
-        else:
-            self._prior = value
-
-    def prior_stat_sum(self):
-        if self.prior is not None:
-            return self.prior(self)
 
     @property
     def is_norm(self):
@@ -464,7 +467,7 @@
         if self._link_label_io is not None:
             output["link"] = self._link_label_io
         if self.prior is not None:
-            output["prior"] = self.prior.to_dict()["prior"]
+            output["prior"] = self.prior.to_dict()
         return output
 
     def autoscale(self):