EIC code displayed by LXR master/​include/​root/​RooStats/​MCMCCalculator.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-12-16 10:29:45

 // @(#)root/roostats:$Id$
 // Authors: Kevin Belasco        17/06/2009
 // Authors: Kyle Cranmer         17/06/2009
 /*************************************************************************
  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
  * All rights reserved.                                                  *
  *                                                                       *
  * For the licensing terms see $ROOTSYS/LICENSE.                         *
  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
  *************************************************************************/
 
 #ifndef ROOSTATS_MCMCCalculator
 #define ROOSTATS_MCMCCalculator
 
 #include "Rtypes.h"
 
 #include "TNamed.h"
 #include "RooAbsPdf.h"
 #include "RooAbsData.h"
 #include "RooArgSet.h"
 #include "RooArgList.h"
 #include "RooStats/ProposalFunction.h"
 #include "RooStats/IntervalCalculator.h"
 #include "RooStats/MCMCInterval.h"
 
 
 namespace RooStats {
 
    class ModelConfig;
 
    class MCMCCalculator : public IntervalCalculator, public TNamed {
 
    public:
       /// default constructor
       MCMCCalculator();
 
       /// Constructor for automatic configuration with basic settings and a
       /// ModelConfig.  Uses a UniformProposal, 10,000 iterations, 40 burn in
       /// steps, 50 bins for each RooRealVar, determines interval by histogram,
       /// and finds a 95% confidence interval.  Any of these basic settings can
       /// be overridden by calling one of the Set...() methods.
       MCMCCalculator(RooAbsData& data, const ModelConfig& model);
 
       /// Main interface to get a ConfInterval
       MCMCInterval* GetInterval() const override;
 
       /// Get the size of the test (eg. rate of Type I error)
       double Size() const override {return fSize;}
       /// Get the Confidence level for the test
       double ConfidenceLevel() const override {return 1.-fSize;}
 
       void SetModel(const ModelConfig & model) override;
 
       /// Set the DataSet if not already there
       void SetData(RooAbsData& data) override { fData = &data; }
 
       /// Set the Pdf if not already there
       virtual void SetPdf(RooAbsPdf& pdf) { fPdf = &pdf; }
 
       /// Set the Prior Pdf if not already there
       virtual void SetPriorPdf(RooAbsPdf& pdf) { fPriorPdf = &pdf; }
 
       /// specify the parameters of interest in the interval
       virtual void SetParameters(const RooArgSet& set) { fPOI.removeAll(); fPOI.add(set); }
 
       /// specify the parameters to store in the Markov chain
       /// By default all the parameters are stored
       virtual void SetChainParameters(const RooArgSet & set) { fChainParams.removeAll(); fChainParams.add(set); }
 
       /// specify the nuisance parameters (eg. the rest of the parameters)
       virtual void SetNuisanceParameters(const RooArgSet& set) {fNuisParams.removeAll(); fNuisParams.add(set);}
 
       /// set the conditional observables which will be used when creating the NLL
       /// so the pdf's will not be normalized on the conditional observables when computing the NLL
       virtual void SetConditionalObservables(const RooArgSet& set) {fConditionalObs.removeAll(); fConditionalObs.add(set);}
 
       /// set the global observables which will be used when creating the NLL
       /// so the constraint pdf's will be normalized correctly on the global observables when computing the NLL
       virtual void SetGlobalObservables(const RooArgSet& set) {fGlobalObs.removeAll(); fGlobalObs.add(set);}
 
       /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
       void SetTestSize(double size) override {fSize = size;}
 
       /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
       void SetConfidenceLevel(double cl) override {fSize = 1.-cl;}
 
       /// set the proposal function for suggesting new points for the MCMC
       virtual void SetProposalFunction(ProposalFunction& proposalFunction)
       { fPropFunc = &proposalFunction; }
 
       /// set the number of iterations to run the metropolis algorithm
       virtual void SetNumIters(Int_t numIters)
       { fNumIters = numIters; }
 
       /// set the number of steps in the chain to discard as burn-in,
       /// starting from the first
       virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
       { fNumBurnInSteps = numBurnInSteps; }
 
       /// set the number of bins to create for each axis when constructing the interval
       virtual void SetNumBins(Int_t numBins) { fNumBins = numBins; }
       /// set which variables to put on each axis
       virtual void SetAxes(RooArgList& axes)
       { fAxes = &axes; }
       /// set whether to use kernel estimation to determine the interval
       virtual void SetUseKeys(bool useKeys) { fUseKeys = useKeys; }
       /// set whether to use sparse histogram (if using histogram at all)
       virtual void SetUseSparseHist(bool useSparseHist)
       { fUseSparseHist = useSparseHist; }
 
       /// set what type of interval to have the MCMCInterval represent
       virtual void SetIntervalType(enum MCMCInterval::IntervalType intervalType)
       { fIntervalType = intervalType; }
 
       /// Set the left side tail fraction. This will automatically configure the
       /// MCMCInterval to find a tail-fraction interval.
       /// Note: that `a' must be in the range 0 <= a <= 1
       /// or the user will be notified of the error
       virtual void SetLeftSideTailFraction(double a);
 
       /// Set the desired level of confidence-level accuracy  for Keys interval
       /// determination.
       //
       /// When determining the cutoff PDF height that gives the
       /// desired confidence level (C_d), the algorithm will consider acceptable
       /// any found confidence level c such that Abs(c - C_d) < epsilon.
       ///
       /// Any value of this "epsilon" > 0 is considered acceptable, though it is
       /// advisable to not use a value too small, because the integration of the
       /// Keys PDF sometimes does not have extremely high accuracy.
       virtual void SetKeysConfidenceAccuracy(double epsilon)
       {
          if (epsilon < 0) {
             coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
                                   << "negative epsilon value" << std::endl;
          } else {
             fEpsilon = epsilon;
          }
       }
 
       /// When the shortest interval using Keys PDF could not be found to have
       /// the desired confidence level +/- the accuracy (see
       /// SetKeysConfidenceAccuracy()), the interval determination algorithm
       /// will have to terminate with an unsatisfactory confidence level when
       /// the bottom and top of the cutoff search range are very close to being
       /// equal.  This scenario comes into play when there seems to be an error
       /// in the accuracy of the Keys PDF integration, so the search range
       /// continues to shrink without converging to a cutoff value that will
       /// give an acceptable confidence level.  To choose how small to allow the
       /// search range to be before terminating, set the fraction delta such
       /// that the search will terminate when topCutoff (a) and bottomCutoff (b)
       /// satisfy this condition:
       ///
       /// std::abs(a - b) < std::abs(delta * (a + b)/2)
       virtual void SetKeysTerminationThreshold(double delta)
       {
          if (delta < 0.) {
             coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
                                   << "negative delta value" << std::endl;
          } else {
             fDelta = delta;
          }
       }
 
    protected:
       double fSize = -1;           ///< size of the test (eg. specified rate of Type I error)
       RooArgSet   fPOI;            ///< parameters of interest for interval
       RooArgSet   fNuisParams;     ///< nuisance parameters for interval (not really used)
       RooArgSet   fChainParams;    ///< parameters to store in the chain (if not specified they are all of them )
       RooArgSet   fConditionalObs; ///< conditional observables
       RooArgSet   fGlobalObs;      ///< global observables
       mutable ProposalFunction* fPropFunc; ///< Proposal function for MCMC integration
       RooAbsPdf * fPdf;      ///< pointer to common PDF (owned by the workspace)
       RooAbsPdf * fPriorPdf; ///< pointer to prior  PDF (owned by the workspace)
       RooAbsData * fData;    ///< pointer to the data (owned by the workspace)
       Int_t fNumIters = 0;   ///< number of iterations to run metropolis algorithm
       Int_t fNumBurnInSteps = 0; ///< number of iterations to discard as burn-in, starting from the first
       Int_t fNumBins = 0;        ///< set the number of bins to create for each
                                  ///< axis when constructing the interval
       RooArgList * fAxes;    ///< which variables to put on each axis
       bool fUseKeys = false; ///< whether to use kernel estimation to determine interval
       bool fUseSparseHist = false; ///< whether to use sparse histogram (if using hist at all)
       double fLeftSideTF = -1;     ///< left side tail-fraction for interval
       double fEpsilon = -1;        ///< acceptable error for Keys interval determination
 
       double fDelta = -1; ///< acceptable error for Keys cutoffs being equal
                           ///< topCutoff (a) considered == bottomCutoff (b) iff
                           ///< (std::abs(a - b) < std::abs(fDelta * (a + b)/2));
                           ///< Theoretically, the Abs is not needed here, but
                           ///< floating-point arithmetic does not always work
                           ///< perfectly, and the Abs doesn't hurt
       enum MCMCInterval::IntervalType fIntervalType = MCMCInterval::kShortest; // type of interval to find
 
       void SetupBasicUsage();
       void SetBins(const RooAbsCollection &coll, Int_t numBins) const
       {
          for (auto *r : dynamic_range_cast<RooRealVar *>(coll)){
             if (r) {
                r->setBins(numBins);
             }
          }
       }
 
       ClassDefOverride(MCMCCalculator,4) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
    };
 }
 
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team