Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TwoSidedFrequentistUpperLimitWithBands.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// TwoSidedFrequentistUpperLimitWithBands
5///
6///
7/// This is a standard demo that can be used with any ROOT file
8/// prepared in the standard way. You specify:
9/// - name for input ROOT file
10/// - name of workspace inside ROOT file that holds model and data
11/// - name of ModelConfig that specifies details for calculator tools
12/// - name of dataset
13///
14/// With default parameters the macro will attempt to run the
15/// standard hist2workspace example and read the ROOT file
16/// that it produces.
17///
18/// You may want to control:
19/// ~~~{.cpp}
20/// double confidenceLevel=0.95;
21/// double additionalToysFac = 1.;
22/// int nPointsToScan = 12;
23/// int nToyMC = 200;
24/// ~~~
25///
26/// This uses a modified version of the profile likelihood ratio as
27/// a test statistic for upper limits (eg. test stat = 0 if muhat>mu).
28///
29/// Based on the observed data, one defines a set of parameter points
30/// to be tested based on the value of the parameter of interest
31/// and the conditional MLE (eg. profiled) values of the nuisance parameters.
32///
33/// At each parameter point, pseudo-experiments are generated using this
34/// fixed reference model and then the test statistic is evaluated.
35/// The auxiliary measurements (global observables) associated with the
36/// constraint terms in nuisance parameters are also fluctuated in the
37/// process of generating the pseudo-experiments in a frequentist manner
38/// forming an 'unconditional ensemble'. One could form a 'conditional'
39/// ensemble in which these auxiliary measurements are fixed. Note that the
40/// nuisance parameters are not randomized, which is a Bayesian procedure.
41/// Note, the nuisance parameters are floating in the fits. For each point,
42/// the threshold that defines the 95% acceptance region is found. This
43/// forms a "Confidence Belt".
44///
45/// After constructing the confidence belt, one can find the confidence
46/// interval for any particular dataset by finding the intersection
47/// of the observed test statistic and the confidence belt. First
48/// this is done on the observed data to get an observed 1-sided upper limt.
49///
50/// Finally, there expected limit and bands (from background-only) are
51/// formed by generating background-only data and finding the upper limit.
52/// The background-only is defined as such that the nuisance parameters are
53/// fixed to their best fit value based on the data with the signal rate fixed to 0.
54/// The bands are done by hand for now, will later be part of the RooStats tools.
55///
56/// On a technical note, this technique IS the generalization of Feldman-Cousins
57/// with nuisance parameters.
58///
59/// Building the confidence belt can be computationally expensive.
60/// Once it is built, one could save it to a file and use it in a separate step.
61///
62/// Note, if you have a boundary on the parameter of interest (eg. cross-section)
63/// the threshold on the two-sided test statistic starts off at moderate values and plateaus.
64///
65/// [#0] PROGRESS:Generation -- generated toys: 500 / 999
66/// NeymanConstruction: Prog: 12/50 total MC = 39 this test stat = 0
67/// SigXsecOverSM=0.69 alpha_syst1=0.136515 alpha_syst3=0.425415 beta_syst2=1.08496 [-1e+30, 0.011215] in interval = 1
68///
69/// this tells you the values of the parameters being used to generate the pseudo-experiments
70/// and the threshold in this case is 0.011215. One would expect for 95% that the threshold
71/// would be ~1.35 once the cross-section is far enough away from 0 that it is essentially
72/// unaffected by the boundary. As one reaches the last points in the scan, the
73/// threshold starts to get artificially high. This is because the range of the parameter in
74/// the fit is the same as the range in the scan. In the future, these should be independently
75/// controlled, but they are not now. As a result the ~50% of pseudo-experiments that have an
76/// upward fluctuation end up with muhat = muMax. Because of this, the upper range of the
77/// parameter should be well above the expected upper limit... but not too high or one will
78/// need a very large value of nPointsToScan to resolve the relevant region. This can be
79/// improved, but this is the first version of this script.
80///
81/// Important note: when the model includes external constraint terms, like a Gaussian
82/// constraint to a nuisance parameter centered around some nominal value there is
83/// a subtlety. The asymptotic results are all based on the assumption that all the
84/// measurements fluctuate... including the nominal values from auxiliary measurements.
85/// If these do not fluctuate, this corresponds to an "conditional ensemble". The
86/// result is that the distribution of the test statistic can become very non-chi^2.
87/// This results in thresholds that become very large.
88///
89/// \macro_image
90/// \macro_output
91/// \macro_code
92///
93/// \authors Kyle Cranmer,Contributions from Aaron Armbruster, Haoshuang Ji, Haichen Wang and Daniel Whiteson
94
95#include "TFile.h"
96#include "TROOT.h"
97#include "TH1F.h"
98#include "TCanvas.h"
99#include "TSystem.h"
100#include <iostream>
101
102#include "RooWorkspace.h"
103#include "RooSimultaneous.h"
104#include "RooAbsData.h"
105
106#include "RooStats/ModelConfig.h"
111
114
115using namespace RooFit;
116using namespace RooStats;
117using std::cout, std::endl;
118
119// -------------------------------------------------------
120
121void TwoSidedFrequentistUpperLimitWithBands(const char *infile = "", const char *workspaceName = "combined",
122 const char *modelConfigName = "ModelConfig",
123 const char *dataName = "obsData")
124{
125
126 double confidenceLevel = 0.95;
127 // degrade/improve number of pseudo-experiments used to define the confidence belt.
128 // value of 1 corresponds to default number of toys in the tail, which is 50/(1-confidenceLevel)
129 double additionalToysFac = 0.5;
130 int nPointsToScan = 20; // number of steps in the parameter of interest
131 int nToyMC = 200; // number of toys used to define the expected limit and band
132
133 // -------------------------------------------------------
134 // First part is just to access a user-defined file
135 // or create the standard example file if it doesn't exist
136 const char *filename = "";
137 if (!strcmp(infile, "")) {
138 filename = "results/example_combined_GaussExample_model.root";
139 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
140 // if file does not exists generate with histfactory
141 if (!fileExist) {
142 // Normally this would be run on the command line
143 cout << "will run standard hist2workspace example" << endl;
144 gROOT->ProcessLine(".! prepareHistFactory .");
145 gROOT->ProcessLine(".! hist2workspace config/example.xml");
146 cout << "\n\n---------------------" << endl;
147 cout << "Done creating example input" << endl;
148 cout << "---------------------\n\n" << endl;
149 }
150
151 } else
153
154 // Try to open the file
156
157 // -------------------------------------------------------
158 // Now get the data and workspace
159
160 // get the workspace out of the file
162
163 // get the modelConfig out of the file
165
166 // get the modelConfig out of the file
167 RooAbsData *data = w->data(dataName);
168
169 cout << "Found data and ModelConfig:" << endl;
170 mc->Print();
171
172 // -------------------------------------------------------
173 // Now get the POI for convenience
174 // you may want to adjust the range of your POI
175 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
176 /* firstPOI->setMin(0);*/
177 /* firstPOI->setMax(10);*/
178
179 // -------------------------------------------------------
180 // create and use the FeldmanCousins tool
181 // to find and plot the 95% confidence interval
182 // on the parameter of interest as specified
183 // in the model config
184 // REMEMBER, we will change the test statistic
185 // so this is NOT a Feldman-Cousins interval
187 fc.SetConfidenceLevel(confidenceLevel);
188 fc.AdditionalNToysFactor(additionalToysFac); // improve sampling that defines confidence belt
189 // fc.UseAdaptiveSampling(true); // speed it up a bit, but don't use for expected limits
190 fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
191 fc.CreateConfBelt(true); // save the information in the belt for plotting
192
193 // -------------------------------------------------------
194 // Feldman-Cousins is a unified limit by definition
195 // but the tool takes care of a few things for us like which values
196 // of the nuisance parameters should be used to generate toys.
197 // so let's just change the test statistic and realize this is
198 // no longer "Feldman-Cousins" but is a fully frequentist Neyman-Construction.
199 // fc.GetTestStatSampler()->SetTestStatistic(&onesided);
200 // ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(true);
201 ToyMCSampler *toymcsampler = (ToyMCSampler *)fc.GetTestStatSampler();
202 ProfileLikelihoodTestStat *testStat = dynamic_cast<ProfileLikelihoodTestStat *>(toymcsampler->GetTestStatistic());
203
204 // Since this tool needs to throw toy MC the PDF needs to be
205 // extended or the tool needs to know how many entries in a dataset
206 // per pseudo experiment.
207 // In the 'number counting form' where the entries in the dataset
208 // are counts, and not values of discriminating variables, the
209 // datasets typically only have one entry and the PDF is not
210 // extended.
211 if (!mc->GetPdf()->canBeExtended()) {
212 if (data->numEntries() == 1)
213 fc.FluctuateNumDataEntries(false);
214 else
215 cout << "Not sure what to do about this model" << endl;
216 }
217
218 if (mc->GetGlobalObservables()) {
219 cout << "will use global observables for unconditional ensemble" << endl;
220 mc->GetGlobalObservables()->Print();
221 toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());
222 }
223
224 // Now get the interval
225 PointSetInterval *interval = fc.GetInterval();
226 ConfidenceBelt *belt = fc.GetConfidenceBelt();
227
228 // print out the interval on the first Parameter of Interest
229 cout << "\n95% interval on " << firstPOI->GetName() << " is : [" << interval->LowerLimit(*firstPOI) << ", "
230 << interval->UpperLimit(*firstPOI) << "] " << endl;
231
232 // get observed UL and value of test statistic evaluated there
234 double observedUL = interval->UpperLimit(*firstPOI);
235 firstPOI->setVal(observedUL);
236 double obsTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*data, tmpPOI);
237
238 // Ask the calculator which points were scanned
239 RooDataSet *parameterScan = (RooDataSet *)fc.GetPointsToScan();
241
242 // make a histogram of parameter vs. threshold
244 new TH1F("histOfThresholds", "", parameterScan->numEntries(), firstPOI->getMin(), firstPOI->getMax());
245 histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
246 histOfThresholds->GetYaxis()->SetTitle("Threshold");
247
248 // loop through the points that were tested and ask confidence belt
249 // what the upper/lower thresholds were.
250 // For FeldmanCousins, the lower cut off is always 0
251 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
252 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
253 // cout <<"get threshold"<<endl;
254 double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
255 double poiVal = tmpPoint->getRealValue(firstPOI->GetName());
257 }
258 TCanvas *c1 = new TCanvas();
259 c1->Divide(2);
260 c1->cd(1);
261 histOfThresholds->SetMinimum(0);
262 histOfThresholds->Draw();
263 c1->cd(2);
264
265 // -------------------------------------------------------
266 // Now we generate the expected bands and power-constraint
267
268 // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
269 std::unique_ptr<RooAbsReal> nll{mc->GetPdf()->createNLL(*data)};
270 std::unique_ptr<RooAbsReal> profile{nll->createProfile(*mc->GetParametersOfInterest())};
271 firstPOI->setVal(0.);
272 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
274 if (mc->GetNuisanceParameters())
275 poiAndNuisance->add(*mc->GetNuisanceParameters());
276 poiAndNuisance->add(*mc->GetParametersOfInterest());
277 w->saveSnapshot("paramsToGenerateData", *poiAndNuisance);
279 cout << "\nWill use these parameter points to generate pseudo data for bkg only" << endl;
280 paramsToGenerateData->Print("v");
281
283 unconditionalObs.add(*mc->GetObservables());
284 unconditionalObs.add(*mc->GetGlobalObservables()); // comment this out for the original conditional ensemble
285
286 double CLb = 0;
287 double CLbinclusive = 0;
288
289 // Now we generate background only and find distribution of upper limits
290 TH1F *histOfUL = new TH1F("histOfUL", "", 100, 0, firstPOI->getMax());
291 histOfUL->GetXaxis()->SetTitle("Upper Limit (background only)");
292 histOfUL->GetYaxis()->SetTitle("Entries");
293 for (int imc = 0; imc < nToyMC; ++imc) {
294
295 // set parameters back to values for generating pseudo data
296 // cout << "\n get current nuis, set vals, print again" << endl;
297 w->loadSnapshot("paramsToGenerateData");
298 // poiAndNuisance->Print("v");
299
300 std::unique_ptr<RooDataSet> toyData;
301 // now generate a toy dataset for the main measurement
302 if (!mc->GetPdf()->canBeExtended()) {
303 if (data->numEntries() == 1)
304 toyData = std::unique_ptr<RooDataSet>{mc->GetPdf()->generate(*mc->GetObservables(), 1)};
305 else
306 cout << "Not sure what to do about this model" << endl;
307 } else {
308 // cout << "generating extended dataset"<<endl;
309 toyData = std::unique_ptr<RooDataSet>{mc->GetPdf()->generate(*mc->GetObservables(), Extended())};
310 }
311
312 // generate global observables
313 std::unique_ptr<RooDataSet> one{mc->GetPdf()->generateSimGlobal(*mc->GetGlobalObservables(), 1)};
314 const RooArgSet *values = one->get();
315 std::unique_ptr<RooArgSet> allVars{mc->GetPdf()->getVariables()};
316 allVars->assign(*values);
317
318 // get test stat at observed UL in observed data
319 firstPOI->setVal(observedUL);
320 double toyTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData, tmpPOI);
321 // toyData->get()->Print("v");
322 // cout <<"obsTSatObsUL " <<obsTSatObsUL << "toyTS " << toyTSatObsUL << endl;
323 if (obsTSatObsUL < toyTSatObsUL) // not sure about <= part yet
324 CLb += (1.) / nToyMC;
325 if (obsTSatObsUL <= toyTSatObsUL) // not sure about <= part yet
326 CLbinclusive += (1.) / nToyMC;
327
328 // loop over points in belt to find upper limit for this toy data
329 double thisUL = 0;
330 for (Int_t i = 0; i < parameterScan->numEntries(); ++i) {
331 tmpPoint = (RooArgSet *)parameterScan->get(i)->clone("temp");
332 double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
333 firstPOI->setVal(tmpPoint->getRealValue(firstPOI->GetName()));
334 // double thisTS = profile->getVal();
335 double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData, tmpPOI);
336
337 // cout << "poi = " << firstPOI->getVal()
338 // << " max is " << arMax << " this profile = " << thisTS << endl;
339 // cout << "thisTS = " << thisTS<<endl;
340 if (thisTS <= arMax) {
341 thisUL = firstPOI->getVal();
342 } else {
343 break;
344 }
345 }
346
347 histOfUL->Fill(thisUL);
348
349 // for few events, data is often the same, and UL is often the same
350 // cout << "thisUL = " << thisUL<<endl;
351 }
352 histOfUL->Draw();
353 c1->SaveAs("two-sided_upper_limit_output.pdf");
354
355 // if you want to see a plot of the sampling distribution for a particular scan point:
356 /*
357 SamplingDistPlot sampPlot;
358 int indexInScan = 0;
359 tmpPoint = (RooArgSet*) parameterScan->get(indexInScan)->clone("temp");
360 firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
361 toymcsampler->SetParametersForTestStat(tmpPOI);
362 SamplingDistribution* samp = toymcsampler->GetSamplingDistribution(*tmpPoint);
363 sampPlot.AddSamplingDistribution(samp);
364 sampPlot.Draw();
365 */
366
367 // Now find bands and power constraint
368 Double_t *bins = histOfUL->GetIntegral();
369 TH1F *cumulative = (TH1F *)histOfUL->Clone("cumulative");
370 cumulative->SetContent(bins);
371 double band2sigDown = 0, band1sigDown = 0, bandMedian = 0, band1sigUp = 0, band2sigUp = 0;
372 for (int i = 1; i <= cumulative->GetNbinsX(); ++i) {
373 if (bins[i] < RooStats::SignificanceToPValue(2))
374 band2sigDown = cumulative->GetBinCenter(i);
375 if (bins[i] < RooStats::SignificanceToPValue(1))
376 band1sigDown = cumulative->GetBinCenter(i);
377 if (bins[i] < 0.5)
378 bandMedian = cumulative->GetBinCenter(i);
379 if (bins[i] < RooStats::SignificanceToPValue(-1))
380 band1sigUp = cumulative->GetBinCenter(i);
381 if (bins[i] < RooStats::SignificanceToPValue(-2))
382 band2sigUp = cumulative->GetBinCenter(i);
383 }
384 cout << "-2 sigma band " << band2sigDown << endl;
385 cout << "-1 sigma band " << band1sigDown << " [Power Constraint)]" << endl;
386 cout << "median of band " << bandMedian << endl;
387 cout << "+1 sigma band " << band1sigUp << endl;
388 cout << "+2 sigma band " << band2sigUp << endl;
389
390 // print out the interval on the first Parameter of Interest
391 cout << "\nobserved 95% upper-limit " << interval->UpperLimit(*firstPOI) << endl;
392 cout << "CLb strict [P(toy>obs|0)] for observed 95% upper-limit " << CLb << endl;
393 cout << "CLb inclusive [P(toy>=obs|0)] for observed 95% upper-limit " << CLbinclusive << endl;
394}
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
#define gROOT
Definition TROOT.h:411
R__EXTERN TSystem * gSystem
Definition TSystem.h:572
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
Container class to hold unbinned data.
Definition RooDataSet.h:32
Variable that can be changed from the outside.
Definition RooRealVar.h:37
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
< A class that holds configuration information for a model using a workspace as a store
Definition ModelConfig.h:34
PointSetInterval is a concrete implementation of the ConfInterval interface.
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23
A ROOT file is an on-disk file, usually with extension .root, that stores objects in a file-system-li...
Definition TFile.h:131
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:3765
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:879
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition TSystem.cxx:1309
return c1
Definition legend1.C:41
double nll(double pdf, double weight, int binnedL, int doBinOffset)
Definition MathFuncs.h:388
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:67
Namespace for the RooStats classes.
Definition CodegenImpl.h:61
double SignificanceToPValue(double Z)
returns p-value corresponding to a 1-sided significance