ROOT logo
// @(#)root/roostats:$Id: MCMCCalculator.h 26805 2009-06-17 14:31:02Z kbelasco $
// Author: Kevin Belasco        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

#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif

#include "RooAbsPdf.h"
#include "RooAbsData.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooWorkspace.h"
#include "RooStats/ProposalFunction.h"
#include "RooStats/IntervalCalculator.h"
#include "RooStats/MCMCInterval.h"

namespace RooStats {

   class MCMCCalculator : public IntervalCalculator {

   public:
      // default constructor
      MCMCCalculator();

      // alternate constructor
      MCMCCalculator(RooWorkspace& ws, RooAbsData& data, RooAbsPdf& pdf,
         RooArgSet& paramsOfInterest, ProposalFunction& proposalFunction,
         Int_t numIters, RooArgList* axes = NULL, Double_t size = 0.05);

      // alternate constructor
      MCMCCalculator(RooAbsData& data, RooAbsPdf& pdf,
         RooArgSet& paramsOfInterest, ProposalFunction& proposalFunction,
         Int_t numIters, RooArgList* axes = NULL, Double_t size = 0.05);

      virtual ~MCMCCalculator()
      {
         if (fOwnsWorkspace)
            delete fWS;
      }
    
      // Main interface to get a ConfInterval
      virtual MCMCInterval* GetInterval() const;

      // Get the size of the test (eg. rate of Type I error)
      virtual Double_t Size() const {return fSize;}
      // Get the Confidence level for the test
      virtual Double_t ConfidenceLevel() const {return 1.-fSize;}

      // set a workspace that owns all the necessary components for the analysis
      virtual void SetWorkspace(RooWorkspace & ws)
      {
         if (!fWS)
            fWS = &ws;
         else {
	   //RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;
            fWS->merge(ws);
	    //RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
         }
      }

      // set the name of the data set
      virtual void SetData(const char* data) { fDataName = data; }

      // Set the DataSet, add to the the workspace if not already there
      virtual void SetData(RooAbsData& data)
      {
         if (!fWS) {
            fWS = new RooWorkspace();
            fOwnsWorkspace = true; 
         }
         if (! fWS->data(data.GetName()) ) {
	   RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;
            fWS->import(data);
	    RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
         }
         SetData(data.GetName());
      }

      // set the name of the pdf
      virtual void SetPdf(const char* name) { fPdfName = name; }

      // Set the Pdf, add to the the workspace if not already there
      virtual void SetPdf(RooAbsPdf& pdf)
      {
         if (!fWS) 
            fWS = new RooWorkspace();
         if (! fWS->pdf( pdf.GetName() ))
         {
            RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;
            fWS->import(pdf);
            RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
         }
         SetPdf(pdf.GetName());
      }

      // specify the parameters of interest in the interval
      virtual void SetParameters(RooArgSet& set) {fPOI = &set;}
      // specify the nuisance parameters (eg. the rest of the parameters)
      virtual void SetNuisanceParameters(RooArgSet& set) {fNuisParams = &set;}
      // set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
      virtual void SetTestSize(Double_t size) {fSize = size;}
      // set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
      virtual void SetConfidenceLevel(Double_t cl) {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; }

   protected:
      Double_t fSize; // size of the test (eg. specified rate of Type I error)
      RooWorkspace* fWS; // owns all the components used by the calculator
      RooArgSet* fPOI; // parameters of interest for interval
      RooArgSet* fNuisParams; // nuisance parameters for interval
      Bool_t fOwnsWorkspace; // whether we own the workspace
      ProposalFunction* fPropFunc; // Proposal function for MCMC integration
      const char* fPdfName; // name of common PDF in workspace
      const char* fDataName; // name of data set in workspace
      Int_t fNumIters; // number of iterations to run metropolis algorithm
      Int_t fNumBurnInSteps; // number of iterations to discard as burn-in, starting from the first
      Int_t fNumBins; // set the number of bins to create for each
                      // axis when constructing the interval
      RooArgList* fAxes; // which variables to put on each axis

      ClassDef(MCMCCalculator,1) // Markov Chain Monte Carlo calculator for Bayesian credible intervals
   };
}


#endif
 MCMCCalculator.h:1
 MCMCCalculator.h:2
 MCMCCalculator.h:3
 MCMCCalculator.h:4
 MCMCCalculator.h:5
 MCMCCalculator.h:6
 MCMCCalculator.h:7
 MCMCCalculator.h:8
 MCMCCalculator.h:9
 MCMCCalculator.h:10
 MCMCCalculator.h:11
 MCMCCalculator.h:12
 MCMCCalculator.h:13
 MCMCCalculator.h:14
 MCMCCalculator.h:15
 MCMCCalculator.h:16
 MCMCCalculator.h:17
 MCMCCalculator.h:18
 MCMCCalculator.h:19
 MCMCCalculator.h:20
 MCMCCalculator.h:21
 MCMCCalculator.h:22
 MCMCCalculator.h:23
 MCMCCalculator.h:24
 MCMCCalculator.h:25
 MCMCCalculator.h:26
 MCMCCalculator.h:27
 MCMCCalculator.h:28
 MCMCCalculator.h:29
 MCMCCalculator.h:30
 MCMCCalculator.h:31
 MCMCCalculator.h:32
 MCMCCalculator.h:33
 MCMCCalculator.h:34
 MCMCCalculator.h:35
 MCMCCalculator.h:36
 MCMCCalculator.h:37
 MCMCCalculator.h:38
 MCMCCalculator.h:39
 MCMCCalculator.h:40
 MCMCCalculator.h:41
 MCMCCalculator.h:42
 MCMCCalculator.h:43
 MCMCCalculator.h:44
 MCMCCalculator.h:45
 MCMCCalculator.h:46
 MCMCCalculator.h:47
 MCMCCalculator.h:48
 MCMCCalculator.h:49
 MCMCCalculator.h:50
 MCMCCalculator.h:51
 MCMCCalculator.h:52
 MCMCCalculator.h:53
 MCMCCalculator.h:54
 MCMCCalculator.h:55
 MCMCCalculator.h:56
 MCMCCalculator.h:57
 MCMCCalculator.h:58
 MCMCCalculator.h:59
 MCMCCalculator.h:60
 MCMCCalculator.h:61
 MCMCCalculator.h:62
 MCMCCalculator.h:63
 MCMCCalculator.h:64
 MCMCCalculator.h:65
 MCMCCalculator.h:66
 MCMCCalculator.h:67
 MCMCCalculator.h:68
 MCMCCalculator.h:69
 MCMCCalculator.h:70
 MCMCCalculator.h:71
 MCMCCalculator.h:72
 MCMCCalculator.h:73
 MCMCCalculator.h:74
 MCMCCalculator.h:75
 MCMCCalculator.h:76
 MCMCCalculator.h:77
 MCMCCalculator.h:78
 MCMCCalculator.h:79
 MCMCCalculator.h:80
 MCMCCalculator.h:81
 MCMCCalculator.h:82
 MCMCCalculator.h:83
 MCMCCalculator.h:84
 MCMCCalculator.h:85
 MCMCCalculator.h:86
 MCMCCalculator.h:87
 MCMCCalculator.h:88
 MCMCCalculator.h:89
 MCMCCalculator.h:90
 MCMCCalculator.h:91
 MCMCCalculator.h:92
 MCMCCalculator.h:93
 MCMCCalculator.h:94
 MCMCCalculator.h:95
 MCMCCalculator.h:96
 MCMCCalculator.h:97
 MCMCCalculator.h:98
 MCMCCalculator.h:99
 MCMCCalculator.h:100
 MCMCCalculator.h:101
 MCMCCalculator.h:102
 MCMCCalculator.h:103
 MCMCCalculator.h:104
 MCMCCalculator.h:105
 MCMCCalculator.h:106
 MCMCCalculator.h:107
 MCMCCalculator.h:108
 MCMCCalculator.h:109
 MCMCCalculator.h:110
 MCMCCalculator.h:111
 MCMCCalculator.h:112
 MCMCCalculator.h:113
 MCMCCalculator.h:114
 MCMCCalculator.h:115
 MCMCCalculator.h:116
 MCMCCalculator.h:117
 MCMCCalculator.h:118
 MCMCCalculator.h:119
 MCMCCalculator.h:120
 MCMCCalculator.h:121
 MCMCCalculator.h:122
 MCMCCalculator.h:123
 MCMCCalculator.h:124
 MCMCCalculator.h:125
 MCMCCalculator.h:126
 MCMCCalculator.h:127
 MCMCCalculator.h:128
 MCMCCalculator.h:129
 MCMCCalculator.h:130
 MCMCCalculator.h:131
 MCMCCalculator.h:132
 MCMCCalculator.h:133
 MCMCCalculator.h:134
 MCMCCalculator.h:135
 MCMCCalculator.h:136
 MCMCCalculator.h:137
 MCMCCalculator.h:138
 MCMCCalculator.h:139
 MCMCCalculator.h:140
 MCMCCalculator.h:141
 MCMCCalculator.h:142
 MCMCCalculator.h:143
 MCMCCalculator.h:144
 MCMCCalculator.h:145
 MCMCCalculator.h:146
 MCMCCalculator.h:147
 MCMCCalculator.h:148
 MCMCCalculator.h:149
 MCMCCalculator.h:150
 MCMCCalculator.h:151