Logo ROOT  
Reference Guide
FeldmanCousins.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer January 2009
3
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12
14
16
18
20
22
25
26#include "RooDataSet.h"
27#include "RooDataHist.h"
28#include "RooGlobalFunc.h"
29#include "RooMsgService.h"
30
31/** \class RooStats::FeldmanCousins
32 \ingroup Roostats
33
34The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a
35specific configuration of the more general NeymanConstruction. It is a concrete
36implementation of the IntervalCalculator interface that, which uses the
37NeymanConstruction in a particular way. As the name suggests, it returns a
38ConfidenceInterval. In particular, it produces a RooStats::PointSetInterval,
39which is a concrete implementation of the ConfInterval interface.
40
41The Neyman Construction is not a uniquely defined statistical technique, it
42requires that one specify an ordering rule or ordering principle, which is
43usually encoded by choosing a specific test statistic and limits of integration
44(corresponding to upper/lower/central limits). As a result, this class must be
45configured with the corresponding information before it can produce an interval.
46
47In the case of the Feldman-Cousins approach, the ordering principle is the
48likelihood ratio -- motivated by the Neyman-Pearson lemma. When nuisance
49parameters are involved, the profile likelihood ratio is the natural
50generalization. One may either choose to perform the construction over the full
51space of the nuisance parameters, or restrict the nuisance parameters to their
52conditional MLE (eg. profiled values).
53
54*/
55
57
58using namespace RooFit;
59using namespace RooStats;
60using namespace std;
61
62////////////////////////////////////////////////////////////////////////////////
63/// standard constructor
64
65FeldmanCousins::FeldmanCousins(RooAbsData& data, ModelConfig& model) :
66 fSize(0.05),
67 fModel(model),
68 fData(data),
69 fTestStatSampler(0),
70 fPointsToTest(0),
71 fPOIToTest(0),
72 fConfBelt(nullptr),
73 fAdaptiveSampling(false),
74 fAdditionalNToysFactor(1.),
75 fNbins(10),
76 fFluctuateData(true),
77 fDoProfileConstruction(true),
78 fSaveBeltToFile(false),
79 fCreateBelt(false)
80{
81}
82
83////////////////////////////////////////////////////////////////////////////////
84/// destructor
85
88 if(fPOIToTest) delete fPOIToTest;
90}
91
92////////////////////////////////////////////////////////////////////////////////
93/// set the model
94
96 fModel = model;
97}
98
99////////////////////////////////////////////////////////////////////////////////
100
103 this->CreateTestStatSampler();
104 return fTestStatSampler;
105}
106
107////////////////////////////////////////////////////////////////////////////////
108/// specify the Test Statistic and create a ToyMC test statistic sampler
109
111 // use the profile likelihood ratio as the test statistic
113
114 // create the ToyMC test statistic sampler
115 fTestStatSampler = new ToyMCSampler(*testStatistic,int(fAdditionalNToysFactor*50./fSize)) ;
120
122 ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point = " << (int) (fAdditionalNToysFactor*50./fSize) << endl;
123 } else{
124 ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point: adaptive" << endl;
125 }
126 if(fFluctuateData){
127 ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will fluctuate about expectation" << endl;
128 } else{
129 ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will not fluctuate, will always be " << fData.numEntries() << endl;
131 }
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// specify the parameter points to perform the construction.
136/// allow ability to profile on some nuisance parameters
137
139 // get ingredients
140 RooAbsPdf* pdf = fModel.GetPdf();
141 if (!pdf ){
142 ooccoutE(&fModel,Generation) << "FeldmanCousins: ModelConfig has no PDF" << endl;
143 return;
144 }
145
146 // get list of all parameters
149 parameters->add(*fModel.GetNuisanceParameters());
150
151
153 // if parameters include nuisance parameters, do profile construction
154 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has nuisance parameters, will do profile construction" << endl;
155
156 // set nbins for the POI
157 for (auto *myarg2 : static_range_cast<RooRealVar *>(*fModel.GetParametersOfInterest())){
158 myarg2->setBins(fNbins);
159 }
160
161 // get dataset for POI scan
162 // RooDataHist* parameterScan = nullptr;
163 RooAbsData* parameterScan = nullptr;
164 if(fPOIToTest)
165 parameterScan = fPOIToTest;
166 else
167 parameterScan = new RooDataHist("parameterScan", "", *fModel.GetParametersOfInterest());
168
169
170 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
171 // make profile construction
174 RooAbsReal* nll = pdf->createNLL(fData,RooFit::CloneData(false));
176
177 RooDataSet* profileConstructionPoints = new RooDataSet("profileConstruction",
178 "profileConstruction",
179 *parameters);
180
181
182 for(Int_t i=0; i<parameterScan->numEntries(); ++i){
183 // here's where we figure out the profiled value of nuisance parameters
184 parameters->assign(*parameterScan->get(i));
185 profile->getVal();
186 profileConstructionPoints->add(*parameters);
187 }
189 delete profile;
190 delete nll;
191 if(!fPOIToTest) delete parameterScan;
192
193 // done
194 fPointsToTest = profileConstructionPoints;
195
196 } else{
197 // Do full construction
198 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;
199
200 for (auto *myarg : static_range_cast<RooRealVar *>(*parameters)){
201 myarg->setBins(fNbins);
202 }
203
204 RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
205 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
206
207 fPointsToTest = parameterScan;
208 }
209
210 delete parameters;
211
212}
213
214////////////////////////////////////////////////////////////////////////////////
215/// Main interface to get a RooStats::ConfInterval.
216/// It constructs a RooStats::PointSetInterval.
217
219 // local variables
220 // RooAbsData* data = fData; //fWS->data(fDataName);
221
222 // fill in implied variables given data
224
225 // create the test statistic sampler (private data member fTestStatSampler)
227 this->CreateTestStatSampler();
228
230
231 if(!fFluctuateData)
233
234 // create parameter points to perform construction (private data member fPointsToTest)
235 this->CreateParameterPoints();
236
237 // Create a Neyman Construction
239 // configure it
241 nc.SetTestSize(fSize); // set size of test
242 // nc.SetParameters( fModel.GetParametersOfInterest);
244 nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
245 nc.SetData(fData);
250 if (fCreateBelt)
252
253 // use it
254 return nc.GetInterval();
255}
size_t fSize
#define ooccoutP(o, a)
Definition: RooMsgService.h:58
#define ooccoutE(o, a)
Definition: RooMsgService.h:60
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
virtual const RooArgSet * get() const
Definition: RooAbsData.h:106
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:994
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:480
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
void add(const RooArgSet &row, double weight=1.0, double weightError=0.0) override
Add one ore more rows of data.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
bool fSaveBeltToFile
controls use if ConfidenceBelt should be saved to a TFile
ConfidenceBelt * fConfBelt
RooAbsData & fData
data set
void CreateParameterPoints() const
initializes fPointsToTest data member (mutable)
void SetModel(const ModelConfig &) override
set the model
bool fFluctuateData
tell ToyMCSampler to fluctuate number of entries in dataset
double fSize
size of the test (eg. specified rate of Type I error)
Int_t fNbins
number of samples per variable
PointSetInterval * GetInterval() const override
Main interface to get a ConfInterval (will be a PointSetInterval)
TestStatSampler * GetTestStatSampler() const
Returns instance of TestStatSampler.
void CreateTestStatSampler() const
initializes fTestStatSampler data member (mutable)
double fAdditionalNToysFactor
give user ability to ask for more toys
RooAbsData * fPOIToTest
value of POI points to perform the construction
bool fCreateBelt
controls use if ConfidenceBelt should be saved to a TFile
bool fDoProfileConstruction
instead of full construction over nuisance parameters, do profile
RooAbsData * fPointsToTest
points to perform the construction
ToyMCSampler * fTestStatSampler
the test statistic sampler
~FeldmanCousins() override
destructor
bool fAdaptiveSampling
controls use of adaptive sampling algorithm
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
void GuessObsAndNuisance(const RooAbsData &data, bool printModelConfig=true)
Makes sensible guesses of observables, parameters of interest and nuisance parameters if one or multi...
Definition: ModelConfig.cxx:66
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
Definition: ModelConfig.h:234
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
Definition: ModelConfig.h:237
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
Definition: ModelConfig.h:246
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Definition: ModelConfig.h:231
NeymanConstruction is a concrete implementation of the NeymanConstruction interface that,...
PointSetInterval * GetInterval() const override
Main interface to get a ConfInterval (will be a PointSetInterval)
void SetTestSize(double size) override
set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
void SaveBeltToFile(bool flag=true)
save the confidence belt to a file
ConfidenceBelt * GetConfidenceBelt()
Get confidence belt. This requires that CreateConfBelt() has been called.
void SetTestStatSampler(TestStatSampler &sampler)
in addition to interface we also need: Set the TestStatSampler (eg.
void SetLeftSideTailFraction(double leftSideFraction=0.)
fLeftSideTailFraction*fSize defines lower edge of acceptance region.
void UseAdaptiveSampling(bool flag=true)
adaptive sampling algorithm to speed up interval calculation
void AdditionalNToysFactor(double fact)
give user ability to ask for more toys
void CreateConfBelt(bool flag=true)
should create confidence belt
void SetParameterPointsToTest(RooAbsData &pointsToTest)
User-defined set of points to test.
void SetData(RooAbsData &data) override
Set the DataSet.
PointSetInterval is a concrete implementation of the ConfInterval interface.
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:67
void SetParametersForTestStat(const RooArgSet &nullpoi) override
Set the Pdf, add to the workspace if not already there.
Definition: ToyMCSampler.h:150
void SetObservables(const RooArgSet &o) override
specify the observables in the dataset (needed to evaluate the test statistic)
Definition: ToyMCSampler.h:167
void SetPdf(RooAbsPdf &pdf) override
Set the Pdf, add to the workspace if not already there.
Definition: ToyMCSampler.h:154
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:144
RooCmdArg CloneData(bool flag)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:59
@ Generation
Definition: RooGlobalFunc.h:61
Namespace for the RooStats classes.
Definition: Asimov.h:19