feat: Poisson confidence intervals for data uncertainty in data/MC pl…

…ots (#466)

* use Poisson confidence intervals to draw uncertainty for data points in data/MC plots
* implementation uses "Garwood" intervals via hist.intervals.poisson_interval

---------

Co-authored-by: Alexander Held <[email protected]>

.pre-commit-config.yaml

@@ -9,13 +9,13 @@ repos:
     -   id: mypy
         name: mypy with Python 3.9
         files: src/cabinetry
-        additional_dependencies: ["numpy>=1.22", "boost-histogram>=1.0.1", "click>=8", "types-tabulate", "types-PyYAML"]
+        additional_dependencies: ["numpy>=1.22", "boost-histogram>=1.0.1", "click>=8", "types-tabulate", "types-PyYAML", "hist>=2.3.0"]
         # numpy 1.25 is no longer compatible with Python 3.8, so use Python >=3.9 for type checking
         args: ["--python-version=3.9"]
     -   id: mypy
         name: mypy with Python 3.11
         files: src/cabinetry
-        additional_dependencies: ["numpy>=1.22", "boost-histogram>=1.0.1", "click>=8", "types-tabulate", "types-PyYAML"]
+        additional_dependencies: ["numpy>=1.22", "boost-histogram>=1.0.1", "click>=8", "types-tabulate", "types-PyYAML", "hist>=2.3.0"]
         args: ["--python-version=3.11"]
 -   repo: https://github.com/pycqa/flake8
     rev: 7.0.0

setup.cfg

@@ -25,6 +25,7 @@ python_requires = >=3.8
 install_requires =
     pyhf[minuit]~=0.7.0  # model.config.suggested_fixed / .par_names API changes, set_poi(None)
     boost_histogram>=1.0.0  # subclassing with family, 1.02 for stdev scaling fix (currently not needed)
+    hist>=2.5.0  # hist.intervals.poisson_interval
     tabulate>=0.8.1  # multiline text
     matplotlib>=3.5.0  # layout kwarg for subplots
     # below are direct dependencies of cabinetry, which are also included via pyhf[iminuit]

src/cabinetry/visualize/plot_model.py

@@ -4,6 +4,7 @@
 import pathlib
 from typing import Any, Dict, List, Optional
 
+import hist.intervals
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
@@ -35,6 +36,9 @@ def data_mc(
 ) -> mpl.figure.Figure:
     """Draws a data/MC histogram with uncertainty bands and ratio panel.
 
+    Uncertainties for data points are drawn using the "Garwood" frequentist coverage
+    interval as provided by ``hist`` via ``hist.intervals.poisson_interval``.
+
     Args:
         histogram_dict_list (List[Dict[str, Any]]): list of samples (with info stored in
             one dict per sample)
@@ -66,7 +70,10 @@ def data_mc(
     for h in histogram_dict_list:
         if h["isData"]:
             data_histogram_yields = h["yields"]
-            data_histogram_stdev = np.sqrt(data_histogram_yields)
+            # frequentist coverage interval
+            data_histogram_interval = hist.intervals.poisson_interval(
+                np.asarray(data_histogram_yields)
+            )
             data_label = h["label"]
         else:
             mc_histograms_yields.append(h["yields"])
@@ -138,7 +145,7 @@ def data_mc(
     data_container = ax1.errorbar(
         bin_centers_data,
         data_histogram_yields,
-        yerr=data_histogram_stdev,
+        yerr=np.abs(data_histogram_yields - data_histogram_interval),
         fmt="o",
         color="k",
     )
@@ -180,12 +187,14 @@ def data_mc(
 
     # data in ratio plot
     data_model_ratio = data_histogram_yields / total_yield
-    data_model_ratio_unc = data_histogram_stdev / total_yield
+    data_model_ratio_unc = (
+        np.abs(data_histogram_yields - data_histogram_interval) / total_yield
+    )
     # mask data in bins where total model yield is 0
     ax2.errorbar(
         bin_centers_data[nonzero_model_yield],
         data_model_ratio[nonzero_model_yield],
-        yerr=data_model_ratio_unc[nonzero_model_yield],
+        yerr=data_model_ratio_unc[:, nonzero_model_yield],
         fmt="o",
         color="k",
     )