Logo ROOT   6.21/01
Reference Guide
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 "TObject.h"
18 #include "RooAbsPdf.h"
19 #include "RooAbsData.h"
20 #include "RooArgSet.h"
21 #include "RooArgList.h"
24 #include "RooStats/MCMCInterval.h"
25 
26 
27 namespace 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.
42  MCMCCalculator(RooAbsData& data, const ModelConfig& model);
43 
44  virtual ~MCMCCalculator() {}
45 
46  /// Main interface to get a ConfInterval
47  virtual MCMCInterval* GetInterval() const;
48 
49  /// Get the size of the test (eg. rate of Type I error)
50  virtual Double_t Size() const {return fSize;}
51  /// Get the Confidence level for the test
52  virtual Double_t ConfidenceLevel() const {return 1.-fSize;}
53 
54  virtual void SetModel(const ModelConfig & model);
55 
56  /// Set the DataSet if not already there
57  virtual void SetData(RooAbsData& data) { 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
81  virtual void SetGlobalObservables(const RooArgSet& set) {fGlobalObs.removeAll(); fGlobalObs.add(set);}
82 
83  /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
84  virtual void SetTestSize(Double_t size) {fSize = size;}
85 
86  /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
87  virtual void SetConfidenceLevel(Double_t cl) {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_t useKeys) { fUseKeys = useKeys; }
109  /// set whether to use sparse histogram (if using histogram at all)
110  virtual void SetUseSparseHist(Bool_t 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_t 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
139  fEpsilon = epsilon;
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)
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_t 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_t fUseKeys; // whether to use kernel estimation to determine interval
183  Bool_t fUseSparseHist; // whether to use sparse histogram (if using hist at all)
184  Double_t fLeftSideTF; // left side tail-fraction for interval
185  Double_t fEpsilon; // acceptable error for Keys interval determination
186 
187  Double_t 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  TIterator* it = coll.createIterator();
199  RooAbsArg* r;
200  while ((r = (RooAbsArg*)it->Next()) != NULL)
201  if (dynamic_cast<RooRealVar*>(r))
202  ((RooRealVar*)r)->setBins(numBins);
203  delete it;
204  }
205 
206  ClassDef(MCMCCalculator,4) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
207  };
208 }
209 
210 
211 #endif
MCMCCalculator()
default constructor
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
#define coutE(a)
Definition: RooMsgService.h:34
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual Double_t Size() const
Get the size of the test (eg. rate of Type I error)
virtual void SetGlobalObservables(const RooArgSet &set)
set the global observables which will be used when creating the NLL so the constraint pdf&#39;s will be n...
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
IntervalCalculator is an interface class for a tools which produce RooStats ConfIntervals.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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 SetParameters(const RooArgSet &set)
specify the parameters of interest in the interval
void SetupBasicUsage()
Constructor for automatic configuration with basic settings.
virtual void SetNumBins(Int_t numBins)
set the number of bins to create for each axis when constructing the interval
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual void SetKeysTerminationThreshold(Double_t delta)
When the shortest interval using Keys PDF could not be found to have the desired confidence level +/-...
Iterator abstract base class.
Definition: TIterator.h:30
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval) ...
#define ClassDef(name, id)
Definition: Rtypes.h:326
virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
set what type of interval to have the MCMCInterval represent
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual void SetConditionalObservables(const RooArgSet &set)
set the conditional observables which will be used when creating the NLL so the pdf&#39;s will not be nor...
virtual void SetPdf(RooAbsPdf &pdf)
Set the Pdf if not already there.
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
virtual void SetNumIters(Int_t numIters)
set the number of iterations to run the metropolis algorithm
virtual void SetTestSize(Double_t size)
set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval) ...
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual void SetAxes(RooArgList &axes)
set which variables to put on each axis
auto * a
Definition: textangle.C:12
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
REAL epsilon
Definition: triangle.c:617
virtual void SetData(RooAbsData &data)
Set the DataSet if not already there.
void SetBins(const RooAbsCollection &coll, Int_t numBins) const
virtual void SetKeysConfidenceAccuracy(Double_t epsilon)
Set the desired level of confidence-level accuracy for Keys interval determination.
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
Namespace for the RooStats classes.
Definition: Asimov.h:20
virtual void SetUseKeys(Bool_t useKeys)
set whether to use kernel estimation to determine the interval
virtual void SetProposalFunction(ProposalFunction &proposalFunction)
set the proposal function for suggesting new points for the MCMC
double Double_t
Definition: RtypesCore.h:55
ProposalFunction * fPropFunc
enum MCMCInterval::IntervalType fIntervalType
virtual void SetNuisanceParameters(const RooArgSet &set)
specify the nuisance parameters (eg. the rest of the parameters)
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects...
virtual void SetChainParameters(const RooArgSet &set)
specify the parameters to store in the Markov chain By default all the parameters are stored ...
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual MCMCInterval * GetInterval() const
Main interface to get a ConfInterval.
virtual TObject * Next()=0
virtual void SetLeftSideTailFraction(Double_t a)
Set the left side tail fraction.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
virtual void SetUseSparseHist(Bool_t useSparseHist)
set whether to use sparse histogram (if using histogram at all)
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
virtual void SetModel(const ModelConfig &model)
Set the Model.
virtual void SetPriorPdf(RooAbsPdf &pdf)
Set the Prior Pdf if not already there.