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"
48#include "RooProdPdf.h"
49
51
52using namespace RooFit;
53using namespace RooStats;
54using namespace std;
55
56////////////////////////////////////////////////////////////////////////////////
57/// default constructor
58
59MCMCCalculator::MCMCCalculator() :
60 fPropFunc(0),
61 fPdf(0),
62 fPriorPdf(0),
63 fData(0),
64 fAxes(0)
65{
66 fNumIters = 0;
68 fNumBins = 0;
69 fUseKeys = false;
70 fUseSparseHist = false;
71 fSize = -1;
73 fLeftSideTF = -1;
74 fEpsilon = -1;
75 fDelta = -1;
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// constructor from a Model Config with a basic settings package configured
80/// by SetupBasicUsage()
81
83 fPropFunc(0),
84 fData(&data),
85 fAxes(0)
86{
87 SetModel(model);
89}
90
92 // set the model
93 fPdf = model.GetPdf();
94 fPriorPdf = model.GetPriorPdf();
99 if (model.GetParametersOfInterest())
101 if (model.GetNuisanceParameters())
103 if (model.GetConditionalObservables())
105 if (model.GetGlobalObservables())
106 fGlobalObs.add( *(model.GetGlobalObservables() ) );
107
108}
109
110////////////////////////////////////////////////////////////////////////////////
111/// Constructor for automatic configuration with basic settings. Uses a
112/// UniformProposal, 10,000 iterations, 40 burn in steps, 50 bins for each
113/// RooRealVar, determines interval by histogram. Finds a 95% confidence
114/// interval.
115
117{
118 fPropFunc = 0;
119 fNumIters = 10000;
120 fNumBurnInSteps = 40;
121 fNumBins = 50;
122 fUseKeys = false;
123 fUseSparseHist = false;
124 SetTestSize(0.05);
126 fLeftSideTF = -1;
127 fEpsilon = -1;
128 fDelta = -1;
129}
130
131////////////////////////////////////////////////////////////////////////////////
132
134{
135 if (a < 0 || a > 1) {
136 coutE(InputArguments) << "MCMCCalculator::SetLeftSideTailFraction: "
137 << "Fraction must be in the range [0, 1]. "
138 << a << "is not allowed." << endl;
139 return;
140 }
141
142 fLeftSideTF = a;
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// Main interface to get a RooStats::ConfInterval.
148
150{
151
152 if (!fData || !fPdf ) return 0;
153 if (fPOI.empty()) return 0;
154
155 if (fSize < 0) {
156 coutE(InputArguments) << "MCMCCalculator::GetInterval: "
157 << "Test size/Confidence level not set. Returning nullptr." << endl;
158 return nullptr;
159 }
160
161 // if a proposal function has not been specified create a default one
162 bool useDefaultPropFunc = (fPropFunc == 0);
163 bool usePriorPdf = (fPriorPdf != 0);
164 if (useDefaultPropFunc) fPropFunc = new UniformProposal();
165
166 // if prior is given create product
167 RooAbsPdf * prodPdf = fPdf;
168 if (usePriorPdf) {
169 TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );
170 prodPdf = new RooProdPdf(prodName,prodName,RooArgList(*fPdf,*fPriorPdf) );
171 }
172
173 RooArgSet* constrainedParams = prodPdf->getParameters(*fData);
175 delete constrainedParams;
176
177 RooArgSet* params = nll->getParameters(*fData);
179 if (fNumBins > 0) {
180 SetBins(*params, fNumBins);
182 if (dynamic_cast<PdfProposal*>(fPropFunc)) {
183 RooArgSet* proposalVars = ((PdfProposal*)fPropFunc)->GetPdf()->
184 getParameters((RooAbsData*)nullptr);
185 SetBins(*proposalVars, fNumBins);
186 }
187 }
188
190 mh.SetFunction(*nll);
193 mh.SetParameters(*params);
197
198 MarkovChain* chain = mh.ConstructChain();
199
200 TString name = TString("MCMCInterval_") + TString(GetName() );
201 MCMCInterval* interval = new MCMCInterval(name, fPOI, *chain);
202 if (fAxes != nullptr)
203 interval->SetAxes(*fAxes);
204 if (fNumBurnInSteps > 0)
206 interval->SetUseKeys(fUseKeys);
211 if (fEpsilon >= 0)
212 interval->SetEpsilon(fEpsilon);
213 if (fDelta >= 0)
214 interval->SetDelta(fDelta);
215 interval->SetConfidenceLevel(1.0 - fSize);
216
217 if (useDefaultPropFunc) delete fPropFunc;
218 if (usePriorPdf) delete prodPdf;
219 delete nll;
220 delete params;
221
222 return interval;
223}
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition: TGX11.cxx:110
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) 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:541
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
bool empty() const
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:996
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:33
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
RooArgSet fChainParams
parameters to store in the chain (if not specified they are all of them )
enum MCMCInterval::IntervalType fIntervalType
RooAbsPdf * fPdf
pointer to common PDF (owned by the workspace)
RooAbsData * fData
pointer to the data (owned by the workspace)
bool fUseKeys
whether to use kernel estimation to determine interval
bool fUseSparseHist
whether to use sparse histogram (if using hist at all)
void SetBins(const RooAbsCollection &coll, Int_t numBins) const
void SetupBasicUsage()
Constructor for automatic configuration with basic settings.
double fDelta
acceptable error for Keys cutoffs being equal
MCMCCalculator()
default constructor
void SetModel(const ModelConfig &model) override
Set the Model.
virtual void SetLeftSideTailFraction(double a)
Set the left side tail fraction.
RooArgSet fNuisParams
nuisance parameters for interval (not really used)
Int_t fNumIters
number of iterations to run metropolis algorithm
Int_t fNumBins
set the number of bins to create for each axis when constructing the interval
RooArgSet fConditionalObs
conditional observables
double fEpsilon
acceptable error for Keys interval determination
void SetTestSize(double size) override
set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
double fLeftSideTF
left side tail-fraction for interval
double fSize
size of the test (eg. specified rate of Type I error)
Int_t fNumBurnInSteps
number of iterations to discard as burn-in, starting from the first
RooArgSet fPOI
parameters of interest for interval
RooArgList * fAxes
which variables to put on each axis
RooArgSet fGlobalObs
global observables
MCMCInterval * GetInterval() const override
Main interface to get a ConfInterval.
RooAbsPdf * fPriorPdf
pointer to prior PDF (owned by the workspace)
ProposalFunction * fPropFunc
Proposal function for MCMC integration.
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:33
virtual void SetUseSparseHist(bool useSparseHist)
set whether to use a sparse histogram.
Definition: MCMCInterval.h:169
virtual void SetDelta(double delta)
kbelasco: The inner-workings of the class really should not be exposed like this in a comment,...
Definition: MCMCInterval.h:263
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
void SetConfidenceLevel(double cl) override
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
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:161
virtual void SetUseKeys(bool useKeys)
set whether to use kernel estimation to determine the interval
Definition: MCMCInterval.h:165
virtual void SetEpsilon(double epsilon)
set the acceptable level or error for Keys interval determination
Definition: MCMCInterval.h:229
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
Definition: MCMCInterval.h:242
virtual void SetLeftSideTailFraction(double a)
set the left-side tail fraction for a tail-fraction interval
Definition: MCMCInterval.h:250
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)
set the proposal function for suggesting new points for the MCMC
virtual void SetChainParameters(const RooArgSet &set)
specify the parameters to store in the chain if not specified all of them will be stored
virtual void SetParameters(const RooArgSet &set)
specify all the parameters of interest in the interval
virtual MarkovChain * ConstructChain()
main purpose of MetropolisHastings - run Metropolis-Hastings algorithm to generate Markov Chain of po...
virtual void SetNumIters(Int_t numIters)
set the number of iterations to run the metropolis algorithm
virtual void SetType(enum FunctionType type)
set the type of the function
virtual void SetFunction(RooAbsReal &function)
set the (likelihood) function
virtual void SetSign(enum FunctionSign sign)
set the sign of the function
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 nullptr if not existing)
Definition: ModelConfig.h:249
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
Definition: ModelConfig.h:252
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
Definition: ModelConfig.h:234
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
Definition: ModelConfig.h:237
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Definition: ModelConfig.h:231
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return nullptr if not existing)
Definition: ModelConfig.h:243
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 ...
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:136
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:67
TArc a
Definition: textangle.C:12