Merge pull request #182 from paulgessinger/normal-interval-eff-error

Add normal interval to plotratio

Author: lgray

coffea/hist/plot.py

@@ -78,6 +78,24 @@ def clopper_pearson_interval(num, denom, coverage=_coverage1sd):
     return interval
 
 
+def normal_interval(pw, tw, pw2, tw2, coverage=_coverage1sd):
+    # see https://root.cern.ch/doc/master/TEfficiency_8cxx_source.html#l02515
+
+    eff = pw / tw
+
+    variance = (pw2 * (1 - 2 * eff) + tw2 * eff**2) / (tw**2)
+    sigma = np.sqrt(variance)
+
+    prob = 0.5 * (1 - coverage)
+    delta = np.zeros_like(sigma)
+    delta[sigma != 0] = scipy.stats.norm.ppf(prob, scale=sigma[sigma != 0])
+
+    lo = eff - np.minimum(eff + delta, np.ones_like(eff))
+    hi = np.maximum(eff - delta, np.zeros_like(eff)) - eff
+
+    return np.array([lo, hi])
+
+
 def plot1d(hist, ax=None, clear=True, overlay=None, stack=False, overflow='none', line_opts=None,
            fill_opts=None, error_opts=None, legend_opts={}, overlay_overflow='none', density=False, binwnorm=None):
     """Create a 1D plot from a 1D or 2D `Hist` object
@@ -352,6 +370,8 @@ def plotratio(num, denom, ax=None, clear=True, overflow='none', error_opts=None,
             rsumw_err = np.abs(clopper_pearson_interval(sumw_num, sumw_num + sumw_denom) - rsumw)
         elif unc == 'num':
             rsumw_err = np.abs(poisson_interval(rsumw, sumw2_num / sumw_denom**2) - rsumw)
+        elif unc == "normal":
+            rsumw_err = np.abs(normal_interval(sumw_num, sumw_denom, sumw2_num, sumw2_denom))
         else:
             raise ValueError("Unrecognized uncertainty option: %r" % unc)
 

tests/test_hist_plot.py

@@ -241,3 +241,47 @@ def test_clopper_pearson_interval():
     threshold = 1e-6
     assert(all((interval[1, :] / ref_hi) - 1 < threshold))
     assert(all((interval[0, :] / ref_lo) - 1 < threshold))
+
+def test_normal_interval():
+    from coffea.hist.plot import normal_interval
+
+    # Reference weighted efficiency and error from ROOTs TEfficiency
+
+    denom = np.array([  89.01457591590004, 2177.066076428943  , 6122.5256890981855 ,
+              0.              ,  100.27757990710668])
+    num = np.array([  75.14287743709515, 2177.066076428943  , 5193.454723043864  ,
+              0.              ,   84.97723540536361])
+    denom_sumw2 = np.array([   94.37919737476827, 10000.              ,  6463.46795877633   ,
+               0.              ,   105.90898005417333])
+    num_sumw2 = np.array([   67.2202147680005 , 10000.              ,  4647.983931785646  ,
+               0.              ,    76.01275761253757])
+    ref_hi = np.array([0.0514643476600107, 0.                , 0.0061403263960343,
+                          np.nan, 0.0480731185500146])
+    ref_lo = np.array([0.0514643476600107, 0.                , 0.0061403263960343,
+                          np.nan, 0.0480731185500146])
+
+    interval = normal_interval(num, denom, num_sumw2, denom_sumw2)
+    threshold = 1e-6
+
+    lo, hi = interval
+
+    assert len(ref_hi) == len(hi)
+    assert len(ref_lo) == len(lo)
+
+    for i in range(len(ref_hi)):
+        if np.isnan(ref_hi[i]):
+            assert np.isnan(ref_hi[i])
+        elif ref_hi[i] == 0.0:
+            assert hi[i] == 0.0
+        else:
+            assert np.abs(hi[i] / ref_hi[i] - 1) < threshold
+
+        if np.isnan(ref_lo[i]):
+            assert np.isnan(ref_lo[i])
+        elif ref_lo[i] == 0.0:
+            assert lo[i] == 0.0
+        else:
+            assert np.abs(lo[i] / ref_lo[i] - 1) < threshold
+
+
+