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 ~MCMCCalculator() override {}
45
46 /// Main interface to get a ConfInterval
47 MCMCInterval* GetInterval() const override;
48
49 /// Get the size of the test (eg. rate of Type I error)
50 double Size() const override {return fSize;}
51 /// Get the Confidence level for the test
52 double ConfidenceLevel() const override {return 1.-fSize;}
53
54 void SetModel(const ModelConfig & model) override;
55
56 /// Set the DataSet if not already there
57 void SetData(RooAbsData& data) override { 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 void SetTestSize(double size) override {fSize = size;}
85
86 /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
87 void SetConfidenceLevel(double cl) override {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 useKeys) { fUseKeys = useKeys; }
109 /// set whether to use sparse histogram (if using histogram at all)
110 virtual void SetUseSparseHist(bool 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 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)
156 virtual void SetKeysTerminationThreshold(double delta)
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 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 fUseKeys; ///< whether to use kernel estimation to determine interval
183 bool fUseSparseHist; ///< whether to use sparse histogram (if using hist at all)
184 double fLeftSideTF; ///< left side tail-fraction for interval
185 double fEpsilon; ///< acceptable error for Keys interval determination
186
187 double 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 for (auto *r : dynamic_range_cast<RooRealVar *>(coll)){
199 if (r) {
200 r->setBins(numBins);
201 }
202 }
203 }
204
205 ClassDefOverride(MCMCCalculator,4) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
206 };
207}
208
209
210#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
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.
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:59
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
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
double epsilon
Definition triangle.c:618