Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 "TNamed.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.
43
44 /// Main interface to get a ConfInterval
45 MCMCInterval* GetInterval() const override;
46
47 /// Get the size of the test (eg. rate of Type I error)
48 double Size() const override {return fSize;}
49 /// Get the Confidence level for the test
50 double ConfidenceLevel() const override {return 1.-fSize;}
51
52 void SetModel(const ModelConfig & model) override;
53
54 /// Set the DataSet if not already there
55 void SetData(RooAbsData& data) override { fData = &data; }
56
57 /// Set the Pdf if not already there
58 virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; }
59
60 /// Set the Prior Pdf if not already there
61 virtual void SetPriorPdf(RooAbsPdf& pdf) { fPriorPdf = &pdf; }
62
63 /// specify the parameters of interest in the interval
64 virtual void SetParameters(const RooArgSet& set) { fPOI.removeAll(); fPOI.add(set); }
65
66 /// specify the parameters to store in the Markov chain
67 /// By default all the parameters are stored
68 virtual void SetChainParameters(const RooArgSet & set) { fChainParams.removeAll(); fChainParams.add(set); }
69
70 /// specify the nuisance parameters (eg. the rest of the parameters)
72
73 /// set the conditional observables which will be used when creating the NLL
74 /// so the pdf's will not be normalized on the conditional observables when computing the NLL
76
77 /// set the global observables which will be used when creating the NLL
78 /// so the constraint pdf's will be normalized correctly on the global observables when computing the NLL
80
81 /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
82 void SetTestSize(double size) override {fSize = size;}
83
84 /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
85 void SetConfidenceLevel(double cl) override {fSize = 1.-cl;}
86
87 /// set the proposal function for suggesting new points for the MCMC
88 virtual void SetProposalFunction(ProposalFunction& proposalFunction)
89 { fPropFunc = &proposalFunction; }
90
91 /// set the number of iterations to run the metropolis algorithm
92 virtual void SetNumIters(Int_t numIters)
93 { fNumIters = numIters; }
94
95 /// set the number of steps in the chain to discard as burn-in,
96 /// starting from the first
97 virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
98 { fNumBurnInSteps = numBurnInSteps; }
99
100 /// set the number of bins to create for each axis when constructing the interval
101 virtual void SetNumBins(Int_t numBins) { fNumBins = numBins; }
102 /// set which variables to put on each axis
103 virtual void SetAxes(RooArgList& axes)
104 { fAxes = &axes; }
105 /// set whether to use kernel estimation to determine the interval
106 virtual void SetUseKeys(bool useKeys) { fUseKeys = useKeys; }
107 /// set whether to use sparse histogram (if using histogram at all)
108 virtual void SetUseSparseHist(bool useSparseHist)
109 { fUseSparseHist = useSparseHist; }
110
111 /// set what type of interval to have the MCMCInterval represent
112 virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
113 { fIntervalType = intervalType; }
114
115 /// Set the left side tail fraction. This will automatically configure the
116 /// MCMCInterval to find a tail-fraction interval.
117 /// Note: that `a' must be in the range 0 <= a <= 1
118 /// or the user will be notified of the error
119 virtual void SetLeftSideTailFraction(double a);
120
121 /// Set the desired level of confidence-level accuracy for Keys interval
122 /// determination.
123 //
124 /// When determining the cutoff PDF height that gives the
125 /// desired confidence level (C_d), the algorithm will consider acceptable
126 /// any found confidence level c such that Abs(c - C_d) < epsilon.
127 ///
128 /// Any value of this "epsilon" > 0 is considered acceptable, though it is
129 /// advisable to not use a value too small, because the integration of the
130 /// Keys PDF sometimes does not have extremely high accuracy.
131 virtual void SetKeysConfidenceAccuracy(double epsilon)
132 {
133 if (epsilon < 0) {
134 coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
135 << "negative epsilon value" << std::endl;
136 } else {
137 fEpsilon = epsilon;
138 }
139 }
140
141 /// When the shortest interval using Keys PDF could not be found to have
142 /// the desired confidence level +/- the accuracy (see
143 /// SetKeysConfidenceAccuracy()), the interval determination algorithm
144 /// will have to terminate with an unsatisfactory confidence level when
145 /// the bottom and top of the cutoff search range are very close to being
146 /// equal. This scenario comes into play when there seems to be an error
147 /// in the accuracy of the Keys PDF integration, so the search range
148 /// continues to shrink without converging to a cutoff value that will
149 /// give an acceptable confidence level. To choose how small to allow the
150 /// search range to be before terminating, set the fraction delta such
151 /// that the search will terminate when topCutoff (a) and bottomCutoff (b)
152 /// satisfy this condition:
153 ///
154 /// std::abs(a - b) < std::abs(delta * (a + b)/2)
155 virtual void SetKeysTerminationThreshold(double delta)
156 {
157 if (delta < 0.) {
158 coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
159 << "negative delta value" << std::endl;
160 } else {
161 fDelta = delta;
162 }
163 }
164
165 protected:
166 double fSize = -1; ///< size of the test (eg. specified rate of Type I error)
167 RooArgSet fPOI; ///< parameters of interest for interval
168 RooArgSet fNuisParams; ///< nuisance parameters for interval (not really used)
169 RooArgSet fChainParams; ///< parameters to store in the chain (if not specified they are all of them )
170 RooArgSet fConditionalObs; ///< conditional observables
171 RooArgSet fGlobalObs; ///< global observables
172 mutable ProposalFunction* fPropFunc; ///< Proposal function for MCMC integration
173 RooAbsPdf * fPdf; ///< pointer to common PDF (owned by the workspace)
174 RooAbsPdf * fPriorPdf; ///< pointer to prior PDF (owned by the workspace)
175 RooAbsData * fData; ///< pointer to the data (owned by the workspace)
176 Int_t fNumIters = 0; ///< number of iterations to run metropolis algorithm
177 Int_t fNumBurnInSteps = 0; ///< number of iterations to discard as burn-in, starting from the first
178 Int_t fNumBins = 0; ///< set the number of bins to create for each
179 ///< axis when constructing the interval
180 RooArgList * fAxes; ///< which variables to put on each axis
181 bool fUseKeys = false; ///< whether to use kernel estimation to determine interval
182 bool fUseSparseHist = false; ///< whether to use sparse histogram (if using hist at all)
183 double fLeftSideTF = -1; ///< left side tail-fraction for interval
184 double fEpsilon = -1; ///< acceptable error for Keys interval determination
185
186 double fDelta = -1; ///< acceptable error for Keys cutoffs being equal
187 ///< topCutoff (a) considered == bottomCutoff (b) iff
188 ///< (std::abs(a - b) < std::abs(fDelta * (a + b)/2));
189 ///< Theoretically, the Abs is not needed here, but
190 ///< floating-point arithmetic does not always work
191 ///< perfectly, and the Abs doesn't hurt
193
194 void SetupBasicUsage();
195 void SetBins(const RooAbsCollection &coll, Int_t numBins) const
196 {
197 for (auto *r : dynamic_range_cast<RooRealVar *>(coll)){
198 if (r) {
199 r->setBins(numBins);
200 }
201 }
202 }
203
204 ClassDefOverride(MCMCCalculator,4) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
205 };
206}
207
208
209#endif
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutE(a)
#define ClassDefOverride(name, id)
Definition Rtypes.h:341
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Abstract container object that can hold multiple RooAbsArg objects.
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
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:55
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.
RooArgSet fChainParams
parameters to store in the chain (if not specified they are all of them )
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 SetUseSparseHist(bool useSparseHist)
set whether to use sparse histogram (if using histogram at all)
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
virtual void SetNumBins(Int_t numBins)
set the number of bins to create for each axis when constructing the 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 (std::...
MCMCCalculator()
default constructor
void SetModel(const ModelConfig &model) override
Set the Model.
virtual void SetLeftSideTailFraction(double a)
Set the left side tail fraction.
virtual void SetKeysConfidenceAccuracy(double 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)
RooArgSet fNuisParams
nuisance parameters for interval (not really used)
void SetConfidenceLevel(double cl) override
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
virtual void SetParameters(const RooArgSet &set)
specify the parameters of interest in the interval
double Size() const override
Get the size of the test (eg. rate of Type I error)
void SetData(RooAbsData &data) override
Set the DataSet if not already there.
Int_t fNumIters
number of iterations to run metropolis algorithm
virtual void SetKeysTerminationThreshold(double delta)
When the shortest interval using Keys PDF could not be found to have the desired confidence level +/-...
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
virtual void SetAxes(RooArgList &axes)
set which variables to put on each axis
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
virtual void SetPdf(RooAbsPdf &pdf)
Set the Pdf if not already there.
RooArgList * fAxes
which variables to put on each axis
RooArgSet fGlobalObs
global observables
MCMCInterval * GetInterval() const override
Main interface to get a ConfInterval.
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 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
RooAbsPdf * fPriorPdf
pointer to prior PDF (owned by the workspace)
ProposalFunction * fPropFunc
Proposal function for MCMC integration.
virtual void SetUseKeys(bool useKeys)
set whether to use kernel estimation to determine the interval
double ConfidenceLevel() const override
Get the Confidence level for the test.
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.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
RooStats::ModelConfig ModelConfig
Namespace for the RooStats classes.
Definition Asimov.h:19