fix: Clopper-Pearson and poisson ratio intervals (#281)

* Fix bug in clopper pearson interval

* Implement correct poisson-ratio interval

Fixes #279

* spelling

* Update tests and fix small issue in poisson interval

* Remove 'magic' from poisson_interval, clarify docstrings

* black

* Add changelog

* docs: Update docstring hyperlinks

* Have hist functions link to the respective docstrings
* Use DOIs or names of work where possible

Co-authored-by: Matthew Feickert <matthew.feickert@cern.ch>

docs/changelog.rst

@@ -11,6 +11,9 @@ Version 2.5.0
   `#266 <https://github.com/scikit-hep/hist/pull/266>`_
   `#278 <https://github.com/scikit-hep/hist/pull/278>`_
 
+* Improve and clarify treatment of confidence intervals in ``intervals`` submodule.
+  `#281 <https://github.com/scikit-hep/hist/pull/281>`_
+
 
 Version 2.4.0
 --------------------

src/hist/intervals.py

@@ -25,22 +25,26 @@ def __dir__() -> tuple[str, ...]:
 
 
 def poisson_interval(
-    values: np.ndarray, variances: np.ndarray, coverage: float | None = None
+    values: np.ndarray,
+    variances: np.ndarray | None = None,
+    coverage: float | None = None,
 ) -> np.ndarray:
     r"""
     The Frequentist coverage interval for Poisson-distributed observations.
 
-    What is calculated is the "Garwood" interval,
-    c.f. https://www.ine.pt/revstat/pdf/rs120203.pdf or
-    http://ms.mcmaster.ca/peter/s743/poissonalpha.html.
-    For weighted data, this approximates the observed count by
-    ``values**2/variances``, which effectively scales the unweighted
-    Poisson interval by the average weight.
-    This may not be the optimal solution: see https://arxiv.org/abs/1309.1287
-    for a proper treatment.
+    What is calculated is the "Garwood" interval, c.f.
+    `V. Patil, H. Kulkarni (Revstat, 2012) <https://www.ine.pt/revstat/pdf/rs120203.pdf>`_
+    or http://ms.mcmaster.ca/peter/s743/poissonalpha.html.
+    If ``variances`` is supplied, the data is assumed to be weighted, and the
+    unweighted count is approximated by ``values**2/variances``, which effectively
+    scales the unweighted Poisson interval by the average weight.
+    This may not be the optimal solution: see
+    `10.1016/j.nima.2014.02.021 <https://doi.org/10.1016/j.nima.2014.02.021>`_
+    (`arXiv:1309.1287 <https://arxiv.org/abs/1309.1287>`_) for a proper treatment.
 
-    When a bin is zero, the scale of the nearest nonzero bin is substituted to
-    scale the nominal upper bound.
+    In cases where the value is zero, an upper limit is well-defined only in the case of
+    unweighted data, so if ``variances`` is supplied, the upper limit for a zero value
+    will be set to ``NaN``.
 
     Args:
         values: Sum of weights.
@@ -55,25 +59,22 @@ def poisson_interval(
     # https://github.com/CoffeaTeam/coffea/blob/8c58807e199a7694bf15e3803dbaf706d34bbfa0/LICENSE
     if coverage is None:
         coverage = stats.norm.cdf(1) - stats.norm.cdf(-1)
-    scale = np.empty_like(values)
-    scale[values != 0] = variances[values != 0] / values[values != 0]
-    if np.sum(values == 0) > 0:
-        missing = np.where(values == 0)
-        available = np.nonzero(values)
-        if len(available[0]) == 0:
-            raise RuntimeError(
-                "All values are zero! Cannot compute meaningful uncertainties.",
-            )
-        nearest = np.sum(
-            [np.square(np.subtract.outer(d, d0)) for d, d0 in zip(available, missing)]
-        ).argmin(axis=0)
-        argnearest = tuple(dim[nearest] for dim in available)
-        scale[missing] = scale[argnearest]
-    counts = values / scale
-    interval_min = scale * stats.chi2.ppf((1 - coverage) / 2, 2 * counts) / 2.0
-    interval_max = scale * stats.chi2.ppf((1 + coverage) / 2, 2 * (counts + 1)) / 2.0
+    if variances is None:
+        interval_min = stats.chi2.ppf((1 - coverage) / 2, 2 * values) / 2.0
+        interval_min[values == 0.0] = 0.0  # chi2.ppf produces NaN for values=0
+        interval_max = stats.chi2.ppf((1 + coverage) / 2, 2 * (values + 1)) / 2.0
+    else:
+        scale = np.ones_like(values)
+        mask = np.isfinite(values) & (values != 0)
+        np.divide(variances, values, out=scale, where=mask)
+        counts = values / scale
+        interval_min = scale * stats.chi2.ppf((1 - coverage) / 2, 2 * counts) / 2.0
+        interval_min[values == 0.0] = 0.0  # chi2.ppf produces NaN for values=0
+        interval_max = (
+            scale * stats.chi2.ppf((1 + coverage) / 2, 2 * (counts + 1)) / 2.0
+        )
+        interval_max[values == 0.0] = np.nan
     interval = np.stack((interval_min, interval_max))
-    interval[interval == np.nan] = 0.0  # chi2.ppf produces nan for counts=0
     return interval
 
 
@@ -105,7 +106,7 @@ def clopper_pearson_interval(
     interval_min = stats.beta.ppf((1 - coverage) / 2, num, denom - num + 1)
     interval_max = stats.beta.ppf((1 + coverage) / 2, num + 1, denom - num)
     interval = np.stack((interval_min, interval_max))
-    interval[:, num == 0.0] = 0.0
+    interval[0, num == 0.0] = 0.0
     interval[1, num == denom] = 1.0
     return interval
 
@@ -124,26 +125,38 @@ def ratio_uncertainty(
         denom: Denominator or number of trials.
         uncertainty_type: Coverage interval type to use in the calculation of
          the uncertainties.
-         ``"poisson"`` (default) implements the Poisson interval for the
-         numerator scaled by the denominator.
-         ``"poisson-ratio"`` implements the Clopper-Pearson interval for Poisson
-         distributed ``num`` and ``denom`` where it is assumed that ``num`` and
-         ``denom`` are independent.
-         ``"efficiency"`` implements the Clopper-Pearson interval for Poisson
-         distributed ``num`` and ``denom``, with ``num`` assumed to be a
-         strict subset of ``denom``.
+
+         * ``"poisson"`` (default) implements the Garwood confidence interval for
+           a Poisson-distributed numerator scaled by the denominator.
+           See :func:`hist.intervals.poisson_interval` for further details.
+         * ``"poisson-ratio"`` implements a confidence interval for the ratio ``num / denom``
+           assuming it is an estimator of the ratio of the expected rates from
+           two independent Poisson distributions.
+           It over-covers to a similar degree as the Clopper-Pearson interval
+           does for the Binomial efficiency parameter estimate.
+         * ``"efficiency"`` implements the Clopper-Pearson confidence interval
+           for the ratio ``num / denom`` assuming it is an estimator of a Binomial
+           efficiency parameter.
+           This is only valid if the entries contributing to ``num`` are a strict
+           subset of those contributing to ``denom``.
 
     Returns:
         The uncertainties for the ratio.
     """
     # Note: As return is a numpy ufuncs the type is "Any"
-    with np.errstate(divide="ignore"):
+    with np.errstate(divide="ignore", invalid="ignore"):
+        # Nota bene: x/0 = inf, 0/0 = nan
         ratio = num / denom
     if uncertainty_type == "poisson":
-        ratio_uncert = np.abs(poisson_interval(ratio, num / np.square(denom)) - ratio)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            ratio_variance = num * np.power(denom, -2.0)
+        ratio_uncert = np.abs(poisson_interval(ratio, ratio_variance) - ratio)
     elif uncertainty_type == "poisson-ratio":
-        # poisson ratio n/m is equivalent to binomial n/(n+m)
-        ratio_uncert = np.abs(clopper_pearson_interval(num, num + denom) - ratio)
+        # Details: see https://github.com/scikit-hep/hist/issues/279
+        p_lim = clopper_pearson_interval(num, num + denom)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            r_lim = p_lim / (1 - p_lim)
+            ratio_uncert = np.abs(r_lim - ratio)
     elif uncertainty_type == "efficiency":
         ratio_uncert = np.abs(clopper_pearson_interval(num, denom) - ratio)
     else:

tests/test_intervals.py

@@ -61,6 +61,12 @@ def test_poisson_interval(hist_fixture):
         ]
     )
 
+    interval_min, interval_max = intervals.poisson_interval(np.arange(4))
+    assert approx(interval_min) == np.array([0.0, 0.17275378, 0.70818544, 1.36729531])
+    assert approx(interval_max) == np.array(
+        [1.84102165, 3.29952656, 4.63785962, 5.91818583]
+    )
+
 
 def test_clopper_pearson_interval(hist_fixture):
     hist_1, _ = hist_fixture
@@ -100,47 +106,79 @@ def test_clopper_pearson_interval(hist_fixture):
     )
 
 
-def test_ratio_uncertainty(hist_fixture):
-    hist_1, hist_2 = hist_fixture
+def test_ratio_uncertainty():
+    num, denom = np.meshgrid(np.array([0, 1, 4, 512]), np.array([0, 1, 4, 512]))
 
     uncertainty_min, uncertainty_max = intervals.ratio_uncertainty(
-        hist_1.values(), hist_2.values(), uncertainty_type="poisson"
+        num, denom, uncertainty_type="poisson"
+    )
+
+    assert approx(uncertainty_min, nan_ok=True) == np.array(
+        [
+            [np.nan, np.nan, np.nan, np.nan],
+            [0.0, 8.27246221e-01, 1.91433919e00, 2.26200365e01],
+            [0.0, 2.06811555e-01, 4.78584797e-01, 5.65500911e00],
+            [0.0, 1.61571528e-03, 3.73894372e-03, 4.41797587e-02],
+        ]
     )
 
-    assert approx(uncertainty_min) == np.array(
+    assert approx(uncertainty_max, nan_ok=True) == np.array(
         [
-            0.1439794096271186,
-            0.12988019998066708,
-            0.0711565635066328,
-            0.045722288708959336,
-            0.04049103990124614,
-            0.038474711321686006,
-            0.045227104349518155,
-            0.06135954973309016,
-            0.12378460125991042,
-            0.19774186117590858,
+            [np.nan, np.nan, np.nan, np.nan],
+            [np.nan, 2.29952656e00, 3.16275317e00, 2.36421589e01],
+            [np.nan, 5.74881640e-01, 7.90688293e-01, 5.91053972e00],
+            [np.nan, 4.49126281e-03, 6.17725229e-03, 4.61760915e-02],
         ]
     )
 
-    assert approx(uncertainty_max) == np.array(
+    uncertainty_min, uncertainty_max = intervals.ratio_uncertainty(
+        num, denom, uncertainty_type="poisson-ratio"
+    )
+
+    assert approx(uncertainty_min, nan_ok=True) == np.array(
         [
-            0.22549817680979262,
-            0.1615766277480729,
-            0.07946632561746425,
-            0.04954668134626106,
-            0.04327624938437291,
-            0.04106267733757407,
-            0.04891233040201837,
-            0.06909296140898324,
-            0.1485919630151803,
-            0.2817958228477908,
+            [np.nan, np.inf, np.inf, np.inf],
+            [0.0, 9.09782858e-01, 3.09251539e00, 3.57174304e02],
+            [0.0, 2.14845433e-01, 6.11992834e-01, 5.67393184e01],
+            [0.0, 1.61631629e-03, 3.75049626e-03, 6.24104251e-02],
+        ]
+    )
+
+    assert approx(uncertainty_max, nan_ok=True) == np.array(
+        [
+            [np.nan, np.nan, np.nan, np.nan],
+            [5.30297438e00, 1.00843679e01, 2.44458061e01, 2.45704433e03],
+            [5.84478627e-01, 8.51947064e-01, 1.57727199e00, 1.18183919e02],
+            [3.60221785e-03, 4.50575120e-03, 6.22048393e-03, 6.65647601e-02],
+        ]
+    )
+
+    with pytest.raises(ValueError):
+        intervals.ratio_uncertainty(num, denom, uncertainty_type="efficiency")
+
+    uncertainty_min, uncertainty_max = intervals.ratio_uncertainty(
+        np.minimum(num, denom), denom, uncertainty_type="efficiency"
+    )
+
+    assert approx(uncertainty_min, nan_ok=True) == np.array(
+        [
+            [np.nan, np.nan, np.nan, np.nan],
+            [0.0, 0.8413447460685429, 0.8413447460685429, 0.8413447460685429],
+            [0.0, 0.207730893696323, 0.36887757085042716, 0.36887757085042716],
+            [0.0, 0.0016157721916044239, 0.003735294987003171, 0.0035892884494188593],
+        ]
+    )
+    assert approx(uncertainty_max, nan_ok=True) == np.array(
+        [
+            [np.nan, np.nan, np.nan, np.nan],
+            [0.8413447460685429, 0.0, 0.0, 0.0],
+            [0.3688775708504272, 0.36840242550395996, 0.0, 0.0],
+            [0.0035892884494188337, 0.004476807721636625, 0.006134065381665161, 0.0],
         ]
     )
 
     with pytest.raises(TypeError):
-        intervals.ratio_uncertainty(
-            hist_1.values(), hist_2.values(), uncertainty_type="fail"
-        )
+        intervals.ratio_uncertainty(num, denom, uncertainty_type="fail")
 
 
 def test_valid_efficiency_ratio_uncertainty(hist_fixture):