Logo ROOT   master
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 
15 #include "RooStats/RooStatsUtils.h"
16 
18 
19 #include "RooStats/ModelConfig.h"
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 
34 The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a
35 specific configuration of the more general NeymanConstruction. It is a concrete
36 implementation of the IntervalCalculator interface that, which uses the
37 NeymanConstruction in a particular way. As the name suggests, it returns a
38 ConfidenceInterval. In particular, it produces a RooStats::PointSetInterval,
39 which is a concrete implementation of the ConfInterval interface.
40 
41 The Neyman Construction is not a uniquely defined statistical technique, it
42 requires that one specify an ordering rule or ordering principle, which is
43 usually 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
45 configured with the corresponding information before it can produce an interval.
46 
47 In the case of the Feldman-Cousins approach, the ordering principle is the
48 likelihood ratio -- motivated by the Neyman-Pearson lemma. When nuisance
49 parameters are involved, the profile likelihood ratio is the natural
50 generalization. One may either choose to perform the construction over the full
51 space of the nuisance parameters, or restrict the nuisance parameters to their
52 conditional MLE (eg. profiled values).
53 
54 */
55 
57 
58 using namespace RooFit;
59 using namespace RooStats;
60 using namespace std;
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// standard constructor
64 
65 FeldmanCousins::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 
87  if(fPointsToTest) delete fPointsToTest;
88  if(fPOIToTest) delete fPOIToTest;
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// set the model
94 
96  fModel = model;
97 }
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 
102  if(!fTestStatSampler)
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)) ;
117  if(fModel.GetObservables())
120 
121  if(!fAdaptiveSampling){
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
147  RooArgSet* parameters = new RooArgSet(*fModel.GetParametersOfInterest());
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 = *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)
230  if(!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 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:916
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
ConfidenceBelt * GetConfidenceBelt()
Get confidence belt. This requires that CreateConfBelt() has been called.
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
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
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:155
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
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) ...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
RooFit::MsgLevel globalKillBelow() const
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:149
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
static RooMsgService & instance()
Return reference to singleton instance.
#define ooccoutP(o, a)
Definition: RooMsgService.h:54
STL namespace.
void SetParameterPointsToTest(RooAbsData &pointsToTest)
User-defined set of points to test.
#define ooccoutE(o, a)
Definition: RooMsgService.h:56
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
void AdditionalNToysFactor(double fact)
give user ability to ask for more toys
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual ~FeldmanCousins()
destructor
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name &#39;name&#39; for this variable.
Definition: RooRealVar.cxx:382
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.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...
Definition: ModelConfig.cxx:66
void CreateTestStatSampler() const
initializes fTestStatSampler data member (mutable)
TObject * Next()
Definition: TCollection.h:249
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
ToyMCSampler * fTestStatSampler
void setGlobalKillBelow(RooFit::MsgLevel level)
virtual void SetData(RooAbsData &data)
Set the DataSet.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
RooCmdArg CloneData(Bool_t flag)
void SetTestStatSampler(TestStatSampler &sampler)
in addition to interface we also need: Set the TestStatSampler (eg.
void UseAdaptiveSampling(bool flag=true)
adaptive sampling algorithm to speed up interval calculation
Namespace for the RooStats classes.
Definition: Asimov.h:19
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
PointSetInterval is a concrete implementation of the ConfInterval interface.
#define ClassImp(name)
Definition: Rtypes.h:361
virtual Int_t numEntries() const
Return the number of bins.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:173
void CreateParameterPoints() const
initializes fPointsToTest data member (mutable)
void CreateConfBelt(bool flag=true)
should create confidence belt
NeymanConstruction is a concrete implementation of the NeymanConstruction interface that...
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:160
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
TestStatSampler * GetTestStatSampler() const
Returns instance of TestStatSampler.
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
void SetLeftSideTailFraction(Double_t leftSideFraction=0.)
fLeftSideTailFraction*fSize defines lower edge of acceptance region.
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:515
void SaveBeltToFile(bool flag=true)
save the confidence belt to a file
virtual void SetModel(const ModelConfig &)
set the model
ConfidenceBelt * fConfBelt
size_t fSize
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
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 &#39;data&#39; argset, to the data set...