Logo ROOT   6.10/09
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 #include "RooStats/RooStatsUtils.h"
26 
27 #include "RooDataSet.h"
28 #include "RooDataHist.h"
29 #include "RooGlobalFunc.h"
30 #include "RooMsgService.h"
31 #include "TFile.h"
32 #include "TTree.h"
33 
34 /** \class RooStats::FeldmanCousins
35  \ingroup Roostats
36 
37 The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a
38 specific configuration of the more general NeymanConstruction. It is a concrete
39 implementation of the IntervalCalculator interface that, which uses the
40 NeymanConstruction in a particular way. As the name suggests, it returns a
41 ConfidenceInterval. In particular, it produces a RooStats::PointSetInterval,
42 which is a concrete implementation of the ConfInterval interface.
43 
44 The Neyman Construction is not a uniquely defined statistical technique, it
45 requires that one specify an ordering rule or ordering principle, which is
46 usually encoded by choosing a specific test statistic and limits of integration
47 (corresponding to upper/lower/central limits). As a result, this class must be
48 configured with the corresponding information before it can produce an interval.
49 
50 In the case of the Feldman-Cousins approach, the ordering principle is the
51 likelihood ratio -- motivated by the Neyman-Pearson lemma. When nuisance
52 parameters are involved, the profile likelihood ratio is the natural
53 generalization. One may either choose to perform the construction over the full
54 space of the nuisance parameters, or restrict the nuisance parameters to their
55 conditional MLE (eg. profiled values).
56 
57 */
58 
60 
61 using namespace RooFit;
62 using namespace RooStats;
63 using namespace std;
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 /// standard constructor
67 
68 FeldmanCousins::FeldmanCousins(RooAbsData& data, ModelConfig& model) :
69  fSize(0.05),
70  fModel(model),
71  fData(data),
72  fTestStatSampler(0),
73  fPointsToTest(0),
74  fPOIToTest(0),
75  fConfBelt(0),
76  fAdaptiveSampling(false),
77  fAdditionalNToysFactor(1.),
78  fNbins(10),
79  fFluctuateData(true),
80  fDoProfileConstruction(true),
81  fSaveBeltToFile(false),
82  fCreateBelt(false)
83 {
84 }
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// destructor
88 
90  if(fPointsToTest) delete fPointsToTest;
91  if(fPOIToTest) delete fPOIToTest;
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// set the model
97 
99  fModel = model;
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 
105  if(!fTestStatSampler)
106  this->CreateTestStatSampler();
107  return fTestStatSampler;
108 }
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 /// specify the Test Statistic and create a ToyMC test statistic sampler
112 
114  // use the profile likelihood ratio as the test statistic
116 
117  // create the ToyMC test statistic sampler
118  fTestStatSampler = new ToyMCSampler(*testStatistic,int(fAdditionalNToysFactor*50./fSize)) ;
119  fTestStatSampler->SetParametersForTestStat(*fModel.GetParametersOfInterest() );
120  if(fModel.GetObservables())
121  fTestStatSampler->SetObservables(*fModel.GetObservables());
122  fTestStatSampler->SetPdf(*fModel.GetPdf());
123 
124  if(!fAdaptiveSampling){
125  ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point = " << (int) (fAdditionalNToysFactor*50./fSize) << endl;
126  } else{
127  ooccoutP(&fModel,Generation) << "FeldmanCousins: ntoys per point: adaptive" << endl;
128  }
129  if(fFluctuateData){
130  ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will fluctuate about expectation" << endl;
131  } else{
132  ooccoutP(&fModel,Generation) << "FeldmanCousins: nEvents per toy will not fluctuate, will always be " << fData.numEntries() << endl;
133  fTestStatSampler->SetNEventsPerToy(fData.numEntries());
134  }
135 }
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// specify the parameter points to perform the construction.
139 /// allow ability to profile on some nuisance parameters
140 
142  // get ingredients
143  RooAbsPdf* pdf = fModel.GetPdf();
144  if (!pdf ){
145  ooccoutE(&fModel,Generation) << "FeldmanCousins: ModelConfig has no PDF" << endl;
146  return;
147  }
148 
149  // get list of all parameters
150  RooArgSet* parameters = new RooArgSet(*fModel.GetParametersOfInterest());
152  parameters->add(*fModel.GetNuisanceParameters());
153 
154 
156  // if parameters include nuisance parameters, do profile construction
157  ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has nuisance parameters, will do profile construction" << endl;
158 
159  // set nbins for the POI
161  RooRealVar *myarg2;
162  while ((myarg2 = dynamic_cast<RooRealVar*>(it2.Next()))) {
163  myarg2->setBins(fNbins);
164  }
165 
166  // get dataset for POI scan
167  // RooDataHist* parameterScan = NULL;
168  RooAbsData* parameterScan = NULL;
169  if(fPOIToTest)
170  parameterScan = fPOIToTest;
171  else
172  parameterScan = new RooDataHist("parameterScan", "", *fModel.GetParametersOfInterest());
173 
174 
175  ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
176  // make profile construction
179  RooAbsReal* nll = pdf->createNLL(fData,RooFit::CloneData(false));
181 
182  RooDataSet* profileConstructionPoints = new RooDataSet("profileConstruction",
183  "profileConstruction",
184  *parameters);
185 
186 
187  for(Int_t i=0; i<parameterScan->numEntries(); ++i){
188  // here's where we figure out the profiled value of nuisance parameters
189  *parameters = *parameterScan->get(i);
190  profile->getVal();
191  profileConstructionPoints->add(*parameters);
192  }
194  delete profile;
195  delete nll;
196  if(!fPOIToTest) delete parameterScan;
197 
198  // done
199  fPointsToTest = profileConstructionPoints;
200 
201  } else{
202  // Do full construction
203  ooccoutP(&fModel,Generation) << "FeldmanCousins: Model has no nuisance parameters" << endl;
204 
205  TIter it = parameters->createIterator();
206  RooRealVar *myarg;
207  while ((myarg = dynamic_cast<RooRealVar*>(it.Next()))) {
208  myarg->setBins(fNbins);
209  }
210 
211  RooDataHist* parameterScan = new RooDataHist("parameterScan", "", *parameters);
212  ooccoutP(&fModel,Generation) << "FeldmanCousins: # points to test = " << parameterScan->numEntries() << endl;
213 
214  fPointsToTest = parameterScan;
215  }
216 
217  delete parameters;
218 
219 }
220 
221 ////////////////////////////////////////////////////////////////////////////////
222 /// Main interface to get a RooStats::ConfInterval.
223 /// It constructs a RooStats::PointSetInterval.
224 
226  // local variables
227  // RooAbsData* data = fData; //fWS->data(fDataName);
228 
229  // fill in implied variables given data
231 
232  // create the test statistic sampler (private data member fTestStatSampler)
233  if(!fTestStatSampler)
234  this->CreateTestStatSampler();
235 
237 
238  if(!fFluctuateData)
240 
241  // create parameter points to perform construction (private data member fPointsToTest)
242  this->CreateParameterPoints();
243 
244  // Create a Neyman Construction
246  // configure it
248  nc.SetTestSize(fSize); // set size of test
249  // nc.SetParameters( fModel.GetParametersOfInterest);
251  nc.SetLeftSideTailFraction(0.); // part of definition of Feldman-Cousins
252  nc.SetData(fData);
258  // use it
259  return nc.GetInterval();
260 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:777
TIterator * createIterator(Bool_t dir=kIterForward) const
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
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:237
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:86
virtual const RooArgSet * get() const
Definition: RooAbsData.h:77
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) ...
RooCmdArg CloneData(Bool_t flag)
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
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.
void GuessObsAndNuisance(const RooAbsData &data)
guesses Observables and ParametersOfInterest if not already set
Definition: ModelConfig.cxx:50
RooFit::MsgLevel globalKillBelow() const
virtual void SetNEventsPerToy(const Int_t nevents)
Definition: ToyMCSampler.h:147
int Int_t
Definition: RtypesCore.h:41
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:52
STL namespace.
#define NULL
Definition: RtypesCore.h:88
void SetParameterPointsToTest(RooAbsData &pointsToTest)
User-defined set of points to test.
#define ooccoutE(o, a)
Definition: RooMsgService.h:54
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
RooDataSet 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
virtual ~FeldmanCousins()
destructor
void setBins(Int_t nBins, const char *name=0)
Definition: RooRealVar.h:77
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
void CreateTestStatSampler() const
initializes fTestStatSampler data member (mutable)
TObject * Next()
Definition: TCollection.h:153
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:37
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:29
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0)
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
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:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:222
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:225
PointSetInterval is a concrete implementation of the ConfInterval interface.
#define ClassImp(name)
Definition: Rtypes.h:336
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:53
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...
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
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:228
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:461
void SaveBeltToFile(bool flag=true)
save the confidence belt to a file
virtual void SetModel(const ModelConfig &)
set the model
ConfidenceBelt * fConfBelt
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269