// @(#)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

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

#ifndef ROOT_TObject
#include "TObject.h"
#endif
#ifndef ROO_ABS_PDF
#include "RooAbsPdf.h"
#endif
#ifndef ROO_ABS_DATA
#include "RooAbsData.h"
#endif
#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif
#ifndef ROO_ARG_LIST
#include "RooArgList.h"
#endif
#ifndef ROOSTATS_ProposalFunction
#include "RooStats/ProposalFunction.h"
#endif
#ifndef ROOSTATS_IntervalCalculator
#include "RooStats/IntervalCalculator.h"
#endif
#ifndef RooStats_MCMCInterval
#include "RooStats/MCMCInterval.h"
#endif

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);

      virtual ~MCMCCalculator() {}

      // 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;}

      virtual void SetModel(const ModelConfig & model); 

      // Set the DataSet if not already there
      virtual void SetData(RooAbsData& data) { 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 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; }
      // set whether to use kernel estimation to determine the interval
      virtual void SetUseKeys(Bool_t useKeys) { fUseKeys = useKeys; }
      // set whether to use sparse histogram (if using histogram at all)
      virtual void SetUseSparseHist(Bool_t 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_t 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_t 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:
      //
      // TMath::Abs(a - b) < TMath::Abs(delta * (a + b)/2)
      virtual void SetKeysTerminationThreshold(Double_t delta)
      {
         if (delta < 0.)
            coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
                                  << "negative delta value" << std::endl;
         else
            fDelta = delta;
      }

   protected:

      Double_t fSize;   // 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
      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; // 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
      Bool_t fUseKeys; // whether to use kernel estimation to determine interval
      Bool_t fUseSparseHist; // whether to use sparse histogram (if using hist at all)
      Double_t fLeftSideTF; // left side tail-fraction for interval
      Double_t fEpsilon; // acceptable error for Keys interval determination

      Double_t fDelta; // acceptable error for Keys cutoffs being equal
                       // topCutoff (a) considered == bottomCutoff (b) iff
                       // (TMath::Abs(a - b) < TMath::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; // type of interval to find

      void SetupBasicUsage();
      void SetBins(const RooAbsCollection& coll, Int_t numBins) const
      {
         TIterator* it = coll.createIterator();
         RooAbsArg* r;
         while ((r = (RooAbsArg*)it->Next()) != NULL)
            if (dynamic_cast<RooRealVar*>(r))
               ((RooRealVar*)r)->setBins(numBins);
         delete it;
      }

      ClassDef(MCMCCalculator,3) // 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
 MCMCCalculator.h:152
 MCMCCalculator.h:153
 MCMCCalculator.h:154
 MCMCCalculator.h:155
 MCMCCalculator.h:156
 MCMCCalculator.h:157
 MCMCCalculator.h:158
 MCMCCalculator.h:159
 MCMCCalculator.h:160
 MCMCCalculator.h:161
 MCMCCalculator.h:162
 MCMCCalculator.h:163
 MCMCCalculator.h:164
 MCMCCalculator.h:165
 MCMCCalculator.h:166
 MCMCCalculator.h:167
 MCMCCalculator.h:168
 MCMCCalculator.h:169
 MCMCCalculator.h:170
 MCMCCalculator.h:171
 MCMCCalculator.h:172
 MCMCCalculator.h:173
 MCMCCalculator.h:174
 MCMCCalculator.h:175
 MCMCCalculator.h:176
 MCMCCalculator.h:177
 MCMCCalculator.h:178
 MCMCCalculator.h:179
 MCMCCalculator.h:180
 MCMCCalculator.h:181
 MCMCCalculator.h:182
 MCMCCalculator.h:183
 MCMCCalculator.h:184
 MCMCCalculator.h:185
 MCMCCalculator.h:186
 MCMCCalculator.h:187
 MCMCCalculator.h:188
 MCMCCalculator.h:189
 MCMCCalculator.h:190
 MCMCCalculator.h:191
 MCMCCalculator.h:192
 MCMCCalculator.h:193
 MCMCCalculator.h:194
 MCMCCalculator.h:195
 MCMCCalculator.h:196
 MCMCCalculator.h:197
 MCMCCalculator.h:198
 MCMCCalculator.h:199
 MCMCCalculator.h:200
 MCMCCalculator.h:201
 MCMCCalculator.h:202
 MCMCCalculator.h:203
 MCMCCalculator.h:204
 MCMCCalculator.h:205
 MCMCCalculator.h:206
 MCMCCalculator.h:207
 MCMCCalculator.h:208
 MCMCCalculator.h:209
 MCMCCalculator.h:210
 MCMCCalculator.h:211
 MCMCCalculator.h:212
 MCMCCalculator.h:213
 MCMCCalculator.h:214
 MCMCCalculator.h:215
 MCMCCalculator.h:216
 MCMCCalculator.h:217
 MCMCCalculator.h:218
 MCMCCalculator.h:219
 MCMCCalculator.h:220
 MCMCCalculator.h:221
 MCMCCalculator.h:222
 MCMCCalculator.h:223