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 std::endl;
55
56////////////////////////////////////////////////////////////////////////////////
57/// default constructor
58
60 fPropFunc(nullptr),
61 fPdf(nullptr),
62 fPriorPdf(nullptr),
63 fData(nullptr),
64 fAxes(nullptr)
65{
66}
67
68////////////////////////////////////////////////////////////////////////////////
69/// constructor from a Model Config with a basic settings package configured
70/// by SetupBasicUsage()
71
73 fPropFunc(nullptr),
74 fData(&data),
75 fAxes(nullptr)
76{
77 SetModel(model);
79}
80
82 // set the model
83 fPdf = model.GetPdf();
84 fPriorPdf = model.GetPriorPdf();
89 if (model.GetParametersOfInterest())
91 if (model.GetNuisanceParameters())
93 if (model.GetConditionalObservables())
95 if (model.GetGlobalObservables())
96 fGlobalObs.add( *(model.GetGlobalObservables() ) );
97
98}
99
100////////////////////////////////////////////////////////////////////////////////
101/// Constructor for automatic configuration with basic settings. Uses a
102/// UniformProposal, 10,000 iterations, 40 burn in steps, 50 bins for each
103/// RooRealVar, determines interval by histogram. Finds a 95% confidence
104/// interval.
105
107{
108 fPropFunc = nullptr;
109 fNumIters = 10000;
110 fNumBurnInSteps = 40;
111 fNumBins = 50;
112 fUseKeys = false;
113 fUseSparseHist = false;
114 SetTestSize(0.05);
116 fLeftSideTF = -1;
117 fEpsilon = -1;
118 fDelta = -1;
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
124{
125 if (a < 0 || a > 1) {
126 coutE(InputArguments) << "MCMCCalculator::SetLeftSideTailFraction: "
127 << "Fraction must be in the range [0, 1]. "
128 << a << "is not allowed." << endl;
129 return;
130 }
131
132 fLeftSideTF = a;
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// Main interface to get a RooStats::ConfInterval.
138
140{
141
142 if (!fData || !fPdf ) return nullptr;
143 if (fPOI.empty()) return nullptr;
144
145 if (fSize < 0) {
146 coutE(InputArguments) << "MCMCCalculator::GetInterval: "
147 << "Test size/Confidence level not set. Returning nullptr." << endl;
148 return nullptr;
149 }
150
151 // if a proposal function has not been specified create a default one
152 bool useDefaultPropFunc = (fPropFunc == nullptr);
153 bool usePriorPdf = (fPriorPdf != nullptr);
154 if (useDefaultPropFunc) fPropFunc = new UniformProposal();
155
156 // if prior is given create product
157 RooAbsPdf * prodPdf = fPdf;
158 if (usePriorPdf) {
159 TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );
160 prodPdf = new RooProdPdf(prodName,prodName,RooArgList(*fPdf,*fPriorPdf) );
161 }
162
163 std::unique_ptr<RooArgSet> constrainedParams{prodPdf->getParameters(*fData)};
164 std::unique_ptr<RooAbsReal> nll{prodPdf->createNLL(*fData, Constrain(*constrainedParams),ConditionalObservables(fConditionalObs),GlobalObservables(fGlobalObs))};
165
166 std::unique_ptr<RooArgSet> params{nll->getParameters(*fData)};
167 RemoveConstantParameters(&*params);
168 if (fNumBins > 0) {
169 SetBins(*params, fNumBins);
171 if (dynamic_cast<PdfProposal*>(fPropFunc)) {
172 std::unique_ptr<RooArgSet> proposalVars{(static_cast<PdfProposal*>(fPropFunc))->GetPdf()->
173 getParameters((RooAbsData*)nullptr)};
174 SetBins(*proposalVars, fNumBins);
175 }
176 }
177
179 mh.SetFunction(*nll);
182 mh.SetParameters(*params);
186
187 MarkovChain* chain = mh.ConstructChain();
188
189 TString name = TString("MCMCInterval_") + TString(GetName() );
190 MCMCInterval* interval = new MCMCInterval(name, fPOI, *chain);
191 if (fAxes != nullptr)
192 interval->SetAxes(*fAxes);
193 if (fNumBurnInSteps > 0)
195 interval->SetUseKeys(fUseKeys);
200 if (fEpsilon >= 0)
201 interval->SetEpsilon(fEpsilon);
202 if (fDelta >= 0)
203 interval->SetDelta(fDelta);
204 interval->SetConfidenceLevel(1.0 - fSize);
205
206 if (useDefaultPropFunc) delete fPropFunc;
207 if (usePriorPdf) delete prodPdf;
208
209 return interval;
210}
#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.
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:163
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 topCutoff (a) considered == bottomCutoff (b) iff (TMath...
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)