diff --git a/docs/changelog.rst b/docs/changelog.rst index a2f849e9..e55bf980 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -11,6 +11,9 @@ Version 2.5.0 `#266 `_ `#278 `_ +* Improve and clarify treatment of confidence intervals in ``intervals`` submodule. + `#281 `_ + Version 2.4.0 -------------------- diff --git a/src/hist/intervals.py b/src/hist/intervals.py index 6d9de11a..41725362 100644 --- a/src/hist/intervals.py +++ b/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) `_ + 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 `_ + (`arXiv: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: diff --git a/tests/test_intervals.py b/tests/test_intervals.py index 0f06dc5f..b2420e15 100644 --- a/tests/test_intervals.py +++ b/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):