Merge pull request #1084 from adonath/flux_points_ul_estimation

Implement flux points upper limit estimation

gammapy/catalog/hess.py

@@ -345,7 +345,7 @@ def _info_flux_points(self):
         flux_points = self.flux_points.table.copy()
 
         energy_cols = ['e_ref', 'e_min', 'e_max']
-        flux_cols = ['dnde', 'dnde_errn', 'dnde_errp']
+        flux_cols = ['dnde', 'dnde_errn', 'dnde_errp', 'dnde_ul']
         flux_points = flux_points[energy_cols + flux_cols]
 
         for _ in energy_cols:

gammapy/spectrum/flux_point.py

@@ -557,7 +557,6 @@ def compute_flux_point(self, energy_group):
             global_model=self.model,
             energy_range=energy_group.energy_range,
         )
-
         energy_ref = self.compute_energy_ref(energy_group)
 
         return self.fit_point(
@@ -617,9 +616,110 @@ def compute_approx_model(global_model, energy_range):
                 par.frozen = True
         return approx_model
 
-    def fit_point(self, model, energy_group, energy_ref):
+    def compute_flux_point_ul(self, fit, best_fit, delta_ts=4, negative=False):
+        """
+        Compute upper limits for flux point values.
+
+
+        Parameters
+        ----------
+        fit : `SpectrumFit`
+            Instance of spectrum fit.
+        best_fit : `SpectrumFitResult`
+            Best fit result.
+        delta_ts : float (4)
+            Difference in log-likelihood for given confidence interval.
+            See Example below.
+        negative : bool
+            Compute limit in negative direction.
+
+
+        Examples
+        --------
+        To compute ~95% confidence upper limits (or 2 sigma) you can use:
+
+            from scipy.stats import chi2, norm
+
+            sigma = 2
+            cl = 1 - 2 * norm.sf(sigma) # using two sided p-value
+            delta_ts = chi2.isf(1 - cl, df=1)
+
+
+        Returns
+        -------
+        dnde_ul : `~astropy.units.Quantity`
+            Flux point upper limit.
+
+        """
+
+        from scipy.optimize import brentq
+
+        # this is a prototype for fast flux point upper limit
+        # calculation using brentq
+        stat_best_fit = best_fit.statval
+        amplitude = best_fit.model.parameters['amplitude'].value / 1E-12
+        amplitude_err = best_fit.model.parameters.error('amplitude') / 1E-12
+
+        if negative:
+            amplitude_max = amplitude
+            amplitude_min = amplitude_max - 1E3 * amplitude_err
+        else:
+            amplitude_max = amplitude + 1E3 * amplitude_err
+            amplitude_min = amplitude
+
+        def ts_diff(x):
+            fit.model.parameters['amplitude'].value = x * 1E-12
+            fit.predict_counts()
+            fit.calc_statval()
+            return (stat_best_fit + delta_ts) - fit.total_stat
+
+        try:
+            result = brentq(ts_diff, amplitude_min, amplitude_max,
+                            maxiter=100, rtol=1e-2)
+            return 1E-12 * result * fit.model.parameters['amplitude'].unit
+        except (RuntimeError, ValueError):
+            # Where the root finding fails NaN is set as amplitude
+            log.debug('Flux point upper limit computation failed.')
+            return np.nan * fit.model.parameters['amplitude'].unit
+
+    def compute_flux_point_sqrt_ts(self, fit, best_fit):
+        """
+        Compute sqrt(TS) for flux point.
+
+
+        Parameters
+        ----------
+        fit : `SpectrumFit`
+            Instance of spectrum fit.
+        best_fit : `SpectrumFitResult`
+            Best fit result.
+
+
+        Returns
+        -------
+        sqrt_ts : float
+            Sqrt(TS) for flux point.
+
+        """
+        amplitude = best_fit.model.parameters['amplitude'].value
+        stat_best_fit = best_fit.statval
+
+        fit.model.parameters['amplitude'].value = 0
+        fit.predict_counts()
+        fit.calc_statval()
+        stat_null = fit.total_stat
+
+        fit.model.parameters['amplitude'].value = amplitude
+        ts = np.abs(stat_null - stat_best_fit)
+        return np.sign(amplitude) * np.sqrt(ts)
+
+    def fit_point(self, model, energy_group, energy_ref, sqrt_ts_threshold=1):
         from .fit import SpectrumFit
 
+        # Set reference and remove min amplitude
+        model.parameters['reference'].value = energy_ref.to('TeV').value
+        model.parameters['amplitude'].parmin = None
+
         fit = SpectrumFit(self.obs, model)
         erange = energy_group.energy_range
 
@@ -640,13 +740,23 @@ def fit_point(self, model, energy_group, energy_ref):
         e_max = energy_group.energy_range.max
         e_min = energy_group.energy_range.min
         dnde, dnde_err = res.model.evaluate_error(energy_ref)
+        sqrt_ts = self.compute_flux_point_sqrt_ts(fit, best_fit=res)
+
+        dnde_ul = self.compute_flux_point_ul(fit, best_fit=res)
+        dnde_errp = self.compute_flux_point_ul(fit, best_fit=res, delta_ts=1.) - dnde
+        dnde_errn = dnde - self.compute_flux_point_ul(fit, best_fit=res, delta_ts=1., negative=True)
 
         return OrderedDict([
             ('e_ref', energy_ref),
             ('e_min', e_min),
             ('e_max', e_max),
             ('dnde', dnde.to(DEFAULT_UNIT['dnde'])),
             ('dnde_err', dnde_err.to(DEFAULT_UNIT['dnde'])),
+            ('dnde_ul', dnde_ul.to(DEFAULT_UNIT['dnde'])),
+            ('is_ul', sqrt_ts < sqrt_ts_threshold),
+            ('sqrt_ts', sqrt_ts),
+            ('dnde_errp', dnde_errp),
+            ('dnde_errn', dnde_errn)
         ])
 
 

gammapy/spectrum/tests/test_flux_point.py

@@ -164,10 +164,13 @@ def test_flux_points(config):
             ),
             obs=obs(),
             seg=seg(obs()),
-            dnde=2.7465439050126e-11 * u.Unit('cm-2 s-1 TeV-1'),
-            dnde_err=4.755502901867284e-12 * u.Unit('cm-2 s-1 TeV-1'),
-            res=-0.11262182922477647,
-            res_err=0.1536450758523701,
+            dnde=2.7465e-11 * u.Unit('cm-2 s-1 TeV-1'),
+            dnde_err=4.7555e-12 * u.Unit('cm-2 s-1 TeV-1'),
+            dnde_errn=4.5333e-12 * u.Unit('cm-2 s-1 TeV-1'),
+            dnde_errp=5.0050e-12 * u.Unit('cm-2 s-1 TeV-1'),
+            dnde_ul=3.7998e-11 * u.Unit('cm-2 s-1 TeV-1'),
+            res=-0.1126,
+            res_err=0.1536,
         )
 
     tester = FluxPointTester(config)
@@ -181,6 +184,7 @@ def test_flux_points(config):
 class FluxPointTester:
     def __init__(self, config):
         self.config = config
+        self.rtol = 0.5E-2  # accuracy of 0.5%
         self.setup()
 
     def setup(self):
@@ -217,11 +221,24 @@ def test_values(self):
 
         actual = flux_points.table['dnde'].quantity[0]
         desired = self.config['dnde']
-        assert_quantity_allclose(actual, desired)
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
 
         actual = flux_points.table['dnde_err'].quantity[0]
         desired = self.config['dnde_err']
-        assert_quantity_allclose(actual, desired)
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
+
+        actual = flux_points.table['dnde_ul'].quantity[0]
+        desired = self.config['dnde_ul']
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
+
+        actual = flux_points.table['dnde_errn'].quantity[0]
+        desired = self.config['dnde_errn']
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
+
+        actual = flux_points.table['dnde_errp'].quantity[0]
+        desired = self.config['dnde_errp']
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
+
 
     def test_spectrum_result(self):
         result = SpectrumResult(
@@ -231,11 +248,11 @@ def test_spectrum_result(self):
 
         actual = result.flux_point_residuals[0][0]
         desired = self.config['res']
-        assert_quantity_allclose(actual, desired)
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
 
         actual = result.flux_point_residuals[1][0]
         desired = self.config['res_err']
-        assert_quantity_allclose(actual, desired)
+        assert_quantity_allclose(actual, desired, rtol=self.rtol)
 
         result.plot(energy_range=[1, 10] * u.TeV)