Repository: gammapy/gammapy
Pull request #5639 (2024-12-10 17:02:04)

Add ParameterSensitivityEstimator and support sensitivity option on FluxpointsEstimator
- Add ParameterSensitivityEstimator a helper class used by the ParameterEstimator to compute sensitivity
 (and expose both classes)
- Add sensitivity option on FluxPointsEstimator

Requires #5638 to use the ts defined in Cowan et al.
Requires #5637 for the new implementation of ts_asimov that should fix an issue with ray.



Patch:
diff --git a/gammapy/datasets/actors.py b/gammapy/datasets/actors.py
index 42fe25089..31b6e36f4 100644
--- a/gammapy/datasets/actors.py
+++ b/gammapy/datasets/actors.py
@@ -84,6 +84,13 @@ class DatasetsActor(Datasets):
         results = self._ray_get([d._update_stat_sum_remote() for d in self._datasets])
         return np.sum(results)
 
+    def _to_asimov_datasets(self):
+        """Create Asimov datasets from the current models."""
+        asimov_datasets = Datasets([d._to_asimov_dataset() for d in self])
+        for d in asimov_datasets:
+            d._models = self.models
+        return asimov_datasets
+
 
 class RayFrontendMixin(object):
     """Ray mixin for a local class that interact with a remote instance."""
diff --git a/gammapy/datasets/map.py b/gammapy/datasets/map.py
index e692adbe2..49211fbd3 100644
--- a/gammapy/datasets/map.py
+++ b/gammapy/datasets/map.py
@@ -1493,14 +1493,8 @@ class MapDataset(Dataset):
         return cash_sum + prior_stat_sum
 
     def _to_asimov_dataset(self):
-        """Create Asimov dataset from the current models.
+        """Create Asimov dataset from the current models."""
 
-        Parameters
-        ----------
-        name : str, optional
-            Name of the new dataset. Default is None.
-
-        """
         npred = self.npred()
         data = np.nan_to_num(npred.data, copy=True, nan=0.0, posinf=0.0, neginf=0.0)
         npred.data = data.astype("float")
@@ -2748,14 +2742,7 @@ class MapDatasetOnOff(MapDataset):
         )
 
     def _to_asimov_dataset(self):
-        """Create Asimov dataset from the current models.
-
-        Parameters
-        ----------
-        name : str, optional
-            Name of the new dataset. Default is None.
-
-        """
+        """Create Asimov dataset from the current models."""
         npred = self.npred()
         data = np.nan_to_num(npred.data, copy=True, nan=0.0, posinf=0.0, neginf=0.0)
         npred.data = data.astype("float")
diff --git a/gammapy/estimators/__init__.py b/gammapy/estimators/__init__.py
index 3a38c2ff2..8ee826bee 100644
--- a/gammapy/estimators/__init__.py
+++ b/gammapy/estimators/__init__.py
@@ -1,11 +1,12 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Estimators."""
+
 from gammapy.utils.registry import Registry
 from .core import Estimator
 from .energydependentmorphology import EnergyDependentMorphologyEstimator
 from .map import ASmoothMapEstimator, ExcessMapEstimator, FluxMaps, TSMapEstimator
 from .metadata import FluxMetaData
-from .parameter import ParameterEstimator
+from .parameter import ParameterEstimator, ParameterSensitivityEstimator
 from .points import (
     FluxPoints,
     FluxPointsEstimator,
@@ -29,6 +30,7 @@ __all__ = [
     "LightCurveEstimator",
     "ParameterEstimator",
     "SensitivityEstimator",
+    "ParameterSensitivityEstimator",
     "TSMapEstimator",
     "EnergyDependentMorphologyEstimator",
     "FluxMetaData",
@@ -43,6 +45,7 @@ ESTIMATOR_REGISTRY = Registry(
         ASmoothMapEstimator,
         LightCurveEstimator,
         SensitivityEstimator,
+        ParameterSensitivityEstimator,
         FluxProfileEstimator,
         ParameterEstimator,
     ]
diff --git a/gammapy/estimators/flux.py b/gammapy/estimators/flux.py
index a056b4065..619233bad 100644
--- a/gammapy/estimators/flux.py
+++ b/gammapy/estimators/flux.py
@@ -31,6 +31,8 @@ class FluxEstimator(ParameterEstimator):
         Sigma to use for asymmetric error computation.
     n_sigma_ul : int
         Sigma to use for upper limit computation.
+    n_sigma_sensitivity : int
+        Sigma to use for sensitivity computation.
     selection_optional : list of str, optional
         Which additional quantities to estimate. Available options are:
 
@@ -66,6 +68,7 @@ class FluxEstimator(ParameterEstimator):
         source=0,
         n_sigma=1,
         n_sigma_ul=2,
+        n_sigma_sensitivity=5,
         selection_optional=None,
         fit=None,
         reoptimize=False,
@@ -79,6 +82,7 @@ class FluxEstimator(ParameterEstimator):
             null_value=0,
             n_sigma=n_sigma,
             n_sigma_ul=n_sigma_ul,
+            n_sigma_sensitivity=n_sigma_sensitivity,
             selection_optional=selection_optional,
             fit=fit,
             reoptimize=reoptimize,
diff --git a/gammapy/estimators/map/core.py b/gammapy/estimators/map/core.py
index 29ae80f48..18f15431a 100644
--- a/gammapy/estimators/map/core.py
+++ b/gammapy/estimators/map/core.py
@@ -61,7 +61,7 @@ OPTIONAL_QUANTITIES = {
     "e2dnde": ["e2dnde_err", "e2dnde_errp", "e2dnde_errn", "e2dnde_ul"],
     "flux": ["flux_err", "flux_errp", "flux_errn", "flux_ul", "flux_sensitivity"],
     "eflux": ["eflux_err", "eflux_errp", "eflux_errn", "eflux_ul"],
-    "likelihood": ["norm_err", "norm_errn", "norm_errp", "norm_ul"],
+    "likelihood": ["norm_err", "norm_errn", "norm_errp", "norm_ul", "norm_sensitivity"],
 }
 
 VALID_QUANTITIES = [
diff --git a/gammapy/estimators/parameter.py b/gammapy/estimators/parameter.py
index fc8219cc1..5aec56735 100644
--- a/gammapy/estimators/parameter.py
+++ b/gammapy/estimators/parameter.py
@@ -1,13 +1,21 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
+import warnings
+
 import logging
 import numpy as np
 from gammapy.datasets import Datasets
 from gammapy.datasets.actors import DatasetsActor
 from gammapy.modeling import Fit
+from gammapy.modeling.selection import TestStatisticNested
+from gammapy.modeling.parameter import restore_parameters_status
+from gammapy.stats.utils import ts_to_sigma
+from gammapy.utils.roots import find_roots
 from .core import Estimator
 
 log = logging.getLogger(__name__)
 
+__all__ = ["ParameterEstimator", "ParameterSensitivityEstimator"]
+
 
 class ParameterEstimator(Estimator):
     """Model parameter estimator.
@@ -22,6 +30,8 @@ class ParameterEstimator(Estimator):
         Sigma to use for asymmetric error computation. Default is 1.
     n_sigma_ul : int
         Sigma to use for upper limit computation. Default is 2.
+    n_sigma_sensitivity : int
+        Sigma to use for sensitivity computation. Default is 5.
     null_value : float
         Which null value to use for the parameter.
     selection_optional : list of str, optional
@@ -60,12 +70,13 @@ class ParameterEstimator(Estimator):
     """
 
     tag = "ParameterEstimator"
-    _available_selection_optional = ["errn-errp", "ul", "scan"]
+    _available_selection_optional = ["errn-errp", "ul", "scan", "sensitivity"]
 
     def __init__(
         self,
         n_sigma=1,
         n_sigma_ul=2,
+        n_sigma_sensitivity=5,
         null_value=1e-150,
         selection_optional=None,
         fit=None,
@@ -73,6 +84,7 @@ class ParameterEstimator(Estimator):
     ):
         self.n_sigma = n_sigma
         self.n_sigma_ul = n_sigma_ul
+        self.n_sigma_sensitivity = n_sigma_sensitivity
         self.null_value = null_value
         self.selection_optional = selection_optional
 
@@ -262,6 +274,27 @@ class ParameterEstimator(Estimator):
         )
         return {f"{parameter.name}_ul": res["errp"] + parameter.value}
 
+    def estimate_sensitivity(self, datasets, parameter):
+        """Estimate norm sensitivity for the flux point.
+
+        Parameters
+        ----------
+        datasets : `~gammapy.datasets.Datasets`
+            Datasets.
+
+        Returns
+        -------
+        result : dict
+            Dictionary with an array with one entry per dataset with the sum of the
+            masked npred.
+        """
+
+        estimator = ParameterSensitivityEstimator(
+            parameter, self.null_value, n_sigma=self.n_sigma_sensitivity
+        )
+        value = estimator.run(datasets)
+        return {f"{parameter.name}_sensitivity": value}
+
     @staticmethod
     def estimate_counts(datasets):
         """Estimate counts for the flux point.
@@ -328,7 +361,6 @@ class ParameterEstimator(Estimator):
         parameter = datasets.parameters[parameter]
 
         with datasets.parameters.restore_status():
-
             if not self.reoptimize:
                 datasets.parameters.freeze_all()
                 parameter.frozen = False
@@ -345,5 +377,107 @@ class ParameterEstimator(Estimator):
             if "scan" in self.selection_optional:
                 result.update(self.estimate_scan(datasets, parameter))
 
+            if "sensitivity" in self.selection_optional:
+                result.update(self.estimate_sensitivity(datasets, parameter))
+
         result.update(self.estimate_counts(datasets))
         return result
+
+
+class ParameterSensitivityEstimator:
+    """Estimate the sensitivity to a given parameter
+
+    Computes the TS distribution in the non-null hypothesis using the
+    log likelihood of the Asimov dataset (i.e. a dataset with counts = npred)
+    and the non central chi2 distribution.
+    Once the TS distribution under the testing hypothesis is known,
+    one can compute the required parameter value
+    to have 50% of measurements above a given significance threshold.
+
+
+    Parameters
+    ----------
+    parameter : `~gammapy.modeling.Parameter`
+       Parameter to test
+    null_value : float or `~gammapy.modeling.Parameter`
+        Value of the parameter for the null hypothesis.
+    n_sigma : int, default=5
+        Number of required significance level.
+    rtol : float
+        Relative precision of the estimate. Used as a stopping criterion.
+        Default is 0.01.
+    max_niter : int
+        Maximal number of iterations used by the root finding algorithm.
+        Default is 100.
+
+    References
+    ----------
+        * `Cowan et al. (2011), "Asymptotic formulae for likelihood-based tests of new physics"
+        <https://arxiv.org/abs/1007.1727>`_
+
+    """
+
+    tag = "ParameterSensitivityEstimator"
+
+    def __init__(
+        self,
+        parameter,
+        null_value,
+        n_sigma=5,
+        n_free_parameters=None,
+        rtol=0.01,
+        max_niter=100,
+    ):
+        self.test = TestStatisticNested(
+            [parameter], [null_value], n_free_parameters=n_free_parameters
+        )
+        self.parameter = parameter
+        self.n_sigma = n_sigma
+        self.rtol = rtol
+        self.max_niter = max_niter
+
+    def _fcn(self, value, datasets):
+        """Call the Test Statistics function."""
+        self.parameter.value = value
+        ts_asimov = self.test.ts_asimov(datasets)
+        return ts_to_sigma(ts_asimov, ts_asimov=ts_asimov) - self.n_sigma
+
+    def parameter_matching_significance(self, datasets):
+        """Parameter value  matching the target significance"""
+
+        if ~np.isfinite(self.parameter.min):
+            vmin = self.parameter.value / 1e3
+        else:
+            vmin = self.parameter.min
+        if ~np.isfinite(self.parameter.max):
+            vmax = self.parameter.value * 1e3
+        else:
+            vmax = self.parameter.max
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            roots, res = find_roots(
+                self._fcn,
+                vmin,
+                vmax,
+                args=(datasets,),
+                nbin=100,
+                maxiter=self.max_niter,
+                rtol=self.rtol,
+                points_scale=self.parameter.interp,
+            )
+        # Where the root finding fails NaN is set as norm
+        roots = roots[roots > 0]
+        if roots.size > 0:
+            return roots[0]
+        else:
+            return np.nan
+
+    def run(self, datasets):
+        """Parameter sensitivity
+        given as the difference between value matching the target significance and the null value.
+        """
+        with restore_parameters_status(self.test.parameters):
+            value = self.parameter_matching_significance(datasets)
+
+        return value - self.test.null_values[0]
diff --git a/gammapy/estimators/points/sed.py b/gammapy/estimators/points/sed.py
index fc56b03ad..3180520f4 100644
--- a/gammapy/estimators/points/sed.py
+++ b/gammapy/estimators/points/sed.py
@@ -41,6 +41,8 @@ class FluxPointsEstimator(FluxEstimator, parallel.ParallelMixin):
         Number of sigma to use for asymmetric error computation. Default is 1.
     n_sigma_ul : int
         Number of sigma to use for upper limit computation. Default is 2.
+    n_sigma_sensitivity : int
+        Sigma to use for sensitivity computation. Default is 5.
     selection_optional : list of str, optional
         Which additional quantities to estimate. Available options are:
 
@@ -48,6 +50,7 @@ class FluxPointsEstimator(FluxEstimator, parallel.ParallelMixin):
             * "errn-errp": estimate asymmetric errors on flux.
             * "ul": estimate upper limits.
             * "scan": estimate fit statistic profiles.
+            * "sensitivity": estimate sensitivity for a given significance.
 
         Default is None so the optional steps are not executed.
     energy_edges : list of `~astropy.units.Quantity`, optional
@@ -229,4 +232,7 @@ class FluxPointsEstimator(FluxEstimator, parallel.ParallelMixin):
             norm_scan = self.norm.copy().scan_values
             result.update({"norm_scan": norm_scan, "stat_scan": np.nan * norm_scan})
 
+        if "sensitivity" in self.selection_optional:
+            result.update({"norm_sensitivity": np.nan})
+
         return result


Test patch:
diff --git a/gammapy/estimators/points/tests/test_sed.py b/gammapy/estimators/points/tests/test_sed.py
index 77f96e030..29c8a4db9 100644
--- a/gammapy/estimators/points/tests/test_sed.py
+++ b/gammapy/estimators/points/tests/test_sed.py
@@ -269,6 +269,9 @@ def test_run_pwl(fpe_pwl, tmpdir):
     actual = table["norm_ul"].data
     assert_allclose(actual, [1.216227, 1.035472, 1.316878], rtol=1e-2)
 
+    actual = table["norm_sensitivity"].data
+    assert_allclose(actual, [0.15901235, 0.13117017, 0.55880332], rtol=1e-2)
+
     actual = table["sqrt_ts"].data
     assert_allclose(actual, [18.568429, 18.054651, 7.057121], rtol=1e-2)
 
diff --git a/gammapy/estimators/points/tests/test_sensitivity.py b/gammapy/estimators/points/tests/test_sensitivity.py
index bc3fecb99..03d4bb83e 100644
--- a/gammapy/estimators/points/tests/test_sensitivity.py
+++ b/gammapy/estimators/points/tests/test_sensitivity.py
@@ -1,11 +1,15 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import pytest
 from numpy.testing import assert_allclose
-from gammapy.datasets import SpectrumDataset, SpectrumDatasetOnOff
-from gammapy.estimators import FluxPoints, SensitivityEstimator
+from gammapy.datasets import SpectrumDataset, SpectrumDatasetOnOff, Datasets
+from gammapy.estimators import (
+    FluxPoints,
+    SensitivityEstimator,
+    ParameterSensitivityEstimator,
+)
 from gammapy.irf import EDispKernelMap
 from gammapy.maps import MapAxis, RegionNDMap
-from gammapy.modeling.models import PowerLawSpectralModel
+from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
 from gammapy.utils.deprecation import GammapyDeprecationWarning
 
 
@@ -26,9 +30,10 @@ def spectrum_dataset():
         region="icrs;circle(0, 0, 0.1)", axes=[e_true], unit="m2 h", data=1e6
     )
 
-    return SpectrumDataset(
+    spectrum_dataset = SpectrumDataset(
         name="test", exposure=exposure, edisp=edisp, background=background
     )
+    return spectrum_dataset
 
 
 def test_cta_sensitivity_estimator(spectrum_dataset, caplog):
@@ -100,3 +105,18 @@ def test_integral_estimation(spectrum_dataset, caplog):
 
     assert_allclose(table["excess"].data.squeeze(), 270540, rtol=1e-3)
     assert_allclose(flux_points.flux.data.squeeze(), 7.52e-9, rtol=1e-3)
+
+
+def test_parameter_sensitivity_estimator(spectrum_dataset):
+    spectral_model = PowerLawSpectralModel()
+    default_value = spectral_model.amplitude.value
+
+    spectrum_dataset.models = SkyModel(spectral_model=spectral_model)
+    datasets = Datasets(spectrum_dataset)
+
+    estimator = ParameterSensitivityEstimator(spectral_model.amplitude, 0, rtol=1e-2)
+
+    value = estimator.run(datasets)
+    assert_allclose(value, 4.570882e-12, rtol=1e-2)
+
+    assert_allclose(spectral_model.amplitude.value, default_value, rtol=1e-2)