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