Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
MCMCInterval.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_MCMCInterval
13#define RooStats_MCMCInterval
14
15#include "Rtypes.h"
16
18#include "RooArgSet.h"
19#include "RooArgList.h"
20#include "RooMsgService.h"
22
23#include <vector>
24
25class RooNDKeysPdf;
26class RooProduct;
27
28
29namespace RooStats {
30
31 class Heaviside;
32
33 class MCMCInterval : public ConfInterval {
34
35
36 public:
37
38 /// default constructor
39 explicit MCMCInterval(const char *name = nullptr);
40
41 /// constructor from parameter of interest and Markov chain object
42 MCMCInterval(const char* name, const RooArgSet& parameters,
43 MarkovChain& chain);
44
45 enum {DEFAULT_NUM_BINS = 50};
47
48 ~MCMCInterval() override;
49
50 /// determine whether this point is in the confidence interval
51 bool IsInInterval(const RooArgSet& point) const override;
52
53 /// set the desired confidence level (see GetActualConfidenceLevel())
54 /// Note: calling this function triggers the algorithm that determines
55 /// the interval, so call this after initializing all other aspects
56 /// of this IntervalCalculator
57 /// Also, calling this function again with a different confidence level
58 /// re-triggers the calculation of the interval
59 void SetConfidenceLevel(double cl) override;
60
61 /// get the desired confidence level (see GetActualConfidenceLevel())
62 double ConfidenceLevel() const override {return fConfidenceLevel;}
63
64 /// return a set containing the parameters of this interval
65 /// the caller owns the returned RooArgSet*
66 RooArgSet* GetParameters() const override;
67
68 /// get the cutoff bin height for being considered in the
69 /// confidence interval
70 virtual double GetHistCutoff();
71
72 /// get the cutoff RooNDKeysPdf value for being considered in the
73 /// confidence interval
74 virtual double GetKeysPdfCutoff();
75 ///virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
76
77 /// get the actual value of the confidence level for this interval.
78 virtual double GetActualConfidenceLevel();
79
80 /// whether the specified confidence level is a floor for the actual
81 /// confidence level (strict), or a ceiling (not strict)
82 virtual void SetHistStrict(bool isHistStrict)
83 { fIsHistStrict = isHistStrict; }
84
85 /// check if parameters are correct. (dummy implementation to start)
86 bool CheckParameters(const RooArgSet& point) const override;
87
88 /// Set the parameters of interest for this interval
89 /// and change other internal data members accordingly
90 virtual void SetParameters(const RooArgSet& parameters);
91
92 /// Set the MarkovChain that this interval is based on
93 virtual void SetChain(MarkovChain& chain) { fChain = &chain; }
94
95 /// Set which parameters go on which axis. The first list element
96 /// goes on the x axis, second (if it exists) on y, third (if it
97 /// exists) on z, etc
98 virtual void SetAxes(RooArgList& axes);
99
100 /// return a list of RooRealVars representing the axes
101 /// you own the returned RooArgList
103 {
104 RooArgList* axes = new RooArgList();
105 for (Int_t i = 0; i < fDimension; i++)
106 axes->addClone(*fAxes[i]);
107 return axes;
108 }
109
110 /// get the lowest value of param that is within the confidence interval
111 virtual double LowerLimit(RooRealVar& param);
112
113 /// determine lower limit of the lower confidence interval
114 virtual double LowerLimitTailFraction(RooRealVar& param);
115
116 /// get the lower limit of param in the shortest confidence interval
117 /// Note that this works better for some distributions (ones with exactly
118 /// one maximum) than others, and sometimes has little value.
119 virtual double LowerLimitShortest(RooRealVar& param);
120
121 /// determine lower limit in the shortest interval by using keys pdf
122 virtual double LowerLimitByKeys(RooRealVar& param);
123
124 /// determine lower limit using histogram
125 virtual double LowerLimitByHist(RooRealVar& param);
126
127 /// determine lower limit using histogram
128 virtual double LowerLimitBySparseHist(RooRealVar& param);
129
130 /// determine lower limit using histogram
131 virtual double LowerLimitByDataHist(RooRealVar& param);
132
133 /// get the highest value of param that is within the confidence interval
134 virtual double UpperLimit(RooRealVar& param);
135
136 /// determine upper limit of the lower confidence interval
137 virtual double UpperLimitTailFraction(RooRealVar& param);
138
139 /// get the upper limit of param in the confidence interval
140 /// Note that this works better for some distributions (ones with exactly
141 /// one maximum) than others, and sometimes has little value.
142 virtual double UpperLimitShortest(RooRealVar& param);
143
144 /// determine upper limit in the shortest interval by using keys pdf
145 virtual double UpperLimitByKeys(RooRealVar& param);
146
147 /// determine upper limit using histogram
148 virtual double UpperLimitByHist(RooRealVar& param);
149
150 /// determine upper limit using histogram
151 virtual double UpperLimitBySparseHist(RooRealVar& param);
152
153 /// determine upper limit using histogram
154 virtual double UpperLimitByDataHist(RooRealVar& param);
155
156 /// Determine the approximate maximum value of the Keys PDF
157 double GetKeysMax();
158
159 /// set the number of steps in the chain to discard as burn-in,
160 /// starting from the first
161 virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
162 { fNumBurnInSteps = numBurnInSteps; }
163
164 /// set whether to use kernel estimation to determine the interval
165 virtual void SetUseKeys(bool useKeys) { fUseKeys = useKeys; }
166
167 /// set whether to use a sparse histogram. you MUST also call
168 /// SetUseKeys(false) to use a histogram.
169 virtual void SetUseSparseHist(bool useSparseHist)
170 { fUseSparseHist = useSparseHist; }
171
172 /// get whether we used kernel estimation to determine the interval
173 virtual bool GetUseKeys() { return fUseKeys; }
174
175 /// get the number of steps in the chain to discard as burn-in,
176
177 /// get the number of steps in the chain to discard as burn-in,
178 /// starting from the first
180
181 /// set the number of bins to use (same for all axes, for now)
182 ///virtual void SetNumBins(Int_t numBins);
183
184 /// Get a clone of the histogram of the posterior
185 virtual TH1* GetPosteriorHist();
186
187 /// Get a clone of the keys pdf of the posterior
189
190 /// Get a clone of the (keyspdf * heaviside) product of the posterior
192
193 /// Get the number of parameters of interest in this interval
194 virtual Int_t GetDimension() const { return fDimension; }
195
196 /// Get the markov chain on which this interval is based
197 /// You do not own the returned MarkovChain*
198 virtual const MarkovChain* GetChain() { return fChain; }
199
200 /// Get a clone of the markov chain on which this interval is based
201 /// as a RooDataSet. You own the returned RooDataSet*
203 { return fChain->GetAsDataSet(whichVars); }
204
205 /// Get the markov chain on which this interval is based
206 /// as a RooDataSet. You do not own the returned RooDataSet*
208 { return fChain->GetAsConstDataSet(); }
209
210 /// Get a clone of the markov chain on which this interval is based
211 /// as a RooDataHist. You own the returned RooDataHist*
213 { return fChain->GetAsDataHist(whichVars); }
214
215 /// Get a clone of the markov chain on which this interval is based
216 /// as a THnSparse. You own the returned THnSparse*
217 virtual THnSparse* GetChainAsSparseHist(RooArgSet* whichVars = nullptr)
218 { return fChain->GetAsSparseHist(whichVars); }
219
220 /// Get a clone of the NLL variable from the markov chain
221 virtual RooRealVar* GetNLLVar() const
222 { return fChain->GetNLLVar(); }
223
224 /// Get a clone of the weight variable from the markov chain
225 virtual RooRealVar* GetWeightVar() const
226 { return fChain->GetWeightVar(); }
227
228 /// set the acceptable level or error for Keys interval determination
229 virtual void SetEpsilon(double epsilon)
230 {
231 if (epsilon < 0) {
232 coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
233 << "negative epsilon value" << std::endl;
234 } else {
235 fEpsilon = epsilon;
236 }
237 }
238
239 /// Set the type of interval to find. This will only have an effect for
240 /// 1-D intervals. If is more than 1 parameter of interest, then a
241 /// "shortest" interval will always be used, since it generalizes directly
242 /// to N dimensions
243 virtual void SetIntervalType(enum IntervalType intervalType)
244 { fIntervalType = intervalType; }
246
247 /// Return the type of this interval
248 virtual enum IntervalType GetIntervalType() { return fIntervalType; }
249
250 /// set the left-side tail fraction for a tail-fraction interval
251 virtual void SetLeftSideTailFraction(double a) {
253 fLeftSideTF = a;
254 }
255
256 /// kbelasco: The inner-workings of the class really should not be exposed
257 /// like this in a comment, but it seems to be the only way to give
258 /// the user any control over this process, if they desire it
259 ///
260 /// Set the fraction delta such that
261 /// topCutoff (a) is considered == bottomCutoff (b) iff
262 /// (std::abs(a - b) < std::abs(fDelta * (a + b)/2))
263 /// when determining the confidence interval by Keys
264 virtual void SetDelta(double delta)
265 {
266 if (delta < 0.) {
267 coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
268 << "negative delta value" << std::endl;
269 } else {
270 fDelta = delta;
271 }
272 }
273
274 private:
275 inline bool AcceptableConfLevel(double confLevel);
276 inline bool WithinDeltaFraction(double a, double b);
277
278 constexpr static const double DEFAULT_EPSILON = 0.01;
279 constexpr static const double DEFAULT_DELTA = 10e-6;
280
281 protected:
282 RooArgSet fParameters; ///< parameters of interest for this interval
283 MarkovChain *fChain = nullptr; ///< the markov chain
284 double fConfidenceLevel = 0.0; ///< Requested confidence level (eg. 0.95 for 95% CL)
285
286 RooDataHist *fDataHist = nullptr; ///< the binned Markov Chain data
287 THnSparse *fSparseHist = nullptr; ///< the binned Markov Chain data
288 double fHistConfLevel = 0.0; ///< the actual conf level determined by hist
289 double fHistCutoff = -1; ///< cutoff bin size to be in interval
290
291 RooNDKeysPdf *fKeysPdf = nullptr; ///< the kernel estimation pdf
292 RooProduct *fProduct = nullptr; ///< the (keysPdf * heaviside) product
293 Heaviside *fHeaviside = nullptr; ///< the Heaviside function
294 RooDataHist *fKeysDataHist = nullptr; ///< data hist representing product
295 RooRealVar *fCutoffVar = nullptr; ///< cutoff variable to use for integrating keys pdf
296 double fKeysConfLevel = 0.0; ///< the actual conf level determined by keys
297 double fKeysCutoff = -1; ///< cutoff keys pdf value to be in interval
298 double fFull = 0.0; ///< Value of intergral of fProduct
299
300 double fLeftSideTF = -1; ///< left side tail-fraction for interval
301 double fTFConfLevel = 0.0; ///< the actual conf level of tail-fraction interval
302 std::vector<Int_t> fVector; ///< vector containing the Markov chain data
303 double fVecWeight = 0; ///< sum of weights of all entries in fVector
304 double fTFLower; ///< lower limit of the tail-fraction interval
305 double fTFUpper; ///< upper limit of the tail-fraction interval
306
307 TH1 *fHist = nullptr; ///< the binned Markov Chain data
308
309 bool fUseKeys = false; ///< whether to use kernel estimation
310 bool fUseSparseHist = false; ///< whether to use sparse hist (vs. RooDataHist)
311 bool fIsHistStrict = true; ///< whether the specified confidence level is a
312 ///< floor for the actual confidence level (strict),
313 ///< or a ceiling (not strict) for determination by
314 ///< histogram
315 Int_t fDimension = 1; ///< number of variables
316 Int_t fNumBurnInSteps = 0; ///< number of steps to discard as burn in, starting
317 ///< from the first
318 RooRealVar **fAxes = nullptr; ///< array of pointers to RooRealVars representing
319 ///< the axes of the histogram
320 ///< fAxes[0] represents x-axis, [1] y, [2] z, etc
321
322 double fEpsilon = DEFAULT_EPSILON; ///< acceptable error for Keys interval determination
323
324 double fDelta = DEFAULT_DELTA; ///< topCutoff (a) considered == bottomCutoff (b) iff
325 ///< (std::abs(a - b) < std::abs(fDelta * (a + b)/2));
326 ///< Theoretically, the Abs is not needed here, but
327 ///< floating-point arithmetic does not always work
328 ///< perfectly, and the Abs doesn't hurt
330
331 // functions
332 virtual void DetermineInterval();
333 virtual void DetermineShortestInterval();
334 virtual void DetermineTailFractionInterval();
335 virtual void DetermineByHist();
336 virtual void DetermineBySparseHist();
337 virtual void DetermineByDataHist();
338 virtual void DetermineByKeys();
339 virtual void CreateHist();
340 virtual void CreateSparseHist();
341 virtual void CreateDataHist();
342 virtual void CreateKeysPdf();
343 virtual void CreateKeysDataHist();
344 virtual void CreateVector(RooRealVar* param);
345 inline virtual double CalcConfLevel(double cutoff, double full);
346
347 ClassDefOverride(MCMCInterval,1) // Concrete implementation of a ConfInterval based on MCMC calculation
348
349 };
350}
351
352#endif
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
#define coutE(a)
#define ClassDefOverride(name, id)
Definition Rtypes.h:341
char name[80]
Definition TGX11.cxx:110
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
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
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:39
Container class to hold unbinned data.
Definition RooDataSet.h:57
Generic N-dimensional implementation of a kernel estimation p.d.f.
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
Variable that can be changed from the outside.
Definition RooRealVar.h:37
ConfInterval is an interface class for a generic interval in the RooStats framework.
Represents the Heaviside function.
Definition Heaviside.h:21
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual void CreateDataHist()
virtual bool GetUseKeys()
get whether we used kernel estimation to determine the interval
virtual void DetermineByDataHist()
virtual void DetermineShortestInterval()
virtual double GetActualConfidenceLevel()
virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
double fKeysConfLevel
the actual conf level determined by keys
virtual void CreateVector(RooRealVar *param)
RooDataHist * fDataHist
the binned Markov Chain data
static constexpr const double DEFAULT_EPSILON
virtual double LowerLimitByDataHist(RooRealVar &param)
determine lower limit using histogram
double fDelta
topCutoff (a) considered == bottomCutoff (b) iff (std::abs(a - b) < std::abs(fDelta * (a + b)/2)); Th...
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
TH1 * fHist
the binned Markov Chain data
virtual double UpperLimitBySparseHist(RooRealVar &param)
determine upper limit using histogram
virtual void SetUseSparseHist(bool useSparseHist)
set whether to use a sparse histogram.
enum IntervalType fIntervalType
virtual double UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
RooRealVar * fCutoffVar
cutoff variable to use for integrating keys pdf
virtual void SetDelta(double delta)
kbelasco: The inner-workings of the class really should not be exposed like this in a comment,...
virtual RooFit::OwningPtr< RooDataSet > GetChainAsDataSet(RooArgSet *whichVars=nullptr)
Get a clone of the markov chain on which this interval is based as a RooDataSet.
virtual const RooDataSet * GetChainAsConstDataSet()
Get the markov chain on which this interval is based as a RooDataSet.
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void DetermineByKeys()
double fTFLower
lower limit of the tail-fraction interval
bool fUseSparseHist
whether to use sparse hist (vs. RooDataHist)
double fVecWeight
sum of weights of all entries in fVector
double fLeftSideTF
left side tail-fraction for interval
bool AcceptableConfLevel(double confLevel)
virtual void DetermineBySparseHist()
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
RooProduct * fProduct
the (keysPdf * heaviside) product
virtual double LowerLimitBySparseHist(RooRealVar &param)
determine lower limit using histogram
bool WithinDeltaFraction(double a, double b)
double fKeysCutoff
cutoff keys pdf value to be in interval
virtual double LowerLimitShortest(RooRealVar &param)
get the lower limit of param in the shortest confidence interval Note that this works better for some...
virtual double LowerLimitTailFraction(RooRealVar &param)
determine lower limit of the lower confidence interval
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
Int_t fDimension
number of variables
void SetConfidenceLevel(double cl) override
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
RooRealVar ** fAxes
array of pointers to RooRealVars representing the axes of the histogram fAxes[0] represents x-axis,...
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList
virtual RooRealVar * GetWeightVar() const
Get a clone of the weight variable from the markov chain.
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 DetermineTailFractionInterval()
double fTFConfLevel
the actual conf level of tail-fraction interval
double fHistConfLevel
the actual conf level determined by hist
virtual void CreateSparseHist()
Heaviside * fHeaviside
the Heaviside function
double fFull
Value of intergral of fProduct.
virtual Int_t GetNumBurnInSteps()
get the number of steps in the chain to discard as burn-in,
virtual double UpperLimitByDataHist(RooRealVar &param)
determine upper limit using histogram
virtual void DetermineInterval()
virtual RooRealVar * GetNLLVar() const
Get a clone of the NLL variable from the markov chain.
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 const MarkovChain * GetChain()
Get the markov chain on which this interval is based You do not own the returned MarkovChain*.
virtual double CalcConfLevel(double cutoff, double full)
virtual double LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
static constexpr const double DEFAULT_DELTA
virtual double LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
virtual double UpperLimitShortest(RooRealVar &param)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
THnSparse * fSparseHist
the binned Markov Chain data
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
virtual void SetHistStrict(bool isHistStrict)
whether the specified confidence level is a floor for the actual confidence level (strict),...
bool CheckParameters(const RooArgSet &point) const override
check if parameters are correct. (dummy implementation to start)
std::vector< Int_t > fVector
vector containing the Markov chain data
virtual double LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
double fTFUpper
upper limit of the tail-fraction interval
double fHistCutoff
cutoff bin size to be in interval
virtual void SetLeftSideTailFraction(double a)
set the left-side tail fraction for a tail-fraction interval
MarkovChain * fChain
the markov chain
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
virtual double UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
bool fIsHistStrict
whether the specified confidence level is a floor for the actual confidence level (strict),...
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual double UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual void CreateKeysPdf()
virtual double UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
Int_t fNumBurnInSteps
number of steps to discard as burn in, starting from the first
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
RooDataHist * fKeysDataHist
data hist representing product
virtual void SetShortestInterval()
RooArgSet fParameters
parameters of interest for this interval
bool IsInInterval(const RooArgSet &point) const override
determine whether this point is in the confidence interval
virtual void DetermineByHist()
double ConfidenceLevel() const override
get the desired confidence level (see GetActualConfidenceLevel())
virtual void SetChain(MarkovChain &chain)
Set the MarkovChain that this interval is based on.
virtual void SetParameters(const RooArgSet &parameters)
Set the parameters of interest for this interval and change other internal data members accordingly.
virtual Int_t GetDimension() const
Get the number of parameters of interest in this interval.
virtual void CreateKeysDataHist()
RooNDKeysPdf * fKeysPdf
the kernel estimation pdf
virtual THnSparse * GetChainAsSparseHist(RooArgSet *whichVars=nullptr)
Get a clone of the markov chain on which this interval is based as a THnSparse.
virtual RooFit::OwningPtr< RooDataHist > GetChainAsDataHist(RooArgSet *whichVars=nullptr)
Get a clone of the markov chain on which this interval is based as a RooDataHist.
double fEpsilon
acceptable error for Keys interval determination
virtual double GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual double GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Stores the steps in a Markov Chain of points.
Definition MarkovChain.h:30
virtual RooRealVar * GetWeightVar() const
get a clone of the weight variable
virtual RooFit::OwningPtr< RooDataHist > GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual RooFit::OwningPtr< RooDataSet > GetAsDataSet(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
virtual THnSparse * GetAsSparseHist(RooAbsCollection *whichVars=nullptr) const
Get a clone of the markov chain on which this interval is based as a sparse histogram.
virtual RooRealVar * GetNLLVar() const
get a clone of the NLL variable
virtual const RooDataSet * GetAsConstDataSet() const
Definition MarkovChain.h:77
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
Efficient multidimensional histogram.
Definition THnSparse.h:36
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35
Namespace for the RooStats classes.
Definition Asimov.h:19