root/RooStats/MetropolisHastings.h

0001 // @(#)root/roostats:$Id$
0002 // Authors: Kevin Belasco        17/06/2009
0003 // Authors: Kyle Cranmer         17/06/2009
0004 /*************************************************************************
0005  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOSTATS_MetropolisHastings
0013 #define ROOSTATS_MetropolisHastings
0014
0015 #include "RooStats/RooStatsUtils.h"
0016 #include "Rtypes.h"
0017 #include "TObject.h"
0018 #include "RooArgSet.h"
0019 #include "RooStats/ProposalFunction.h"
0020 #include "RooStats/MarkovChain.h"
0021
0022 namespace RooStats {
0023
0024    class MetropolisHastings :  public TObject {
0025
0026
0027    public:
0028       enum FunctionSign {kNegative, kPositive, kSignUnset};
0029       enum FunctionType {kRegular, kLog, kTypeUnset};
0030
0031       /// default constructor
0032       MetropolisHastings() = default;
0033
0034       /// alternate constructor
0035       MetropolisHastings(RooAbsReal& function, const RooArgSet& paramsOfInterest,
0036             ProposalFunction& proposalFunction, Int_t numIters);
0037
0038       /// main purpose of MetropolisHastings - run Metropolis-Hastings
0039       /// algorithm to generate Markov Chain of points in the parameter space
0040       virtual MarkovChain* ConstructChain();
0041
0042       /// specify the parameters to store in the chain
0043       /// if not specified all of them will be stored
0044       virtual void SetChainParameters(const RooArgSet& set)
0045       { fChainParams.removeAll();  fChainParams.add(set);  RemoveConstantParameters(&fChainParams); }
0046       /// specify all the parameters of interest in the interval
0047       virtual void SetParameters(const RooArgSet& set)
0048       { fParameters.removeAll();  fParameters.add(set);  RemoveConstantParameters(&fParameters); }
0049       /// set the proposal function for suggesting new points for the MCMC
0050       virtual void SetProposalFunction(ProposalFunction& proposalFunction)
0051       { fPropFunc = &proposalFunction; }
0052       /// set the number of iterations to run the metropolis algorithm
0053       virtual void SetNumIters(Int_t numIters)
0054       { fNumIters = numIters; }
0055       /// set the number of steps in the chain to discard as burn-in,
0056       /// starting from the first
0057       virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
0058       { fNumBurnInSteps = numBurnInSteps; }
0059       /// set the (likelihood) function
0060       virtual void SetFunction(RooAbsReal& function) { fFunction = &function; }
0061       /// set the sign of the function
0062       virtual void SetSign(enum FunctionSign sign) { fSign = sign; }
0063       /// set the type of the function
0064       virtual void SetType(enum FunctionType type) { fType = type; }
0065
0066
0067    protected:
0068       RooAbsReal *fFunction = nullptr; ///< function that will generate likelihood values
0069       RooArgSet fParameters;       ///< RooRealVars that define all parameter space
0070       RooArgSet fChainParams;      ///< RooRealVars that are stored in the chain
0071       ProposalFunction *fPropFunc = nullptr; ///< Proposal function for MCMC integration
0072       Int_t fNumIters = 0;                   ///< number of iterations to run metropolis algorithm
0073       Int_t fNumBurnInSteps = 0;             ///< number of iterations to discard as burn-in, starting from the first
0074       enum FunctionSign fSign = kSignUnset;  ///< whether the likelihood is negative (like NLL) or positive
0075       enum FunctionType fType = kTypeUnset;  ///< whether the likelihood is on a regular, log, (or other) scale
0076
0077       // whether we should take the step, based on the value of d, fSign, fType
0078       virtual bool ShouldTakeStep(double d);
0079       virtual double CalcNLL(double xL);
0080
0081       ClassDefOverride(MetropolisHastings,2) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
0082    };
0083 }
0084
0085
0086 #endif