ROOT   Reference Guide
Searching...
No Matches
FeldmanCousins.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer January 2009
3
4/*************************************************************************
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 std::endl;
61
62////////////////////////////////////////////////////////////////////////////////
63/// standard constructor
64
66 fSize(0.05),
67 fModel(model),
68 fData(data),
69 fTestStatSampler(nullptr),
70 fPointsToTest(nullptr),
71 fPOIToTest(nullptr),
72 fConfBelt(nullptr),
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
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 std::unique_ptr<RooAbsReal> nll{pdf->createNLL(fData,RooFit::CloneData(false))};
175 std::unique_ptr<RooAbsReal> profile{nll->createProfile(*fModel.GetParametersOfInterest())};
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();
187 }
189 if(!fPOIToTest) delete parameterScan;
190
191 // done
192 fPointsToTest = profileConstructionPoints;
193
194 } else{
195 // Do full construction
196 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;
197
198 for (auto *myarg : static_range_cast<RooRealVar *>(*parameters)){
199 myarg->setBins(fNbins);
200 }
201
202 RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
203 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
204
205 fPointsToTest = parameterScan;
206 }
207
208 delete parameters;
209
210}
211
212////////////////////////////////////////////////////////////////////////////////
213/// Main interface to get a RooStats::ConfInterval.
214/// It constructs a RooStats::PointSetInterval.
215
217 // local variables
218 // RooAbsData* data = fData; //fWS->data(fDataName);
219
220 // fill in implied variables given data
222
223 // create the test statistic sampler (private data member fTestStatSampler)
225 this->CreateTestStatSampler();
226
228
229 if(!fFluctuateData)
231
232 // create parameter points to perform construction (private data member fPointsToTest)
233 this->CreateParameterPoints();
234
235 // Create a Neyman Construction
237 // configure it
239 nc.SetTestSize(fSize); // set size of test
240 // nc.SetParameters( fModel.GetParametersOfInterest);
242 nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
243 nc.SetData(fData);
248 if (fCreateBelt)
250
251 // use it
252 return nc.GetInterval();
253}
dim_t fSize
#define ooccoutP(o, a)
#define ooccoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
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...
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:163
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:39
Container class to hold unbinned data.
Definition RooDataSet.h:57
void add(const RooArgSet &row, double weight, double weightError)
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
FeldmanCousins(RooAbsData &data, ModelConfig &model)
Common constructor.
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)
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
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:35
void GuessObsAndNuisance(const RooAbsData &data, bool printModelConfig=true)
Makes sensible guesses of observables, parameters of interest and nuisance parameters if one or multi...
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
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.
adaptive sampling algorithm to speed up interval calculation
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.
void SetParametersForTestStat(const RooArgSet &nullpoi) override