Multiple quantiles for angular_resolution

docs/changes/290.feature.rst

@@ -0,0 +1,2 @@
+Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them.
+If only one quantile is passed, the output is the same as before.

pyirf/benchmarks/angular_resolution.py

@@ -10,16 +10,19 @@
 
 
 def angular_resolution(
-    events, energy_bins,
+    events,
+    energy_bins,
     energy_type="true",
-    quantile=ONE_SIGMA_QUANTILE,
+    quantile=[ONE_SIGMA_QUANTILE],
 ):
     """
     Calculate the angular resolution.
 
     This implementation corresponds to the 68% containment of the angular
     distance distribution.
 
+    Passing a list of quantiles results in all the quantiles being calculated.
+
     Parameters
     ----------
     events : astropy.table.QTable
@@ -29,8 +32,8 @@ def angular_resolution(
     energy_type: str
         Either "true" or "reco" energy.
         Default is "true".
-    quantile : float
-        Which quantile to use for the angular resolution,
+    quantile : list(float)
+        Which quantile(s) to use for the angular resolution,
         by default, the containment of the 1-sigma region
         of the normal distribution (~68%) is used.
 
@@ -52,8 +55,17 @@ def angular_resolution(
     result[f"{energy_key}_center"] = 0.5 * (energy_bins[:-1] + energy_bins[1:])
     result["n_events"] = 0
 
-    key = "angular_resolution"
-    result[key] = u.Quantity(np.nan, table["theta"].unit)
+    # keep backwards compatibility
+    if not isinstance(quantile, list) and not isinstance(quantile, np.ndarray):
+        quantile = [quantile]
+
+    if len(quantile) < 2:
+        keys = ["angular_resolution"]
+    else:
+        keys = [f"containment_quantile_{value * 100:.0f}" for value in quantile]
+
+    for key in keys:
+        result[key] = u.Quantity(np.nan, table["theta"].unit)
 
     # if we get an empty input (no selected events available)
     # we return the table filled with NaNs
@@ -63,6 +75,8 @@ def angular_resolution(
     # use groupby operations to calculate the percentile in each bin
     by_bin = table[valid].group_by(bin_index[valid])
     for bin_idx, group in zip(by_bin.groups.keys, by_bin.groups):
-        result[key][bin_idx] = np.nanquantile(group["theta"], quantile)
-        result["n_events"][bin_idx] = len(group)
+        for key, value in zip(keys, quantile):
+            result[key][bin_idx] = np.nanquantile(group["theta"], value)
+            result["n_events"][bin_idx] = len(group)
+
     return result

pyirf/benchmarks/tests/test_angular_resolution.py

@@ -8,18 +8,18 @@
 def test_empty_angular_resolution():
     from pyirf.benchmarks import angular_resolution
 
-    events = QTable({
-        'true_energy': [] * u.TeV,
-        'theta': [] * u.deg,
-    })
-
-    table = angular_resolution(
-        events,
-        [1, 10, 100] * u.TeV
+    events = QTable(
+        {
+            "true_energy": [] * u.TeV,
+            "theta": [] * u.deg,
+        }
     )
 
+    table = angular_resolution(events, [1, 10, 100] * u.TeV)
+
     assert np.all(np.isnan(table["angular_resolution"]))
 
+
 @pytest.mark.parametrize("unit", (u.deg, u.rad))
 def test_angular_resolution(unit):
     from pyirf.benchmarks import angular_resolution
@@ -32,28 +32,32 @@ def test_angular_resolution(unit):
 
     true_resolution = np.append(np.full(N, TRUE_RES_1), np.full(N, TRUE_RES_2))
 
-
     rng = np.random.default_rng(0)
 
-    events = QTable({
-        'true_energy': np.concatenate([
-            [0.5], # below bin 1 to test with underflow
-            np.full(N - 1, 5.0),
-            np.full(N - 1, 50.0),
-            [500], # above bin 2 to test with overflow
-        ]) * u.TeV,
-        'theta': np.abs(rng.normal(0, true_resolution)) * u.deg
-    })
+    events = QTable(
+        {
+            "true_energy": np.concatenate(
+                [
+                    [0.5],  # below bin 1 to test with underflow
+                    np.full(N - 1, 5.0),
+                    np.full(N - 1, 50.0),
+                    [500],  # above bin 2 to test with overflow
+                ]
+            )
+            * u.TeV,
+            "theta": np.abs(rng.normal(0, true_resolution)) * u.deg,
+        }
+    )
 
-    events['theta'] = events['theta'].to(unit)
+    events["theta"] = events["theta"].to(unit)
 
     # add nans to test if nans are ignored
     events["true_energy"].value[N // 2] = np.nan
     events["true_energy"].value[(2 * N) // 2] = np.nan
 
     bins = [1, 10, 100] * u.TeV
     table = angular_resolution(events, bins)
-    ang_res = table['angular_resolution'].to(u.deg)
+    ang_res = table["angular_resolution"].to(u.deg)
     assert len(ang_res) == 2
     assert u.isclose(ang_res[0], TRUE_RES_1 * u.deg, rtol=0.05)
     assert u.isclose(ang_res[1], TRUE_RES_2 * u.deg, rtol=0.05)
@@ -64,9 +68,26 @@ def test_angular_resolution(unit):
     # 2 sigma coverage interval
     quantile = norm(0, 1).cdf(2) - norm(0, 1).cdf(-2)
     table = angular_resolution(events, bins, quantile=quantile)
-    ang_res = table['angular_resolution'].to(u.deg)
-
+    ang_res = table["angular_resolution"].to(u.deg)
     assert len(ang_res) == 2
-
     assert u.isclose(ang_res[0], 2 * TRUE_RES_1 * u.deg, rtol=0.05)
     assert u.isclose(ang_res[1], 2 * TRUE_RES_2 * u.deg, rtol=0.05)
+
+    # 25%, 50%, 90% coverage interval
+    table = angular_resolution(events, bins, quantile=[0.25, 0.5, 0.9])
+    cov_25 = table["containment_quantile_25"].to(u.deg)
+    assert len(cov_25) == 2
+    assert u.isclose(
+        cov_25[0], norm(0, TRUE_RES_1).interval(0.25)[1] * u.deg, rtol=0.05
+    )
+    assert u.isclose(
+        cov_25[1], norm(0, TRUE_RES_2).interval(0.25)[1] * u.deg, rtol=0.05
+    )
+
+    cov_50 = table["containment_quantile_50"].to(u.deg)
+    assert u.isclose(cov_50[0], norm(0, TRUE_RES_1).interval(0.5)[1] * u.deg, rtol=0.05)
+    assert u.isclose(cov_50[1], norm(0, TRUE_RES_2).interval(0.5)[1] * u.deg, rtol=0.05)
+
+    cov_90 = table["containment_quantile_90"].to(u.deg)
+    assert u.isclose(cov_90[0], norm(0, TRUE_RES_1).interval(0.9)[1] * u.deg, rtol=0.05)
+    assert u.isclose(cov_90[1], norm(0, TRUE_RES_2).interval(0.9)[1] * u.deg, rtol=0.05)