Repository: gammapy/gammapy
Pull request #5669 (2025-01-15 13:25:29)

Add Log uniform prior
This PR is a brick of the larger Nested Sampling discussed in the PIG https://github.com/gammapy/gammapy/pull/5499

We introduce here the `LogUniformPrior` that is critically needed for sampling for example normalizations that span mulitple order of magnitude values.

**Implementation in this PR**

We add the `LogUniformPrior` in a similar way to the `UniformPrior` or `GaussianPrior` with an `evaluate` and `inverse_cdf` methods.
For now for the Nested sampling, only the inverse cdf is needed.   The `evaluate` method has been implemented but raised several questions (see next section)

I needed to change the dictionary used for the tests to accomodate the 3 priors in a similar way. For example the LogUniformPrior can not be evaluated at 0. 
I also removed the tests on the weight for `GaussianPrior` and `LogUniformPrior` as the weights should not affect the Gaussian (see next section).

**General issues raised while preparing this PR**
While preparing this PR with @registerrier we stumbled on two concepts developed in the Priors:

- weights : While it makes sense to have a weight attribute to the UniformPrior (a likelihood weight/penalty outside the min,max range), it is not a good idea for the Gaussian Prior. Multiplying the Gaussian by a weight value would defacto modify the sigma of this Gaussian... So weight should be used only for the UniformPrior. right now weight is applied at the class level upon call so this should be modified to be multiplied only locally in the prior I think.
- evaluate : the way the current priors are implemented, you get a zero likelihood penalty at the "good" values (in the min,max range for the Uniform or at mu for Gaussian) and then a penalty outside this range. However it is not clear what this value should be in the case of the LogUniformPrior. It cannot be a constant value between min, max. Using scipy.stat.loguniform, one can get the probability density function (pdf) which will but more penalty for higher values than low values. But this implies that there will always be a non-zero penalty and hence the Likelihood with/without prior will be different.

This discussion is probably outside the scope of this PR and could be moved to an issue.



Patch:
diff --git a/gammapy/modeling/models/__init__.py b/gammapy/modeling/models/__init__.py
index 560bfcc15..ba1c305e2 100644
--- a/gammapy/modeling/models/__init__.py
+++ b/gammapy/modeling/models/__init__.py
@@ -9,7 +9,7 @@ from .cube import (
     TemplateNPredModel,
     create_fermi_isotropic_diffuse_model,
 )
-from .prior import GaussianPrior, Prior, UniformPrior
+from .prior import GaussianPrior, Prior, UniformPrior, LogUniformPrior
 from .spatial import (
     ConstantFluxSpatialModel,
     ConstantSpatialModel,
@@ -112,6 +112,7 @@ __all__ = [
     "Prior",
     "GaussianPrior",
     "UniformPrior",
+    "LogUniformPrior",
     "scale_plot_flux",
     "ScaleSpectralModel",
     "Shell2SpatialModel",
diff --git a/gammapy/modeling/models/prior.py b/gammapy/modeling/models/prior.py
index c09d9158e..c6bd7a750 100644
--- a/gammapy/modeling/models/prior.py
+++ b/gammapy/modeling/models/prior.py
@@ -4,11 +4,11 @@
 import logging
 import numpy as np
 import astropy.units as u
-from scipy.stats import norm, uniform
+from scipy.stats import norm, uniform, loguniform
 from gammapy.modeling import PriorParameter, PriorParameters
 from .core import ModelBase
 
-__all__ = ["GaussianPrior", "UniformPrior", "Prior"]
+__all__ = ["GaussianPrior", "UniformPrior", "LogUniformPrior", "Prior"]
 
 log = logging.getLogger(__name__)
 
@@ -196,3 +196,37 @@ class UniformPrior(Prior):
         """Return inverse CDF for prior."""
         rv = uniform(self.min.value, self.max.value - self.min.value)
         return rv.ppf(value)
+
+
+class LogUniformPrior(Prior):
+    """LogUniform Prior.
+
+     Equivalent to a uniform prior on the log of the parameter
+
+     Parameters
+     ----------
+     min : float
+         minimum value.
+    max : float
+         maximum value.
+    """
+
+    tag = ["LogUniformPrior"]
+    _type = "prior"
+    min = PriorParameter(name="min", value=1e-14, unit="")
+    max = PriorParameter(name="max", value=1e-10, unit="")
+
+    @staticmethod
+    def evaluate(value, min, max):
+        """
+        Evaluate the likelihood penalization term (hence -2*).
+        Note that this is currently a different scaling that the Uniform or Gaussian priors.
+        With current implementation the TS of a source with/without LogUniform prior would be different... TBD
+        """
+        rv = loguniform(min, max)
+        return -2 * rv.logpdf(value)
+
+    def _inverse_cdf(self, value):
+        """Return inverse CDF for prior."""
+        rv = loguniform(self.min.value, self.max.value)
+        return rv.ppf(value)


Test patch:
diff --git a/gammapy/modeling/models/tests/test_prior.py b/gammapy/modeling/models/tests/test_prior.py
index 25e1d6a15..110a22a31 100644
--- a/gammapy/modeling/models/tests/test_prior.py
+++ b/gammapy/modeling/models/tests/test_prior.py
@@ -3,51 +3,81 @@ import pytest
 import numpy as np
 from numpy.testing import assert_allclose
 import astropy.units as u
-from gammapy.modeling.models import PRIOR_REGISTRY, GaussianPrior, Model, UniformPrior
+from gammapy.modeling.models import (
+    PRIOR_REGISTRY,
+    GaussianPrior,
+    Model,
+    UniformPrior,
+    LogUniformPrior,
+)
 from gammapy.utils.testing import assert_quantity_allclose
 
 TEST_PRIORS = [
     dict(
         name="gaussian",
         model=GaussianPrior(mu=4.0, sigma=1.0),
+        prior_0=0.0 * u.Unit(""),
+        prior_1=1.0 * u.Unit(""),
         val_at_0=16.0,
         val_at_1=9.0,
-        val_with_weight_2=32.0,
         inverse_cdf_at_0=-np.inf,
         inverse_cdf_at_1=np.inf,
     ),
     dict(
-        name="uni",
+        name="uniform",
         model=UniformPrior(min=0.0, max=10),
+        prior_0=0.0 * u.Unit(""),
+        prior_1=1.0 * u.Unit(""),
         val_at_0=1.0,
         val_at_1=0.0,
         val_with_weight_2=2.0,
         inverse_cdf_at_0=0.0,
         inverse_cdf_at_1=10.0,
     ),
+    dict(
+        name="loguniform",
+        model=LogUniformPrior(min=1e-14, max=1e-10),
+        prior_0=1e-10 * u.Unit(""),
+        prior_1=1e-14 * u.Unit(""),
+        val_at_0=-41.61104824714522,
+        val_at_1=-60.03172899109759,
+        inverse_cdf_at_0=1e-14,
+        inverse_cdf_at_1=1e-10,
+    ),
 ]
 
 
 @pytest.mark.parametrize("prior", TEST_PRIORS)
-def test_priors(prior):
+def test_prior_evaluation(prior):
     model = prior["model"]
+    # Test the evaluation of the prior at specific points
+    assert_allclose(model(prior["prior_0"]), prior["val_at_0"], rtol=1e-7)
+    assert_allclose(model(prior["prior_1"]), prior["val_at_1"], rtol=1e-7)
+
+    # Test the inverse_cdf at specific points
+    value_0 = model._inverse_cdf(0)
+    value_1 = model._inverse_cdf(1)
+    assert_allclose(value_0, prior["inverse_cdf_at_0"], rtol=1e-7)
+    assert_allclose(value_1, prior["inverse_cdf_at_1"], rtol=1e-7)
+
+
+@pytest.mark.parametrize("prior", TEST_PRIORS)
+def test_prior_parameters(prior):
+    model = prior["model"]
+    # Check that all parameters of the model have type 'prior'
     for p in model.parameters:
         assert p.type == "prior"
 
-    value_0 = model(0.0 * u.Unit(""))
-    value_1 = model(1.0 * u.Unit(""))
-    assert_allclose(value_0, prior["val_at_0"], rtol=1e-7)
-    assert_allclose(value_1, prior["val_at_1"], rtol=1e-7)
 
+@pytest.mark.parametrize("prior", TEST_PRIORS)
+def test_uniform_prior_weight(prior):
+    prior = TEST_PRIORS[1]
+    model = prior["model"]
+    # Test the uniform prior with a specific weight
     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_0 = model._inverse_cdf(0.0)
-    value_1 = model._inverse_cdf(1.0)
-    assert_allclose(value_0, prior["inverse_cdf_at_0"], rtol=1e-7)
-    assert_allclose(value_1, prior["inverse_cdf_at_1"], rtol=1e-7)
-
 
 def test_to_from_dict():
     prior = TEST_PRIORS[1]