Logo ROOT   master
Reference Guide
MCMCCalculator.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Authors: Kevin Belasco 17/06/2009
3 // Authors: Kyle Cranmer 17/06/2009
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 
13 /** \class RooStats::MCMCCalculator
14  \ingroup Roostats
15 
16  Bayesian Calculator estimating an interval or a credible region using the
17  Markov-Chain Monte Carlo method to integrate the likelihood function with the
18  prior to obtain the posterior function.
19 
20  By using the Markov-Chain Monte Carlo methods this calculator can work with
21  model which require the integration of a large number of parameters.
22 
23  MCMCCalculator is a concrete implementation of IntervalCalculator. It uses a
24  MetropolisHastings object to construct a Markov Chain of data points in the
25  parameter space. From this Markov Chain, this class can generate a
26  MCMCInterval as per user specification.
27 
28  The interface allows one to pass the model, data, and parameters via a
29  workspace and then specify them with names.
30 
31  After configuring the calculator, one only needs to ask GetInterval(), which
32  will return an ConfInterval (MCMCInterval in this case).
33  */
34 
35 #include "Rtypes.h"
36 #include "RooGlobalFunc.h"
37 #include "RooAbsReal.h"
38 #include "RooArgSet.h"
39 #include "RooArgList.h"
40 #include "RooStats/ModelConfig.h"
41 #include "RooStats/RooStatsUtils.h"
44 #include "RooStats/MarkovChain.h"
45 #include "RooStats/MCMCInterval.h"
46 #include "TIterator.h"
48 #include "RooStats/PdfProposal.h"
49 #include "RooProdPdf.h"
50 
52 
53 using namespace RooFit;
54 using namespace RooStats;
55 using namespace std;
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// default constructor
59 
60 MCMCCalculator::MCMCCalculator() :
61  fPropFunc(0),
62  fPdf(0),
63  fPriorPdf(0),
64  fData(0),
65  fAxes(0)
66 {
67  fNumIters = 0;
68  fNumBurnInSteps = 0;
69  fNumBins = 0;
70  fUseKeys = kFALSE;
72  fSize = -1;
74  fLeftSideTF = -1;
75  fEpsilon = -1;
76  fDelta = -1;
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 /// constructor from a Model Config with a basic settings package configured
81 /// by SetupBasicUsage()
82 
84  fPropFunc(0),
85  fData(&data),
86  fAxes(0)
87 {
88  SetModel(model);
90 }
91 
93  // set the model
94  fPdf = model.GetPdf();
95  fPriorPdf = model.GetPriorPdf();
96  fPOI.removeAll();
100  if (model.GetParametersOfInterest())
101  fPOI.add(*model.GetParametersOfInterest());
102  if (model.GetNuisanceParameters())
104  if (model.GetConditionalObservables())
106  if (model.GetGlobalObservables())
107  fGlobalObs.add( *(model.GetGlobalObservables() ) );
108 
109 }
110 
111 ////////////////////////////////////////////////////////////////////////////////
112 /// Constructor for automatic configuration with basic settings. Uses a
113 /// UniformProposal, 10,000 iterations, 40 burn in steps, 50 bins for each
114 /// RooRealVar, determines interval by histogram. Finds a 95% confidence
115 /// interval.
116 
118 {
119  fPropFunc = 0;
120  fNumIters = 10000;
121  fNumBurnInSteps = 40;
122  fNumBins = 50;
123  fUseKeys = kFALSE;
125  SetTestSize(0.05);
127  fLeftSideTF = -1;
128  fEpsilon = -1;
129  fDelta = -1;
130 }
131 
132 ////////////////////////////////////////////////////////////////////////////////
133 
135 {
136  if (a < 0 || a > 1) {
137  coutE(InputArguments) << "MCMCCalculator::SetLeftSideTailFraction: "
138  << "Fraction must be in the range [0, 1]. "
139  << a << "is not allowed." << endl;
140  return;
141  }
142 
143  fLeftSideTF = a;
145 }
146 
147 ////////////////////////////////////////////////////////////////////////////////
148 /// Main interface to get a RooStats::ConfInterval.
149 
151 {
152 
153  if (!fData || !fPdf ) return 0;
154  if (fPOI.getSize() == 0) return 0;
155 
156  if (fSize < 0) {
157  coutE(InputArguments) << "MCMCCalculator::GetInterval: "
158  << "Test size/Confidence level not set. Returning NULL." << endl;
159  return NULL;
160  }
161 
162  // if a proposal function has not been specified create a default one
163  bool useDefaultPropFunc = (fPropFunc == 0);
164  bool usePriorPdf = (fPriorPdf != 0);
165  if (useDefaultPropFunc) fPropFunc = new UniformProposal();
166 
167  // if prior is given create product
168  RooAbsPdf * prodPdf = fPdf;
169  if (usePriorPdf) {
170  TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );
171  prodPdf = new RooProdPdf(prodName,prodName,RooArgList(*fPdf,*fPriorPdf) );
172  }
173 
174  RooArgSet* constrainedParams = prodPdf->getParameters(*fData);
176  delete constrainedParams;
177 
178  RooArgSet* params = nll->getParameters(*fData);
179  RemoveConstantParameters(params);
180  if (fNumBins > 0) {
181  SetBins(*params, fNumBins);
183  if (dynamic_cast<PdfProposal*>(fPropFunc)) {
184  RooArgSet* proposalVars = ((PdfProposal*)fPropFunc)->GetPdf()->
185  getParameters((RooAbsData*)NULL);
186  SetBins(*proposalVars, fNumBins);
187  }
188  }
189 
191  mh.SetFunction(*nll);
194  mh.SetParameters(*params);
198 
199  MarkovChain* chain = mh.ConstructChain();
200 
201  TString name = TString("MCMCInterval_") + TString(GetName() );
202  MCMCInterval* interval = new MCMCInterval(name, fPOI, *chain);
203  if (fAxes != NULL)
204  interval->SetAxes(*fAxes);
205  if (fNumBurnInSteps > 0)
207  interval->SetUseKeys(fUseKeys);
208  interval->SetUseSparseHist(fUseSparseHist);
209  interval->SetIntervalType(fIntervalType);
212  if (fEpsilon >= 0)
213  interval->SetEpsilon(fEpsilon);
214  if (fDelta >= 0)
215  interval->SetDelta(fDelta);
216  interval->SetConfidenceLevel(1.0 - fSize);
217 
218  if (useDefaultPropFunc) delete fPropFunc;
219  if (usePriorPdf) delete prodPdf;
220  delete nll;
221  delete params;
222 
223  return interval;
224 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:918
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void SetSign(enum FunctionSign sign)
MCMCCalculator()
default constructor
#define coutE(a)
Definition: RooMsgService.h:33
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
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 SetNumBurnInSteps(Int_t numBurnInSteps)
set the number of steps in the chain to discard as burn-in, starting from the first ...
Definition: MCMCInterval.h:158
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooProdPdf is an efficient implementation of a product of PDFs of the form .
Definition: RooProdPdf.h:31
UniformProposal is a concrete implementation of the ProposalFunction interface for use with a Markov ...
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void SetType(enum FunctionType type)
virtual MarkovChain * ConstructChain()
void SetupBasicUsage()
Constructor for automatic configuration with basic settings.
Basic string class.
Definition: TString.h:131
STL namespace.
RooCmdArg GlobalObservables(const RooArgSet &globs)
This class uses the Metropolis-Hastings algorithm to construct a Markov Chain of data points using Mo...
virtual void SetParameters(const RooArgSet &set)
virtual void SetNumIters(Int_t numIters)
virtual void SetEpsilon(Double_t epsilon)
set the acceptable level or error for Keys interval determination
Definition: MCMCInterval.h:226
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void SetFunction(RooAbsReal &function)
RooCmdArg Constrain(const RooArgSet &params)
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
Definition: ModelConfig.h:246
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return NULL if not existing)
Definition: ModelConfig.h:252
Int_t getSize() const
virtual void SetProposalFunction(ProposalFunction &proposalFunction)
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) ...
auto * a
Definition: textangle.C:12
RooCmdArg ConditionalObservables(const RooArgSet &set)
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
Definition: MCMCInterval.h:239
virtual void SetChainParameters(const RooArgSet &set)
virtual void SetLeftSideTailFraction(Double_t a)
set the left-side tail fraction for a tail-fraction interval
Definition: MCMCInterval.h:247
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual void SetConfidenceLevel(Double_t cl)
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
const Bool_t kFALSE
Definition: RtypesCore.h:90
void SetBins(const RooAbsCollection &coll, Int_t numBins) const
virtual void SetUseKeys(Bool_t useKeys)
set whether to use kernel estimation to determine the interval
Definition: MCMCInterval.h:162
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
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
#define ClassImp(name)
Definition: Rtypes.h:361
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:545
ProposalFunction * fPropFunc
enum MCMCInterval::IntervalType fIntervalType
PdfProposal is a concrete implementation of the ProposalFunction interface.
Definition: PdfProposal.h:30
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual MCMCInterval * GetInterval() const
Main interface to get a ConfInterval.
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:255
virtual void SetDelta(Double_t delta)
kbelasco: The inner-workings of the class really should not be exposed like this in a comment...
Definition: MCMCInterval.h:260
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
virtual void SetLeftSideTailFraction(Double_t a)
Set the left side tail fraction.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
char name[80]
Definition: TGX11.cxx:109
virtual void SetModel(const ModelConfig &model)
Set the Model.
virtual void SetUseSparseHist(Bool_t useSparseHist)
set whether to use a sparse histogram.
Definition: MCMCInterval.h:166