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#include "TFile.h"
31#include "TTree.h"
32
33/** \class RooStats::FeldmanCousins
34 \ingroup Roostats
35
36The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a
37specific configuration of the more general NeymanConstruction. It is a concrete
38implementation of the IntervalCalculator interface that, which uses the
39NeymanConstruction in a particular way. As the name suggests, it returns a
40ConfidenceInterval. In particular, it produces a RooStats::PointSetInterval,
41which is a concrete implementation of the ConfInterval interface.
42
43The Neyman Construction is not a uniquely defined statistical technique, it
44requires that one specify an ordering rule or ordering principle, which is
45usually encoded by choosing a specific test statistic and limits of integration
46(corresponding to upper/lower/central limits). As a result, this class must be
47configured with the corresponding information before it can produce an interval.
48
49In the case of the Feldman-Cousins approach, the ordering principle is the
50likelihood ratio -- motivated by the Neyman-Pearson lemma. When nuisance
51parameters are involved, the profile likelihood ratio is the natural
52generalization. One may either choose to perform the construction over the full
53space of the nuisance parameters, or restrict the nuisance parameters to their
54conditional MLE (eg. profiled values).
55
56*/
57
59
60using namespace RooFit;
61using namespace RooStats;
62using namespace std;
63
64////////////////////////////////////////////////////////////////////////////////
65/// standard constructor
66
67FeldmanCousins::FeldmanCousins(RooAbsData& data, ModelConfig& model) :
68 fSize(0.05),
69 fModel(model),
70 fData(data),
71 fTestStatSampler(0),
72 fPointsToTest(0),
73 fPOIToTest(0),
74 fConfBelt(nullptr),
75 fAdaptiveSampling(false),
76 fAdditionalNToysFactor(1.),
77 fNbins(10),
78 fFluctuateData(true),
79 fDoProfileConstruction(true),
80 fSaveBeltToFile(false),
81 fCreateBelt(false)
82{
83}
84
85////////////////////////////////////////////////////////////////////////////////
86/// destructor
87
90 if(fPOIToTest) delete fPOIToTest;
92}
93
94////////////////////////////////////////////////////////////////////////////////
95/// set the model
96
98 fModel = model;
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
105 this->CreateTestStatSampler();
106 return fTestStatSampler;
107}
108
109////////////////////////////////////////////////////////////////////////////////
110/// specify the Test Statistic and create a ToyMC test statistic sampler
111
113 // use the profile likelihood ratio as the test statistic
115
116 // create the ToyMC test statistic sampler
117 fTestStatSampler = new ToyMCSampler(*testStatistic,int(fAdditionalNToysFactor*50./fSize)) ;
122
124 ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point = " << (int) (fAdditionalNToysFactor*50./fSize) << endl;
125 } else{
126 ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point: adaptive" << endl;
127 }
128 if(fFluctuateData){
129 ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will fluctuate about expectation" << endl;
130 } else{
131 ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will not fluctuate, will always be " << fData.numEntries() << endl;
133 }
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// specify the parameter points to perform the construction.
138/// allow ability to profile on some nuisance parameters
139
141 // get ingredients
142 RooAbsPdf* pdf = fModel.GetPdf();
143 if (!pdf ){
144 ooccoutE(&fModel,Generation) << "FeldmanCousins: ModelConfig has no PDF" << endl;
145 return;
146 }
147
148 // get list of all parameters
151 parameters->add(*fModel.GetNuisanceParameters());
152
153
155 // if parameters include nuisance parameters, do profile construction
156 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has nuisance parameters, will do profile construction" << endl;
157
158 // set nbins for the POI
160 RooRealVar *myarg2;
161 while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
162 myarg2->setBins(fNbins);
163 }
164
165 // get dataset for POI scan
166 // RooDataHist* parameterScan = NULL;
167 RooAbsData* parameterScan = NULL;
168 if(fPOIToTest)
169 parameterScan = fPOIToTest;
170 else
171 parameterScan = new RooDataHist("parameterScan", "", *fModel.GetParametersOfInterest());
172
173
174 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
175 // make profile construction
178 RooAbsReal* nll = pdf->createNLL(fData,RooFit::CloneData(false));
180
181 RooDataSet* profileConstructionPoints = new RooDataSet("profileConstruction",
182 "profileConstruction",
183 *parameters);
184
185
186 for(Int_t i=0; i<parameterScan->numEntries(); ++i){
187 // here's where we figure out the profiled value of nuisance parameters
188 *parameters = *parameterScan->get(i);
189 profile->getVal();
190 profileConstructionPoints->add(*parameters);
191 }
193 delete profile;
194 delete nll;
195 if(!fPOIToTest) delete parameterScan;
196
197 // done
198 fPointsToTest = profileConstructionPoints;
199
200 } else{
201 // Do full construction
202 ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;
203
204 TIter it = parameters->createIterator();
205 RooRealVar *myarg;
206 while ((myarg = dynamic_cast<RooRealVar*>(it.Next()))) {
207 myarg->setBins(fNbins);
208 }
209
210 RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
211 ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
212
213 fPointsToTest = parameterScan;
214 }
215
216 delete parameters;
217
218}
219
220////////////////////////////////////////////////////////////////////////////////
221/// Main interface to get a RooStats::ConfInterval.
222/// It constructs a RooStats::PointSetInterval.
223
225 // local variables
226 // RooAbsData* data = fData; //fWS->data(fDataName);
227
228 // fill in implied variables given data
230
231 // create the test statistic sampler (private data member fTestStatSampler)
233 this->CreateTestStatSampler();
234
236
237 if(!fFluctuateData)
239
240 // create parameter points to perform construction (private data member fPointsToTest)
241 this->CreateParameterPoints();
242
243 // Create a Neyman Construction
245 // configure it
247 nc.SetTestSize(fSize); // set size of test
248 // nc.SetParameters( fModel.GetParametersOfInterest);
250 nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
251 nc.SetData(fData);
256 if (fCreateBelt)
258
259 // use it
260 return nc.GetInterval();
261}
size_t fSize
#define ooccoutP(o, a)
Definition: RooMsgService.h:54
#define ooccoutE(o, a)
Definition: RooMsgService.h:56
int Int_t
Definition: RtypesCore.h:43
#define ClassImp(name)
Definition: Rtypes.h:361
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:44
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:918
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:515
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual Int_t numEntries() const
Return the number of bins.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
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:35
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:418
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)
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...
Definition: ModelConfig.cxx:68
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
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.
Definition: ToyMCSampler.h:72
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:173
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:160
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:150
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:156
TObject * Next()
Definition: TCollection.h:249
RooCmdArg CloneData(Bool_t flag)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
@ Generation
Definition: RooGlobalFunc.h:67
Namespace for the RooStats classes.
Definition: Asimov.h:19