Logo ROOT  
Reference Guide
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Loading...
Searching...
No Matches
HybridStandardForm.py
2# \ingroup tutorial_roostats
3# \notebook -js
4# A hypothesis testing example based on number counting with background uncertainty.
5#
6# A hypothesis testing example based on number counting
7# with background uncertainty.
8#
9# NOTE: This example is like HybridInstructional, but the model is more clearly
10# generalizable to an analysis with shapes. There is a lot of flexibility
11# for how one models a problem in RooFit/RooStats. Models come in a few
12# common forms:
13# - standard form: extended PDF of some discriminating variable m:
14# eg: P(m) ~ S*fs(m) + B*fb(m), with S+B events expected
15# in this case the dataset has N rows corresponding to N events
16# and the extended term is Pois(N|S+B)
17#
18# - fractional form: non-extended PDF of some discriminating variable m:
19# eg: P(m) ~ s*fs(m) + (1-s)*fb(m), where s is a signal fraction
20# in this case the dataset has N rows corresponding to N events
21# and there is no extended term
22#
23# - number counting form: in which there is no discriminating variable
24# and the counts are modeled directly (see HybridInstructional)
25# eg: P(N) = Pois(N|S+B)
26# in this case the dataset has 1 row corresponding to N events
27# and the extended term is the PDF itself.
28#
29# Here we convert the number counting form into the standard form by
30# introducing a dummy discriminating variable m with a uniform distribution.
31#
32# This example:
33# - demonstrates the usage of the HybridCalcultor (Part 4-6)
34# - demonstrates the numerical integration of RooFit (Part 2)
35# - validates the RooStats against an example with a known analytic answer
36# - demonstrates usage of different test statistics
37# - explains subtle choices in the prior used for hybrid methods
38# - demonstrates usage of different priors for the nuisance parameters
39#
40# The basic setup here is that a main measurement has observed x events with an
41# expectation of s+b. One can choose an ad hoc prior for the uncertainty on b,
42# or try to base it on an auxiliary measurement. In this case, the auxiliary
43# measurement (aka control measurement, sideband) is another counting experiment
44# with measurement y and expectation tau*b. With an 'original prior' on b,
45# called \f$ \eta(b) \f$ then one can obtain a posterior from the auxiliary measurement
46# \f$ \pi(b) = \eta(b) * Pois(y|tau*b) \f$. This is a principled choice for a prior
47# on b in the main measurement of x, which can then be treated in a hybrid
48# Bayesian/Frequentist way. Additionally, one can try to treat the two
49# measurements simultaneously, which is detailed in Part 6 of the tutorial.
50#
51# This tutorial is related to the FourBin.C tutorial in the modeling, but
52# focuses on hypothesis testing instead of interval estimation.
53#
54# More background on this 'prototype problem' can be found in the
55# following papers:
56#
57# - Evaluation of three methods for calculating statistical significance
58# when incorporating a systematic uncertainty into a test of the
59# background-only hypothesis for a Poisson process
60# Authors: Robert D. Cousins, James T. Linnemann, Jordan Tucker
61# http:#arxiv.org/abs/physics/0702156
62# NIM A 595 (2008) 480--501
63#
64# - Statistical Challenges for Searches for New Physics at the LHC
65# Author: Kyle Cranmer
66# http:#arxiv.org/abs/physics/0511028
67#
68# - Measures of Significance in HEP and Astrophysics
69# Author: J. T. Linnemann
70# http:#arxiv.org/abs/physics/0312059
71#
72# \macro_image
73# \macro_output
74# \macro_code
75#
76# \authors Kyle Cranmer, Wouter Verkerke, Sven Kreiss, and Jolly Chen (Python translation)
77
78import ROOT
79
80# User configuration parameters
81ntoys = 6000
82nToysRatio = 20 # ratio Ntoys Null/ntoys ALT
83
84# -------------------------------------------------------
85# A New Test Statistic Class for this example.
86# It simply returns the sum of the values in a particular
87# column of a dataset.
88# You can ignore this class and focus on the macro below
89
91 """
92using namespace RooFit;
93using namespace RooStats;
94
95class BinCountTestStat : public TestStatistic {
96public:
97 BinCountTestStat(void) : fColumnName("tmp") {}
98 BinCountTestStat(string columnName) : fColumnName(columnName) {}
99
100 virtual Double_t Evaluate(RooAbsData &data, RooArgSet & /*nullPOI*/)
101 {
102 // This is the main method in the interface
103 Double_t value = 0.0;
104 for (int i = 0; i < data.numEntries(); i++) {
105 value += data.get(i)->getRealValue(fColumnName.c_str());
106 }
107 return value;
108 }
109 virtual const TString GetVarName() const { return fColumnName; }
110
111private:
112 string fColumnName;
113
114protected:
115 ClassDef(BinCountTestStat, 1)
116};
117"""
118)
119
120
121# -----------------------------
122# The Actual Tutorial Macro
123# -----------------------------
124
125# This tutorial has 6 parts
126# Table of Contents
127# Setup
128# 1. Make the model for the 'prototype problem'
129# Special cases
130# 2. NOT RELEVANT HERE
131# 3. Use RooStats analytic solution for this problem
132# RooStats HybridCalculator -- can be generalized
133# 4. RooStats ToyMC version of 2. & 3.
134# 5. RooStats ToyMC with an equivalent test statistic
135# 6. RooStats ToyMC with simultaneous control & main measurement
136
137t = ROOT.TStopwatch()
138t.Start()
139c = ROOT.TCanvas()
140c.Divide(2, 2)
141
142# -----------------------------------------------------
143# P A R T 1 : D I R E C T I N T E G R A T I O N
144# ====================================================
145# Make model for prototype on/off problem
146# Pois(x | s+b) * Pois(y | tau b )
147# for Z_Gamma, use uniform prior on b.
148w = ROOT.RooWorkspace("w")
149
150# replace the pdf in 'number counting form'
151# w.factory("Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))")
152# with one in standard form. Now x is encoded in event count
153w.factory("Uniform::f(m[0,1])") # m is a dummy discriminating variable
154w.factory("ExtendPdf::px(f,sum::splusb(s[0,0,100],b[100,0.1,300]))")
155w.factory("Poisson::py(y[100,0.1,500],prod::taub(tau[1.],b))")
156w.factory("PROD::model(px,py)")
157w.factory("Uniform::prior_b(b)")
158
159# We will control the output level in a few places to avoid
160# verbose progress messages. We start by keeping track
161# of the current threshold on messages.
162msglevel = ROOT.RooMsgService.instance().globalKillBelow()
163
164# -----------------------------------------------
165# P A R T 3 : A N A L Y T I C R E S U L T
166# ==============================================
167# In this special case, the integrals are known analytically
168# and they are implemented in RooStats::NumberCountingUtils
169
170# analytic Z_Bi
173print("-----------------------------------------")
174print("Part 3")
175print(f"Z_Bi p-value (analytic): {p_Bi}")
176print(f"Z_Bi significance (analytic) {Z_Bi}")
177t.Stop()
178t.Print()
179t.Reset()
180t.Start()
181
182# --------------------------------------------------------------
183# P A R T 4 : U S I N G H Y B R I D C A L C U L A T O R
184# ==============================================================
185# Now we demonstrate the RooStats HybridCalculator.
186#
187# Like all RooStats calculators it needs the data and a ModelConfig
188# for the relevant hypotheses. Since we are doing hypothesis testing
189# we need a ModelConfig for the null (background only) and the alternate
190# (signal+background) hypotheses. We also need to specify the PDF,
191# the parameters of interest, and the observables. Furthermore, since
192# the parameter of interest is floating, we need to specify which values
193# of the parameter corresponds to the null and alternate (eg. s=0 and s=50)
194#
195# define some sets of variables obs={x} and poi={s}
196# note here, x is the only observable in the main measurement
197# and y is treated as a separate measurement, which is used
198# to produce the prior that will be used in this calculation
199# to randomize the nuisance parameters.
200w.defineSet("obs", "m")
201w.defineSet("poi", "s")
202
203# create a toy dataset with the x=150
204# data = ROOT.RooDataSet("d", "d", w.set("obs"))
205# data.add(w.set("obs"))
206data = w.pdf("px").generate(w.set("obs"), 150)
207
208# Part 3a : Setup ModelConfigs
209# -------------------------------------------------------
210# create the null (background-only) ModelConfig with s=0
211b_model = ROOT.RooStats.ModelConfig("B_model", w)
215w.var("s").setVal(0.0) # important!
217
218# create the alternate (signal+background) ModelConfig with s=50
219sb_model = ROOT.RooStats.ModelConfig("S+B_model", w)
223w.var("s").setVal(50.0) # important!
225
226# Part 3b : Choose Test Statistic
227# --------------------------------------------------------------
228# To make an equivalent calculation we need to use x as the test
229# statistic. This is not a built-in test statistic in RooStats
230# so we define it above. The new class inherits from the
231# RooStats.TestStatistic interface, and simply returns the value
232# of x in the dataset.
233
234eventCount = ROOT.RooStats.NumEventsTestStat(w.pdf("px"))
235
236# Part 3c : Define Prior used to randomize nuisance parameters
237# -------------------------------------------------------------
238#
239# The prior used for the hybrid calculator is the posterior
240# from the auxiliary measurement y. The model for the aux.
241# measurement is Pois(y|tau*b), thus the likelihood function
242# is proportional to (has the form of) a Gamma distribution.
243# if the 'original prior' $\eta(b)$ is uniform, then from
244# Bayes's theorem we have the posterior:
245# $\pi(b) = Pois(y|tau*b) * \eta(b)$
246# If $\eta(b)$ is flat, then we arrive at a Gamma distribution.
247# Since RooFit will normalize the PDF we can actually supply
248# py=Pois(y,tau*b) that will be equivalent to multiplying by a Uniform.
249#
250# Alternatively, we could explicitly use a gamma distribution:
251#
252# `w.factory("Gamma::gamma(b,sum::temp(y,1),1,0)")`
253#
254# or we can use some other ad hoc prior that do not naturally
255# follow from the known form of the auxiliary measurement.
256# The common choice is the equivalent Gaussian:
257w.factory("Gaussian::gauss_prior(b,y, expr::sqrty('sqrt(y)',y))")
258# this corresponds to the "Z_N" calculation.
259#
260# or one could use the analogous log-normal prior
261w.factory("Lognormal::lognorm_prior(b,y, expr::kappa('1+1./sqrt(y)',y))")
262#
263# Ideally, the HybridCalculator would be able to inspect the full
264# model Pois(x | s+b) * Pois(y | tau b ) and be given the original
265# prior $\eta(b)$ to form $\pi(b) = Pois(y|tau*b) * \eta(b)$.
266# This is not yet implemented because in the general case
267# it is not easy to identify the terms in the PDF that correspond
268# to the auxiliary measurement. So for now, it must be set
269# explicitly with:
270# - ForcePriorNuisanceNull()
271# - ForcePriorNuisanceAlt()
272# the name "ForcePriorNuisance" was chosen because we anticipate
273# this to be auto-detected, but will leave the option open
274# to force to a different prior for the nuisance parameters.
275
276# Part 3d : Construct and configure the HybridCalculator
277# -------------------------------------------------------
278
279hc1 = ROOT.RooStats.HybridCalculator(data, sb_model, b_model)
280toymcs1 = hc1.GetTestStatSampler()
281# toymcs1.SetNEventsPerToy(1) # because the model is in number counting form
282toymcs1.SetTestStatistic(eventCount) # set the test statistic
283# toymcs1.SetGenerateBinned()
284hc1.SetToys(ntoys, ntoys // nToysRatio)
287# if you wanted to use the ad hoc Gaussian prior instead
288# ~~~
289# hc1.ForcePriorNuisanceAlt(w.pdf("gauss_prior"))
290# hc1.ForcePriorNuisanceNull(w.pdf("gauss_prior"))
291# ~~~
292# if you wanted to use the ad hoc log-normal prior instead
293# ~~~
294# hc1.ForcePriorNuisanceAlt(w.pdf("lognorm_prior"))
295# hc1.ForcePriorNuisanceNull(w.pdf("lognorm_prior"))
296# ~~~
297
298# these lines save current msg level and then kill any messages below ERROR
300# Get the result
301r1 = hc1.GetHypoTest()
302ROOT.RooMsgService.instance().setGlobalKillBelow(msglevel) # set it back
303print("-----------------------------------------")
304print("Part 4")
305r1.Print()
306t.Stop()
307t.Print()
308t.Reset()
309t.Start()
310
311c.cd(2)
312p1 = ROOT.RooStats.HypoTestPlot(r1, 30) # 30 bins, TS is discrete
313p1.Draw()
314
315quit() # keep the running time short by default
316# -------------------------------------------------------------------------
317# # P A R T 5 : U S I N G H Y B R I D C A L C U L A T O R W I T H
318# # A N A L T E R N A T I V E T E S T S T A T I S T I C
319#
320# A likelihood ratio test statistics should be 1-to-1 with the count x
321# when the value of b is fixed in the likelihood. This is implemented
322# by the SimpleLikelihoodRatioTestStat
323
327
328# HYBRID CALCULATOR
329hc2 = ROOT.RooStats.HybridCalculator(data, sb_model, b_model)
331toymcs2 = hc2.GetTestStatSampler()
332# toymcs2.SetNEventsPerToy(1)
334# toymcs2.SetGenerateBinned()
335hc2.SetToys(ntoys, ntoys // nToysRatio)
338# if you wanted to use the ad hoc Gaussian prior instead
339# ~~~
340# hc2.ForcePriorNuisanceAlt(w.pdf("gauss_prior"))
341# hc2.ForcePriorNuisanceNull(w.pdf("gauss_prior"))
342# ~~~
343# if you wanted to use the ad hoc log-normal prior instead
344# ~~~
345# hc2.ForcePriorNuisanceAlt(w.pdf("lognorm_prior"))
346# hc2.ForcePriorNuisanceNull(w.pdf("lognorm_prior"))
347# ~~~
348
349# these lines save current msg level and then kill any messages below ERROR
351# Get the result
352r2 = hc2.GetHypoTest()
353print("-----------------------------------------")
354print("Part 5")
355r2.Print()
356t.Stop()
357t.Print()
358t.Reset()
359t.Start()
360ROOT.RooMsgService.instance().setGlobalKillBelow(msglevel)
361
362c.cd(3)
363p2 = ROOT.RooStats.HypoTestPlot(r2, 30) # 30 bins
364p2.Draw()
365
366quit() # so standard tutorial runs faster
367
368# ---------------------------------------------
369# OUTPUT W/O PROOF (2.66 GHz Intel Core i7)
370# ============================================
371
372# -----------------------------------------
373# Part 3
374# Z_Bi p-value (analytic): 0.00094165
375# Z_Bi significance (analytic): 3.10804
376# Real time 0:00:00, CP time 0.610
377
378# Results HybridCalculator_result:
379# - Null p-value = 0.00103333 +/- 0.000179406
380# - Significance = 3.08048 sigma
381# - Number of S+B toys: 1000
382# - Number of B toys: 30000
383# - Test statistic evaluated on data: 150
384# - CL_b: 0.998967 +/- 0.000185496
385# - CL_s+b: 0.495 +/- 0.0158106
386# - CL_s: 0.495512 +/- 0.0158272
387# Real time 0:04:43, CP time 283.780
388
389# With PROOF
390# -----------------------------------------
391# Part 5
392
393# Results HybridCalculator_result:
394# - Null p-value = 0.00105 +/- 0.000206022
395# - Significance = 3.07571 sigma
396# - Number of S+B toys: 1000
397# - Number of B toys: 20000
398# - Test statistic evaluated on data: 10.8198
399# - CL_b: 0.99895 +/- 0.000229008
400# - CL_s+b: 0.491 +/- 0.0158088
401# - CL_s: 0.491516 +/- 0.0158258
402# Real time 0:02:22, CP time 0.990
403
404# -------------------------------------------------------
405# Comparison
406# -------------------------------------------------------
407# LEPStatToolsForLHC
408# https:#plone4.fnal.gov:4430/P0/phystat/packages/0703002
409# Uses Gaussian prior
410# CL_b = 6.218476e-04, Significance = 3.228665 sigma
411#
412# -------------------------------------------------------
413# Comparison
414# -------------------------------------------------------
415# Asymptotic
416# From the value of the profile likelihood ratio (5.0338)
417# The significance can be estimated using Wilks's theorem
418# significance = sqrt(2*profileLR) = 3.1729 sigma
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.