Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
60 fPropFunc(nullptr),
61 fPdf(nullptr),
62 fPriorPdf(nullptr),
63 fData(nullptr),
64 fAxes(nullptr)
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(nullptr),
84 fData(&data),
85 fAxes(nullptr)
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 = nullptr;
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 nullptr;
153 if (fPOI.empty()) return nullptr;
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 == nullptr);
163 bool usePriorPdf = (fPriorPdf != nullptr);
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 std::unique_ptr<RooArgSet> constrainedParams{prodPdf->getParameters(*fData)};
174 std::unique_ptr<RooAbsReal> nll{prodPdf->createNLL(*fData, Constrain(*constrainedParams),ConditionalObservables(fConditionalObs),GlobalObservables(fGlobalObs))};
175
176 std::unique_ptr<RooArgSet> params{nll->getParameters(*fData)};
177 RemoveConstantParameters(&*params);
178 if (fNumBins > 0) {
179 SetBins(*params, fNumBins);
181 if (dynamic_cast<PdfProposal*>(fPropFunc)) {
182 std::unique_ptr<RooArgSet> proposalVars{((PdfProposal*)fPropFunc)->GetPdf()->
183 getParameters((RooAbsData*)nullptr)};
184 SetBins(*proposalVars, fNumBins);
185 }
186 }
187
189 mh.SetFunction(*nll);
192 mh.SetParameters(*params);
196
197 MarkovChain* chain = mh.ConstructChain();
198
199 TString name = TString("MCMCInterval_") + TString(GetName() );
200 MCMCInterval* interval = new MCMCInterval(name, fPOI, *chain);
201 if (fAxes != nullptr)
202 interval->SetAxes(*fAxes);
203 if (fNumBurnInSteps > 0)
205 interval->SetUseKeys(fUseKeys);
210 if (fEpsilon >= 0)
211 interval->SetEpsilon(fEpsilon);
212 if (fDelta >= 0)
213 interval->SetDelta(fDelta);
214 interval->SetConfidenceLevel(1.0 - fSize);
215
216 if (useDefaultPropFunc) delete fPropFunc;
217 if (usePriorPdf) delete prodPdf;
218
219 return interval;
220}
#define a(i)
Definition RSha256.hxx:99
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
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
RooFit::OwningPtr< 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...
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
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.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:162
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
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.
virtual void SetUseSparseHist(bool useSparseHist)
set whether to use a sparse histogram.
virtual void SetDelta(double delta)
kbelasco: The inner-workings of the class really should not be exposed like this in a comment,...
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
virtual void SetUseKeys(bool useKeys)
set whether to use kernel estimation to determine the interval
virtual void SetEpsilon(double epsilon)
set the acceptable level or error for Keys interval determination
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
virtual void SetLeftSideTailFraction(double a)
set the left-side tail fraction for a tail-fraction interval
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:35
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return nullptr if not existing)
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return nullptr if not existing)
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:139
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 JSONIO.h:26
@ InputArguments
Namespace for the RooStats classes.
Definition Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)