Merge pull request #89 from cta-observatory/fix_lima

Fix brentq returning the lower bound of 0 for flat li ma function

pyirf/sensitivity.py

@@ -25,7 +25,6 @@ def relative_sensitivity(
     alpha,
     target_significance=5,
     significance_function=li_ma_significance,
-    initial_guess=0.01,
 ):
     """
     Calculate the relative sensitivity defined as the flux
@@ -62,20 +61,18 @@ def relative_sensitivity(
         "Analysis methods for results in gamma-ray astronomy."
         The Astrophysical Journal 272 (1983): 317-324.
         Formula (17)
-    initial_guess: float
-        Initial guess for the root finder
     """
-    n_background = n_off * alpha
-    n_signal = n_on - n_background
-
     if np.isnan(n_on) or np.isnan(n_off):
         return np.nan
 
-    if n_on < 1 or n_off < 1:
-        return np.nan
+    if n_on < 0 or n_off < 0:
+        raise ValueError(f'n_on and n_off must be positive, got {n_on}, {n_off}')
+
+    n_background = n_off * alpha
+    n_signal = n_on - n_background
 
     if n_signal <= 0:
-        return np.nan
+        return np.inf
 
     def equation(relative_flux):
         n_on = n_signal * relative_flux + n_background
@@ -84,20 +81,27 @@ def equation(relative_flux):
 
     try:
         # brentq needs a lower and an upper bound
-        # lower can be trivially set to zero, but the upper bound is more tricky
         # we will use the simple, analytically  solvable significance formula and scale it
         # with 10 to be sure it's above the Li and Ma solution
         # so rel * n_signal / sqrt(n_background) = target_significance
-        upper_bound = 10 * target_significance * np.sqrt(n_background) / n_signal
-        result = brentq(equation, 0, upper_bound,)
-    except (RuntimeError, ValueError):
+        if n_off > 1:
+            relative_flux_naive = target_significance * np.sqrt(n_background) / n_signal
+            upper_bound = 10 * relative_flux_naive
+            lower_bound = 0.01 * relative_flux_naive
+        else:
+            upper_bound = 100
+            lower_bound = 1e-4
+
+        relative_flux = brentq(equation, lower_bound, upper_bound)
+
+    except (RuntimeError, ValueError) as e:
         log.warn(
             "Could not calculate relative significance for"
-            f" n_signal={n_signal:.1f}, n_off={n_off:.1f}, returning nan"
+            f" n_signal={n_signal:.1f}, n_off={n_off:.1f}, returning nan {e}"
         )
         return np.nan
 
-    return result
+    return relative_flux
 
 
 def calculate_sensitivity(

pyirf/statistics.py

@@ -41,12 +41,14 @@ def li_ma_significance(n_on, n_off, alpha=0.2):
         t2 = n_off * np.log((1 + alpha) * p_off)
 
         # lim x+->0  (x log(x)) = 0
+        t1[n_on == 0] = 0
         t2[n_off == 0] = 0
+
         ts = t1 + t2
 
         significance = np.sqrt(ts * 2)
 
-    significance[n_on < alpha * n_off] = 0
+    significance[n_on < (alpha * n_off)] = 0
 
     if scalar:
         return significance[0]

pyirf/tests/test_sensitivity.py

@@ -3,6 +3,30 @@
 import astropy.units as u
 
 
+def test_relative_sensitivity():
+    from pyirf.sensitivity import relative_sensitivity
+
+    # some general case
+    n_on = 100
+    n_off = 200
+    alpha = 0.2
+    assert 0.1 < relative_sensitivity(n_on, n_off, alpha) < 1.0
+
+    # numbers yield lima = 5 relatively precisely, so sensitivity should be 1
+    assert np.isclose(relative_sensitivity(81, 202, 0.2), 1, rtol=0.01)
+
+    # test different target significance
+    # numbers yield lima = 8 relatively precisely, so sensitivity should be 1
+    result = relative_sensitivity(93, 151, 0.2, target_significance=8)
+    assert np.isclose(result, 1, rtol=0.01)
+
+    # no signal => inf
+    assert np.isinf(relative_sensitivity(10, 100, 0.2))
+
+    # no background, should work
+    assert relative_sensitivity(10, 0, 0.2) > 0
+
+
 def test_estimate_background():
     from pyirf.sensitivity import estimate_background
     N = 1000