Logo ROOT   6.10/09
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.
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 size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
80  virtual void SetTestSize(Double_t size) {fSize = size;}
81 
82  /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
83  virtual void SetConfidenceLevel(Double_t cl) {fSize = 1.-cl;}
84 
85  /// set the proposal function for suggesting new points for the MCMC
86  virtual void SetProposalFunction(ProposalFunction& proposalFunction)
87  { fPropFunc = &proposalFunction; }
88 
89  /// set the number of iterations to run the metropolis algorithm
90  virtual void SetNumIters(Int_t numIters)
91  { fNumIters = numIters; }
92 
93  /// set the number of steps in the chain to discard as burn-in,
94  /// starting from the first
95  virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
96  { fNumBurnInSteps = numBurnInSteps; }
97 
98  /// set the number of bins to create for each axis when constructing the interval
99  virtual void SetNumBins(Int_t numBins) { fNumBins = numBins; }
100  /// set which variables to put on each axis
101  virtual void SetAxes(RooArgList& axes)
102  { fAxes = &axes; }
103  /// set whether to use kernel estimation to determine the interval
104  virtual void SetUseKeys(Bool_t useKeys) { fUseKeys = useKeys; }
105  /// set whether to use sparse histogram (if using histogram at all)
106  virtual void SetUseSparseHist(Bool_t useSparseHist)
107  { fUseSparseHist = useSparseHist; }
108 
109  /// set what type of interval to have the MCMCInterval represent
110  virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
111  { fIntervalType = intervalType; }
112 
113  /// Set the left side tail fraction. This will automatically configure the
114  /// MCMCInterval to find a tail-fraction interval.
115  /// Note: that `a' must be in the range 0 <= a <= 1
116  /// or the user will be notified of the error
117  virtual void SetLeftSideTailFraction(Double_t a);
118 
119  /// Set the desired level of confidence-level accuracy for Keys interval
120  /// determination.
121  //
122  /// When determining the cutoff PDF height that gives the
123  /// desired confidence level (C_d), the algorithm will consider acceptable
124  /// any found confidence level c such that Abs(c - C_d) < epsilon.
125  ///
126  /// Any value of this "epsilon" > 0 is considered acceptable, though it is
127  /// advisable to not use a value too small, because the integration of the
128  /// Keys PDF sometimes does not have extremely high accuracy.
130  {
131  if (epsilon < 0)
132  coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
133  << "negative epsilon value" << std::endl;
134  else
135  fEpsilon = epsilon;
136  }
137 
138  /// When the shortest interval using Keys PDF could not be found to have
139  /// the desired confidence level +/- the accuracy (see
140  /// SetKeysConfidenceAccuracy()), the interval determination algorithm
141  /// will have to terminate with an unsatisfactory confidence level when
142  /// the bottom and top of the cutoff search range are very close to being
143  /// equal. This scenario comes into play when there seems to be an error
144  /// in the accuracy of the Keys PDF integration, so the search range
145  /// continues to shrink without converging to a cutoff value that will
146  /// give an acceptable confidence level. To choose how small to allow the
147  /// search range to be before terminating, set the fraction delta such
148  /// that the search will terminate when topCutoff (a) and bottomCutoff (b)
149  /// satisfy this condition:
150  ///
151  /// TMath::Abs(a - b) < TMath::Abs(delta * (a + b)/2)
153  {
154  if (delta < 0.)
155  coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
156  << "negative delta value" << std::endl;
157  else
158  fDelta = delta;
159  }
160 
161  protected:
162 
163  Double_t fSize; // size of the test (eg. specified rate of Type I error)
164  RooArgSet fPOI; // parameters of interest for interval
165  RooArgSet fNuisParams; // nuisance parameters for interval (not really used)
166  RooArgSet fChainParams; // parameters to store in the chain (if not specified they are all of them )
167  RooArgSet fConditionalObs; // conditional observables
168  mutable ProposalFunction* fPropFunc; // Proposal function for MCMC integration
169  RooAbsPdf * fPdf; // pointer to common PDF (owned by the workspace)
170  RooAbsPdf * fPriorPdf; // pointer to prior PDF (owned by the workspace)
171  RooAbsData * fData; // pointer to the data (owned by the workspace)
172  Int_t fNumIters; // number of iterations to run metropolis algorithm
173  Int_t fNumBurnInSteps; // number of iterations to discard as burn-in, starting from the first
174  Int_t fNumBins; // set the number of bins to create for each
175  // axis when constructing the interval
176  RooArgList * fAxes; // which variables to put on each axis
177  Bool_t fUseKeys; // whether to use kernel estimation to determine interval
178  Bool_t fUseSparseHist; // whether to use sparse histogram (if using hist at all)
179  Double_t fLeftSideTF; // left side tail-fraction for interval
180  Double_t fEpsilon; // acceptable error for Keys interval determination
181 
182  Double_t fDelta; // acceptable error for Keys cutoffs being equal
183  // topCutoff (a) considered == bottomCutoff (b) iff
184  // (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2));
185  // Theoretically, the Abs is not needed here, but
186  // floating-point arithmetic does not always work
187  // perfectly, and the Abs doesn't hurt
188  enum MCMCInterval::IntervalType fIntervalType; // type of interval to find
189 
190  void SetupBasicUsage();
191  void SetBins(const RooAbsCollection& coll, Int_t numBins) const
192  {
193  TIterator* it = coll.createIterator();
194  RooAbsArg* r;
195  while ((r = (RooAbsArg*)it->Next()) != NULL)
196  if (dynamic_cast<RooRealVar*>(r))
197  ((RooRealVar*)r)->setBins(numBins);
198  delete it;
199  }
200 
201  ClassDef(MCMCCalculator,3) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
202  };
203 }
204 
205 
206 #endif
MCMCCalculator()
default constructor
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
TIterator * createIterator(Bool_t dir=kIterForward) const
#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 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:86
IntervalCalculator is an interface class for a tools which produce RooStats ConfIntervals.
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
TArc * a
Definition: textangle.C:12
#define NULL
Definition: RtypesCore.h:88
virtual void SetKeysTerminationThreshold(Double_t delta)
When the shortest interval using Keys PDF could not be found to have the desired confidence level +/-...
Float_t delta
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:297
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.
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) ...
TRandom2 r(17)
virtual void SetAxes(RooArgList &axes)
set which variables to put on each axis
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
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 is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
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.
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:66
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.