From 77389e8983d7c8e3f35ced011d762b2737cd7f86 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 17 Oct 2019 11:07:27 +0200 Subject: [PATCH 1/5] add normal interrval to plotratio --- coffea/hist/plot.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/coffea/hist/plot.py b/coffea/hist/plot.py index ae8a428c1..0c72375ed 100644 --- a/coffea/hist/plot.py +++ b/coffea/hist/plot.py @@ -77,6 +77,23 @@ def clopper_pearson_interval(num, denom, coverage=_coverage1sd): interval[1, num == denom] = 1. 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-_coverage1sd) + 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): @@ -352,6 +369,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) From c6971fad5eff4ae25511c286eea16984d3026205 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 17 Oct 2019 16:27:14 +0200 Subject: [PATCH 2/5] flake8 fixes --- coffea/hist/plot.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/coffea/hist/plot.py b/coffea/hist/plot.py index 0c72375ed..fce8119f3 100644 --- a/coffea/hist/plot.py +++ b/coffea/hist/plot.py @@ -77,15 +77,16 @@ def clopper_pearson_interval(num, denom, coverage=_coverage1sd): interval[1, num == denom] = 1. 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 + eff = pw / tw - variance = ( pw2 * (1 - 2*eff) + tw2 * eff**2 ) / (tw**2) + variance = (pw2 * (1 - 2 * eff) + tw2 * eff**2) / (tw**2) sigma = np.sqrt(variance) - prob = 0.5*(1-_coverage1sd) + prob = 0.5 * (1 - _coverage1sd) delta = np.zeros_like(sigma) delta[sigma != 0] = -scipy.stats.norm.ppf(prob, scale=sigma[sigma != 0]) From 217ed3ea3451a53e0b4d87243743e6baa3bd7963 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 18 Oct 2019 10:13:18 +0200 Subject: [PATCH 3/5] s/_coverage1sd/coverage/g --- coffea/hist/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coffea/hist/plot.py b/coffea/hist/plot.py index fce8119f3..b0e12ba02 100644 --- a/coffea/hist/plot.py +++ b/coffea/hist/plot.py @@ -86,7 +86,7 @@ def normal_interval(pw, tw, pw2, tw2, coverage=_coverage1sd): variance = (pw2 * (1 - 2 * eff) + tw2 * eff**2) / (tw**2) sigma = np.sqrt(variance) - prob = 0.5 * (1 - _coverage1sd) + prob = 0.5 * (1 - coverage) delta = np.zeros_like(sigma) delta[sigma != 0] = -scipy.stats.norm.ppf(prob, scale=sigma[sigma != 0]) From a4fd53a592ad42076bbc10f81d53d8bee6c92d2b Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 18 Oct 2019 13:18:03 +0200 Subject: [PATCH 4/5] flip sign --- coffea/hist/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coffea/hist/plot.py b/coffea/hist/plot.py index b0e12ba02..3f326ad00 100644 --- a/coffea/hist/plot.py +++ b/coffea/hist/plot.py @@ -88,7 +88,7 @@ def normal_interval(pw, tw, pw2, tw2, coverage=_coverage1sd): prob = 0.5 * (1 - coverage) delta = np.zeros_like(sigma) - delta[sigma != 0] = -scipy.stats.norm.ppf(prob, scale=sigma[sigma != 0]) + 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 From a8a0ddc0c1426bd387375259ee86385174ee3125 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 18 Oct 2019 13:18:16 +0200 Subject: [PATCH 5/5] add test for normal_interval --- tests/test_hist_plot.py | 44 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/tests/test_hist_plot.py b/tests/test_hist_plot.py index a2cae2c6b..ce98916a4 100644 --- a/tests/test_hist_plot.py +++ b/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 + + +