Logo ROOT  
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"
46#include "TIterator.h"
49#include "RooProdPdf.h"
50
52
53using namespace RooFit;
54using namespace RooStats;
55using namespace std;
56
57////////////////////////////////////////////////////////////////////////////////
58/// default constructor
59
60MCMCCalculator::MCMCCalculator() :
61 fPropFunc(0),
62 fPdf(0),
63 fPriorPdf(0),
64 fData(0),
65 fAxes(0)
66{
67 fNumIters = 0;
69 fNumBins = 0;
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();
100 if (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;
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);
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);
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}
#define coutE(a)
Definition: RooMsgService.h:34
const Bool_t kFALSE
Definition: RtypesCore.h:88
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
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't match any of...
Definition: RooAbsArg.cxx:548
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
Int_t getSize() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:39
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:957
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
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
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
enum MCMCInterval::IntervalType fIntervalType
virtual MCMCInterval * GetInterval() const
Main interface to get a ConfInterval.
void SetBins(const RooAbsCollection &coll, Int_t numBins) const
void SetupBasicUsage()
Constructor for automatic configuration with basic settings.
MCMCCalculator()
default constructor
virtual void SetLeftSideTailFraction(Double_t a)
Set the left side tail fraction.
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)
virtual void SetModel(const ModelConfig &model)
Set the Model.
ProposalFunction * fPropFunc
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
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
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
virtual void SetLeftSideTailFraction(Double_t a)
set the left-side tail fraction for a tail-fraction interval
Definition: MCMCInterval.h:247
virtual void SetConfidenceLevel(Double_t cl)
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
Definition: MCMCInterval.h:239
virtual void SetEpsilon(Double_t epsilon)
set the acceptable level or error for Keys interval determination
Definition: MCMCInterval.h:226
virtual void SetUseKeys(Bool_t useKeys)
set whether to use kernel estimation to determine the interval
Definition: MCMCInterval.h:162
virtual void SetUseSparseHist(Bool_t useSparseHist)
set whether to use a sparse histogram.
Definition: MCMCInterval.h:166
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
This class uses the Metropolis-Hastings algorithm to construct a Markov Chain of data points using Mo...
virtual void SetProposalFunction(ProposalFunction &proposalFunction)
virtual void SetChainParameters(const RooArgSet &set)
virtual void SetParameters(const RooArgSet &set)
virtual MarkovChain * ConstructChain()
virtual void SetNumIters(Int_t numIters)
virtual void SetType(enum FunctionType type)
virtual void SetFunction(RooAbsReal &function)
virtual void SetSign(enum FunctionSign sign)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return NULL if not existing)
Definition: ModelConfig.h:252
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:255
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
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
Definition: ModelConfig.h:246
PdfProposal is a concrete implementation of the ProposalFunction interface.
Definition: PdfProposal.h:30
UniformProposal is a concrete implementation of the ProposalFunction interface for use with a Markov ...
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(const RooArgSet &globs)
@ InputArguments
Definition: RooGlobalFunc.h:68
RooCmdArg ConditionalObservables(const RooArgSet &set)
Namespace for the RooStats classes.
Definition: Asimov.h:20
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:68
auto * a
Definition: textangle.C:12