Logo ROOT  
Reference Guide
MCMCCalculator.h
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#ifndef ROOSTATS_MCMCCalculator
13#define ROOSTATS_MCMCCalculator
14
15#include "Rtypes.h"
16
17#include "TObject.h"
18#include "RooAbsPdf.h"
19#include "RooAbsData.h"
20#include "RooArgSet.h"
21#include "RooArgList.h"
25
26
27namespace RooStats {
28
29 class ModelConfig;
30
31 class MCMCCalculator : public IntervalCalculator, public TNamed {
32
33 public:
34 /// default constructor
36
37 /// Constructor for automatic configuration with basic settings and a
38 /// ModelConfig. Uses a UniformProposal, 10,000 iterations, 40 burn in
39 /// steps, 50 bins for each RooRealVar, determines interval by histogram,
40 /// and finds a 95% confidence interval. Any of these basic settings can
41 /// be overridden by calling one of the Set...() methods.
42 MCMCCalculator(RooAbsData& data, const ModelConfig& model);
43
44 virtual ~MCMCCalculator() {}
45
46 /// Main interface to get a ConfInterval
47 virtual MCMCInterval* GetInterval() const;
48
49 /// Get the size of the test (eg. rate of Type I error)
50 virtual Double_t Size() const {return fSize;}
51 /// Get the Confidence level for the test
52 virtual Double_t ConfidenceLevel() const {return 1.-fSize;}
53
54 virtual void SetModel(const ModelConfig & model);
55
56 /// Set the DataSet if not already there
57 virtual void SetData(RooAbsData& data) { fData = &data; }
58
59 /// Set the Pdf if not already there
60 virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; }
61
62 /// Set the Prior Pdf if not already there
63 virtual void SetPriorPdf(RooAbsPdf& pdf) { fPriorPdf = &pdf; }
64
65 /// specify the parameters of interest in the interval
66 virtual void SetParameters(const RooArgSet& set) { fPOI.removeAll(); fPOI.add(set); }
67
68 /// specify the parameters to store in the Markov chain
69 /// By default all the parameters are stored
70 virtual void SetChainParameters(const RooArgSet & set) { fChainParams.removeAll(); fChainParams.add(set); }
71
72 /// specify the nuisance parameters (eg. the rest of the parameters)
74
75 /// set the conditional observables which will be used when creating the NLL
76 /// so the pdf's will not be normalized on the conditional observables when computing the NLL
78
79 /// set the global observables which will be used when creating the NLL
80 /// so the constraint pdf's will be normalized correctly on the global observables when computing the NLL
82
83 /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
84 virtual void SetTestSize(Double_t size) {fSize = size;}
85
86 /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
87 virtual void SetConfidenceLevel(Double_t cl) {fSize = 1.-cl;}
88
89 /// set the proposal function for suggesting new points for the MCMC
90 virtual void SetProposalFunction(ProposalFunction& proposalFunction)
91 { fPropFunc = &proposalFunction; }
92
93 /// set the number of iterations to run the metropolis algorithm
94 virtual void SetNumIters(Int_t numIters)
95 { fNumIters = numIters; }
96
97 /// set the number of steps in the chain to discard as burn-in,
98 /// starting from the first
99 virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
100 { fNumBurnInSteps = numBurnInSteps; }
101
102 /// set the number of bins to create for each axis when constructing the interval
103 virtual void SetNumBins(Int_t numBins) { fNumBins = numBins; }
104 /// set which variables to put on each axis
105 virtual void SetAxes(RooArgList& axes)
106 { fAxes = &axes; }
107 /// set whether to use kernel estimation to determine the interval
108 virtual void SetUseKeys(Bool_t useKeys) { fUseKeys = useKeys; }
109 /// set whether to use sparse histogram (if using histogram at all)
110 virtual void SetUseSparseHist(Bool_t useSparseHist)
111 { fUseSparseHist = useSparseHist; }
112
113 /// set what type of interval to have the MCMCInterval represent
114 virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
115 { fIntervalType = intervalType; }
116
117 /// Set the left side tail fraction. This will automatically configure the
118 /// MCMCInterval to find a tail-fraction interval.
119 /// Note: that `a' must be in the range 0 <= a <= 1
120 /// or the user will be notified of the error
121 virtual void SetLeftSideTailFraction(Double_t a);
122
123 /// Set the desired level of confidence-level accuracy for Keys interval
124 /// determination.
125 //
126 /// When determining the cutoff PDF height that gives the
127 /// desired confidence level (C_d), the algorithm will consider acceptable
128 /// any found confidence level c such that Abs(c - C_d) < epsilon.
129 ///
130 /// Any value of this "epsilon" > 0 is considered acceptable, though it is
131 /// advisable to not use a value too small, because the integration of the
132 /// Keys PDF sometimes does not have extremely high accuracy.
134 {
135 if (epsilon < 0)
136 coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
137 << "negative epsilon value" << std::endl;
138 else
140 }
141
142 /// When the shortest interval using Keys PDF could not be found to have
143 /// the desired confidence level +/- the accuracy (see
144 /// SetKeysConfidenceAccuracy()), the interval determination algorithm
145 /// will have to terminate with an unsatisfactory confidence level when
146 /// the bottom and top of the cutoff search range are very close to being
147 /// equal. This scenario comes into play when there seems to be an error
148 /// in the accuracy of the Keys PDF integration, so the search range
149 /// continues to shrink without converging to a cutoff value that will
150 /// give an acceptable confidence level. To choose how small to allow the
151 /// search range to be before terminating, set the fraction delta such
152 /// that the search will terminate when topCutoff (a) and bottomCutoff (b)
153 /// satisfy this condition:
154 ///
155 /// TMath::Abs(a - b) < TMath::Abs(delta * (a + b)/2)
157 {
158 if (delta < 0.)
159 coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
160 << "negative delta value" << std::endl;
161 else
162 fDelta = delta;
163 }
164
165 protected:
166
167 Double_t fSize; // size of the test (eg. specified rate of Type I error)
168 RooArgSet fPOI; // parameters of interest for interval
169 RooArgSet fNuisParams; // nuisance parameters for interval (not really used)
170 RooArgSet fChainParams; // parameters to store in the chain (if not specified they are all of them )
171 RooArgSet fConditionalObs; // conditional observables
172 RooArgSet fGlobalObs; // global observables
173 mutable ProposalFunction* fPropFunc; // Proposal function for MCMC integration
174 RooAbsPdf * fPdf; // pointer to common PDF (owned by the workspace)
175 RooAbsPdf * fPriorPdf; // pointer to prior PDF (owned by the workspace)
176 RooAbsData * fData; // pointer to the data (owned by the workspace)
177 Int_t fNumIters; // number of iterations to run metropolis algorithm
178 Int_t fNumBurnInSteps; // number of iterations to discard as burn-in, starting from the first
179 Int_t fNumBins; // set the number of bins to create for each
180 // axis when constructing the interval
181 RooArgList * fAxes; // which variables to put on each axis
182 Bool_t fUseKeys; // whether to use kernel estimation to determine interval
183 Bool_t fUseSparseHist; // whether to use sparse histogram (if using hist at all)
184 Double_t fLeftSideTF; // left side tail-fraction for interval
185 Double_t fEpsilon; // acceptable error for Keys interval determination
186
187 Double_t fDelta; // acceptable error for Keys cutoffs being equal
188 // topCutoff (a) considered == bottomCutoff (b) iff
189 // (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2));
190 // Theoretically, the Abs is not needed here, but
191 // floating-point arithmetic does not always work
192 // perfectly, and the Abs doesn't hurt
193 enum MCMCInterval::IntervalType fIntervalType; // type of interval to find
194
195 void SetupBasicUsage();
196 void SetBins(const RooAbsCollection& coll, Int_t numBins) const
197 {
198 TIterator* it = coll.createIterator();
199 RooAbsArg* r;
200 while ((r = (RooAbsArg*)it->Next()) != NULL)
201 if (dynamic_cast<RooRealVar*>(r))
202 ((RooRealVar*)r)->setBins(numBins);
203 delete it;
204 }
205
206 ClassDef(MCMCCalculator,4) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
207 };
208}
209
210
211#endif
ROOT::R::TRInterface & r
Definition: Object.C:4
#define coutE(a)
Definition: RooMsgService.h:33
double Double_t
Definition: RtypesCore.h:57
#define ClassDef(name, id)
Definition: Rtypes.h:322
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
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
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
IntervalCalculator is an interface class for a tools which produce RooStats ConfIntervals.
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
virtual void SetPriorPdf(RooAbsPdf &pdf)
Set the Prior Pdf if not already there.
virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
set the number of steps in the chain to discard as burn-in, starting from the first
enum MCMCInterval::IntervalType fIntervalType
virtual void SetNumBins(Int_t numBins)
set the number of bins to create for each axis when constructing the interval
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 SetUseKeys(Bool_t useKeys)
set whether to use kernel estimation to determine the interval
virtual void SetLeftSideTailFraction(Double_t a)
Set the left side tail fraction.
virtual void SetKeysConfidenceAccuracy(Double_t epsilon)
Set the desired level of confidence-level accuracy for Keys interval determination.
virtual void SetNumIters(Int_t numIters)
set the number of iterations to run the metropolis algorithm
virtual void SetNuisanceParameters(const RooArgSet &set)
specify the nuisance parameters (eg. the rest of the parameters)
virtual void SetParameters(const RooArgSet &set)
specify the parameters of interest in the interval
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 SetKeysTerminationThreshold(Double_t delta)
When the shortest interval using Keys PDF could not be found to have the desired confidence level +/-...
virtual Double_t Size() const
Get the size of the test (eg. rate of Type I error)
virtual void SetData(RooAbsData &data)
Set the DataSet if not already there.
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
virtual void SetModel(const ModelConfig &model)
Set the Model.
virtual void SetAxes(RooArgList &axes)
set which variables to put on each axis
virtual void SetPdf(RooAbsPdf &pdf)
Set the Pdf if not already there.
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
virtual void SetConditionalObservables(const RooArgSet &set)
set the conditional observables which will be used when creating the NLL so the pdf's will not be nor...
virtual void SetProposalFunction(ProposalFunction &proposalFunction)
set the proposal function for suggesting new points for the MCMC
virtual void SetUseSparseHist(Bool_t useSparseHist)
set whether to use sparse histogram (if using histogram at all)
virtual void SetGlobalObservables(const RooArgSet &set)
set the global observables which will be used when creating the NLL so the constraint pdf's will be n...
virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
set what type of interval to have the MCMCInterval represent
ProposalFunction * fPropFunc
virtual void SetChainParameters(const RooArgSet &set)
specify the parameters to store in the Markov chain By default all the parameters are stored
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:19
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617