Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
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
158 RooRealVar *myarg2;
159 while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
160 myarg2->setBins(fNbins);
161 }
162
163 // get dataset for POI scan
164 // RooDataHist* parameterScan = NULL;
165 RooAbsData* parameterScan = NULL;
166 if(fPOIToTest)
167 parameterScan = fPOIToTest;
168 else
169 parameterScan = new RooDataHist("parameterScan", "", *fModel.GetParametersOfInterest());
170
171
172 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
173 // make profile construction
176 RooAbsReal* nll = pdf->createNLL(fData,RooFit::CloneData(false));
178
179 RooDataSet* profileConstructionPoints = new RooDataSet("profileConstruction",
180 "profileConstruction",
181 *parameters);
182
183
184 for(Int_t i=0; i<parameterScan->numEntries(); ++i){
185 // here's where we figure out the profiled value of nuisance parameters
186 parameters->assign(*parameterScan->get(i));
187 profile->getVal();
188 profileConstructionPoints->add(*parameters);
189 }
191 delete profile;
192 delete nll;
193 if(!fPOIToTest) delete parameterScan;
194
195 // done
196 fPointsToTest = profileConstructionPoints;
197
198 } else{
199 // Do full construction
200 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;
201
202 TIter it = parameters->createIterator();
203 RooRealVar *myarg;
204 while ((myarg = dynamic_cast<RooRealVar*>(it.Next()))) {
205 myarg->setBins(fNbins);
206 }
207
208 RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
209 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
210
211 fPointsToTest = parameterScan;
212 }
213
214 delete parameters;
215
216}
217
218////////////////////////////////////////////////////////////////////////////////
219/// Main interface to get a RooStats::ConfInterval.
220/// It constructs a RooStats::PointSetInterval.
221
223 // local variables
224 // RooAbsData* data = fData; //fWS->data(fDataName);
225
226 // fill in implied variables given data
228
229 // create the test statistic sampler (private data member fTestStatSampler)
231 this->CreateTestStatSampler();
232
234
235 if(!fFluctuateData)
237
238 // create parameter points to perform construction (private data member fPointsToTest)
239 this->CreateParameterPoints();
240
241 // Create a Neyman Construction
243 // configure it
245 nc.SetTestSize(fSize); // set size of test
246 // nc.SetParameters( fModel.GetParametersOfInterest);
248 nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
249 nc.SetData(fData);
254 if (fCreateBelt)
256
257 // use it
258 return nc.GetInterval();
259}
size_t fSize
#define ooccoutP(o, a)
#define ooccoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:364
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
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...
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:82
virtual const RooArgSet * get() const
Definition RooAbsData.h:128
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:64
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:94
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:35
The RooDataHist is a container class to hold N-dimensional binned data.
Definition RooDataHist.h:45
Int_t numEntries() const override
Return the number of bins.
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:36
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
ConfidenceBelt * fConfBelt
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
void CreateParameterPoints() const
initializes fPointsToTest data member (mutable)
FeldmanCousins(RooAbsData &data, ModelConfig &model)
Common constructor.
TestStatSampler * GetTestStatSampler() const
Returns instance of TestStatSampler.
void CreateTestStatSampler() const
initializes fTestStatSampler data member (mutable)
virtual ~FeldmanCousins()
destructor
ToyMCSampler * fTestStatSampler
virtual void SetModel(const ModelConfig &)
set the model
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...
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
NeymanConstruction is a concrete implementation of the NeymanConstruction interface that,...
virtual void SetTestSize(Double_t size)
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.
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
void SetTestStatSampler(TestStatSampler &sampler)
in addition to interface we also need: Set the TestStatSampler (eg.
void SetLeftSideTailFraction(Double_t 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.
virtual void SetData(RooAbsData &data)
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.
virtual void SetObservables(const RooArgSet &o)
virtual void SetPdf(RooAbsPdf &pdf)
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
TObject * Next()
RooCmdArg CloneData(Bool_t 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.
Namespace for the RooStats classes.
Definition Asimov.h:19