Logo ROOT   6.10/09
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 
17 #include "RooStats/ConfInterval.h"
18 #include "RooArgSet.h"
19 #include "RooArgList.h"
20 #include "RooStats/MarkovChain.h"
21 
22 class RooNDKeysPdf;
23 class RooProduct;
24 
25 
26 namespace 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
99  virtual RooArgList* GetAxes()
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
128  virtual Double_t LowerLimitByDataHist(RooRealVar& param);
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
151  virtual Double_t UpperLimitByDataHist(RooRealVar& param);
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*
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*
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*
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
232  fEpsilon = epsilon;
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
virtual Bool_t IsInInterval(const RooArgSet &point) const
determine whether this point is in the confidence interval
MarkovChain * fChain
Definition: MCMCInterval.h:276
virtual const RooDataSet * GetChainAsConstDataSet()
Get the markov chain on which this interval is based as a RooDataSet.
Definition: MCMCInterval.h:204
#define coutE(a)
Definition: RooMsgService.h:34
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 RooRealVar * GetNLLVar() const
Get a clone of the NLL variable from the markov chain.
Definition: MCMCInterval.h:218
virtual Double_t UpperLimitBySparseHist(RooRealVar &param)
determine upper limit using histogram
virtual Double_t 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); ...
virtual RooAbsArg * addClone(const RooAbsArg &var, Bool_t silent=kFALSE)
Add a clone of the specified argument to list.
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 RooArgSet * GetParameters() const
return a set containing the parameters of this interval the caller owns the returned RooArgSet* ...
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
virtual void CreateKeysPdf()
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 void DetermineByKeys()
virtual void CreateVector(RooRealVar *param)
virtual Double_t GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
std::vector< Int_t > fVector
Definition: MCMCInterval.h:295
RooRealVar * fCutoffVar
Definition: MCMCInterval.h:288
virtual Bool_t GetUseKeys()
get whether we used kernel estimation to determine the interval
Definition: MCMCInterval.h:170
virtual Double_t UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
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 RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList ...
Definition: MCMCInterval.h:99
virtual void DetermineInterval()
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
#define NULL
Definition: RtypesCore.h:88
virtual RooRealVar * GetWeightVar() const
Get a clone of the weight variable from the markov chain.
Definition: MCMCInterval.h:222
virtual void CreateKeysDataHist()
Float_t delta
virtual void DetermineShortestInterval()
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual void CreateHist()
virtual Double_t LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual void SetEpsilon(Double_t epsilon)
set the acceptable level or error for Keys interval determination
Definition: MCMCInterval.h:226
RooNDKeysPdf * fKeysPdf
Definition: MCMCInterval.h:284
Bool_t WithinDeltaFraction(Double_t a, Double_t b)
#define ClassDef(name, id)
Definition: Rtypes.h:297
Efficient multidimensional histogram.
Definition: THnSparse.h:36
virtual Double_t UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual void SetShortestInterval()
Definition: MCMCInterval.h:241
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 Int_t GetDimension() const
Get the number of parameters of interest in this interval.
Definition: MCMCInterval.h:191
virtual void SetChain(MarkovChain &chain)
Set the MarkovChain that this interval is based on.
Definition: MCMCInterval.h:90
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
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
RooDataHist * fKeysDataHist
Definition: MCMCInterval.h:287
virtual void DetermineByDataHist()
MCMCInterval(const char *name=0)
default constructor
enum IntervalType fIntervalType
Definition: MCMCInterval.h:323
virtual Double_t GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval ...
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual Double_t LowerLimitBySparseHist(RooRealVar &param)
determine lower limit using histogram
virtual const RooDataSet * GetAsConstDataSet() const
Definition: MarkovChain.h:77
virtual void SetIntervalType(enum IntervalType intervalType)
Set the type of interval to find.
Definition: MCMCInterval.h:239
virtual void DetermineByHist()
virtual void SetLeftSideTailFraction(Double_t a)
set the left-side tail fraction for a tail-fraction interval
Definition: MCMCInterval.h:247
RooDataHist * fDataHist
Definition: MCMCInterval.h:279
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
virtual void SetConfidenceLevel(Double_t cl)
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Represents the Heaviside function.
Definition: Heaviside.h:18
Generic N-dimensional implementation of a kernel estimation p.d.f.
Definition: RooNDKeysPdf.h:45
REAL epsilon
Definition: triangle.c:617
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
virtual Double_t GetActualConfidenceLevel()
virtual Double_t GetKeysPdfCutoff() { return fKeysCutoff; }
virtual void SetUseKeys(Bool_t useKeys)
set whether to use kernel estimation to determine the interval
Definition: MCMCInterval.h:162
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
Namespace for the RooStats classes.
Definition: Asimov.h:20
double Double_t
Definition: RtypesCore.h:55
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 void CreateDataHist()
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
The TH1 histogram class.
Definition: TH1.h:56
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
Bool_t AcceptableConfLevel(Double_t confLevel)
virtual RooRealVar * GetWeightVar() const
get a clone of the weight variable
Definition: MarkovChain.h:104
Double_t GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
virtual void CreateSparseHist()
virtual void DetermineTailFractionInterval()
virtual Double_t ConfidenceLevel() const
get the desired confidence level (see GetActualConfidenceLevel())
Definition: MCMCInterval.h:59
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
Definition: MCMCInterval.h:244
virtual RooRealVar * GetNLLVar() const
get a clone of the NLL variable
Definition: MarkovChain.h:100
virtual Double_t UpperLimitByDataHist(RooRealVar &param)
determine upper limit using histogram
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 SetParameters(const RooArgSet &parameters)
Set the parameters of interest for this interval and change other internal data members accordingly...
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Double_t UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
virtual Double_t LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void DetermineBySparseHist()
virtual Double_t UpperLimitShortest(RooRealVar &param)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
virtual Double_t CalcConfLevel(Double_t cutoff, Double_t full)
virtual Int_t GetNumBurnInSteps()
get the number of steps in the chain to discard as burn-in,
Definition: MCMCInterval.h:176
Bool_t CheckParameters(const RooArgSet &point) const
check if parameters are correct. (dummy implementation to start)
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 Double_t LowerLimitByDataHist(RooRealVar &param)
determine lower limit using histogram
virtual void SetUseSparseHist(Bool_t useSparseHist)
set whether to use a sparse histogram.
Definition: MCMCInterval.h:166