// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
 * 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_BayesianCalculator
#define ROOSTATS_BayesianCalculator

#include "TNamed.h"

#include "Math/IFunctionfwd.h"

#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif

#ifndef ROOSTATS_IntervalCalculator
#include "RooStats/IntervalCalculator.h"
#endif

#ifndef ROOSTATS_SimpleInterval
#include "RooStats/SimpleInterval.h"
#endif

class RooAbsData; 
class RooAbsPdf; 
class RooPlot; 
class RooAbsReal;
class TF1;

namespace RooStats {

   class ModelConfig; 
   class SimpleInterval; 

   class BayesianCalculator : public IntervalCalculator, public TNamed {

   public:

      // constructor
      BayesianCalculator( );

      BayesianCalculator( RooAbsData& data,
                          RooAbsPdf& pdf,
                          const RooArgSet& POI,
                          RooAbsPdf& priorPdf,
                          const RooArgSet* nuisanceParameters = 0 );

      BayesianCalculator( RooAbsData& data,
                          ModelConfig& model );

      // destructor
      virtual ~BayesianCalculator();

      // get the plot with option to get it normalized 
      RooPlot* GetPosteriorPlot(bool norm = false, double precision = 0.01) const; 

      // return posterior pdf (object is managed by the BayesianCalculator class)
      RooAbsPdf* GetPosteriorPdf() const; 
      // return posterior function (object is managed by the BayesianCalculator class)
      RooAbsReal* GetPosteriorFunction() const; 

      // compute the interval. By Default a central interval is computed 
      // By using SetLeftTileFraction can control if central/ upper/lower interval
      // For shortest interval use SetShortestInterval(true)
      virtual SimpleInterval* GetInterval() const ; 

      virtual void SetData( RooAbsData & data ) {
         fData = &data;
         ClearAll();
      }


      // set the model via the ModelConfig
      virtual void SetModel( const ModelConfig& model ); 

      // specify the parameters of interest in the interval
      virtual void SetParameters(const RooArgSet& set) { fPOI.removeAll(); fPOI.add(set); }

      // specify the nuisance parameters (eg. the rest of the parameters)
      virtual void SetNuisanceParameters(const RooArgSet& set) {fNuisanceParameters.removeAll(); fNuisanceParameters.add(set);}

      // Set only the Prior Pdf 
      virtual void SetPriorPdf(RooAbsPdf& pdf) { fPriorPdf = &pdf; }

      // 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;
         fValidInterval = false; 
      }
      // set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
      virtual void SetConfidenceLevel( Double_t cl ) { SetTestSize(1.-cl); }
      // 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 the fraction of probability content on the left tail
      // Central limits use 0.5 (default case)  
      // for upper limits it is 0 and 1 for lower limit
      // For shortest intervals a negative value (i.e. -1) must be given
      void SetLeftSideTailFraction(Double_t leftSideFraction )  {fLeftSideFraction = leftSideFraction;} 

      // set the Bayesian calculator to compute the shorest interval (default is central interval) 
      // to switch off SetLeftSideTailFraction to the rght value
      void SetShortestInterval() { fLeftSideFraction = -1; }

      // set the precision of the Root Finder 
      void SetBrfPrecision( double precision ) { fBrfPrecision = precision; }

      // use directly the approximate posterior function obtained by binning it in nbins
      // by default the cdf is used by integrating the posterior
      // if a value of nbin <= 0 the cdf function will be used
      void SetScanOfPosterior(int nbin = 100) { fNScanBins = nbin; }

      // set the number of iterations when running a MC integration algorithm
      // If not set use default algorithmic values 
      // In case of ToyMC sampling of the nuisance the value is 100
      // In case of using the GSL MCintegrations types the default value is 
      // defined in ROOT::Math::IntegratorMultiDimOptions::DefaultNCalls()
      virtual void SetNumIters(Int_t numIters)  { fNumIterations = numIters; }

      // set the integration type (possible type are) : 
      void SetIntegrationType(const char * type); 

      // return the mode (most probable value of the posterior function) 
      double GetMode() const; 

      // force the nuisance pdf when using the toy mc sampling
      void ForceNuisancePdf(RooAbsPdf & pdf) { fNuisancePdf = &pdf; }

   protected:

      void ClearAll() const; 

      void ApproximatePosterior() const; 

      void ComputeIntervalFromApproxPosterior(double c1, double c2) const;

      void ComputeIntervalFromCdf(double c1, double c2) const; 

      void ComputeIntervalUsingRooFit(double c1, double c2) const;

      void ComputeShortestInterval() const; 
   
   private:

      // plan to replace the above: return a SimpleInterval integrating 
      // over all other parameters except the one specified as argument
      //virtual SimpleInterval* GetInterval( RooRealVar* parameter  ) const { return 0; }
    
      RooAbsData* fData;                          // data set 
      RooAbsPdf* fPdf;                           // model pdf  (could contain the nuisance pdf as constraint term)
      RooArgSet fPOI;                            // POI
      RooAbsPdf* fPriorPdf;                      // prior pdf (typically for the POI)
      RooAbsPdf* fNuisancePdf;                   // nuisance pdf (needed when using nuisance sampling technique)
      RooArgSet fNuisanceParameters;             // nuisance parameters
      RooArgSet fConditionalObs    ;             // conditional observables

      mutable RooAbsPdf* fProductPdf;              // internal pointer to model * prior
      mutable RooAbsReal* fLogLike;                // internal pointer to log likelihood function
      mutable RooAbsReal* fLikelihood;             // internal pointer to likelihood function 
      mutable RooAbsReal* fIntegratedLikelihood;   // integrated likelihood function, i.e - unnormalized posterior function  
      mutable RooAbsPdf* fPosteriorPdf;             // normalized (on the poi) posterior pdf 
      mutable ROOT::Math::IGenFunction * fPosteriorFunction;   // function representing the posterior
      mutable TF1 * fApproxPosterior;    // TF1 representing the scanned posterior function
      mutable Double_t  fLower;    // computer lower interval bound
      mutable Double_t  fUpper;    // upper interval bound
      mutable Double_t  fNLLMin;   // minimum value of Nll 
      double fSize;  // size used for getting the interval
      double fLeftSideFraction;    // fraction of probability content on left side of interval
      double fBrfPrecision;     // root finder precision
      mutable int fNScanBins;            // number of bins to scan, if = -1 no scan is done (default)
      int fNumIterations;        // number of iterations (when using ToyMC)
      mutable Bool_t    fValidInterval; 
      


      TString fIntegrationType; 

   protected:

      ClassDef(BayesianCalculator,2)  // BayesianCalculator class

   };
}

#endif
 BayesianCalculator.h:1
 BayesianCalculator.h:2
 BayesianCalculator.h:3
 BayesianCalculator.h:4
 BayesianCalculator.h:5
 BayesianCalculator.h:6
 BayesianCalculator.h:7
 BayesianCalculator.h:8
 BayesianCalculator.h:9
 BayesianCalculator.h:10
 BayesianCalculator.h:11
 BayesianCalculator.h:12
 BayesianCalculator.h:13
 BayesianCalculator.h:14
 BayesianCalculator.h:15
 BayesianCalculator.h:16
 BayesianCalculator.h:17
 BayesianCalculator.h:18
 BayesianCalculator.h:19
 BayesianCalculator.h:20
 BayesianCalculator.h:21
 BayesianCalculator.h:22
 BayesianCalculator.h:23
 BayesianCalculator.h:24
 BayesianCalculator.h:25
 BayesianCalculator.h:26
 BayesianCalculator.h:27
 BayesianCalculator.h:28
 BayesianCalculator.h:29
 BayesianCalculator.h:30
 BayesianCalculator.h:31
 BayesianCalculator.h:32
 BayesianCalculator.h:33
 BayesianCalculator.h:34
 BayesianCalculator.h:35
 BayesianCalculator.h:36
 BayesianCalculator.h:37
 BayesianCalculator.h:38
 BayesianCalculator.h:39
 BayesianCalculator.h:40
 BayesianCalculator.h:41
 BayesianCalculator.h:42
 BayesianCalculator.h:43
 BayesianCalculator.h:44
 BayesianCalculator.h:45
 BayesianCalculator.h:46
 BayesianCalculator.h:47
 BayesianCalculator.h:48
 BayesianCalculator.h:49
 BayesianCalculator.h:50
 BayesianCalculator.h:51
 BayesianCalculator.h:52
 BayesianCalculator.h:53
 BayesianCalculator.h:54
 BayesianCalculator.h:55
 BayesianCalculator.h:56
 BayesianCalculator.h:57
 BayesianCalculator.h:58
 BayesianCalculator.h:59
 BayesianCalculator.h:60
 BayesianCalculator.h:61
 BayesianCalculator.h:62
 BayesianCalculator.h:63
 BayesianCalculator.h:64
 BayesianCalculator.h:65
 BayesianCalculator.h:66
 BayesianCalculator.h:67
 BayesianCalculator.h:68
 BayesianCalculator.h:69
 BayesianCalculator.h:70
 BayesianCalculator.h:71
 BayesianCalculator.h:72
 BayesianCalculator.h:73
 BayesianCalculator.h:74
 BayesianCalculator.h:75
 BayesianCalculator.h:76
 BayesianCalculator.h:77
 BayesianCalculator.h:78
 BayesianCalculator.h:79
 BayesianCalculator.h:80
 BayesianCalculator.h:81
 BayesianCalculator.h:82
 BayesianCalculator.h:83
 BayesianCalculator.h:84
 BayesianCalculator.h:85
 BayesianCalculator.h:86
 BayesianCalculator.h:87
 BayesianCalculator.h:88
 BayesianCalculator.h:89
 BayesianCalculator.h:90
 BayesianCalculator.h:91
 BayesianCalculator.h:92
 BayesianCalculator.h:93
 BayesianCalculator.h:94
 BayesianCalculator.h:95
 BayesianCalculator.h:96
 BayesianCalculator.h:97
 BayesianCalculator.h:98
 BayesianCalculator.h:99
 BayesianCalculator.h:100
 BayesianCalculator.h:101
 BayesianCalculator.h:102
 BayesianCalculator.h:103
 BayesianCalculator.h:104
 BayesianCalculator.h:105
 BayesianCalculator.h:106
 BayesianCalculator.h:107
 BayesianCalculator.h:108
 BayesianCalculator.h:109
 BayesianCalculator.h:110
 BayesianCalculator.h:111
 BayesianCalculator.h:112
 BayesianCalculator.h:113
 BayesianCalculator.h:114
 BayesianCalculator.h:115
 BayesianCalculator.h:116
 BayesianCalculator.h:117
 BayesianCalculator.h:118
 BayesianCalculator.h:119
 BayesianCalculator.h:120
 BayesianCalculator.h:121
 BayesianCalculator.h:122
 BayesianCalculator.h:123
 BayesianCalculator.h:124
 BayesianCalculator.h:125
 BayesianCalculator.h:126
 BayesianCalculator.h:127
 BayesianCalculator.h:128
 BayesianCalculator.h:129
 BayesianCalculator.h:130
 BayesianCalculator.h:131
 BayesianCalculator.h:132
 BayesianCalculator.h:133
 BayesianCalculator.h:134
 BayesianCalculator.h:135
 BayesianCalculator.h:136
 BayesianCalculator.h:137
 BayesianCalculator.h:138
 BayesianCalculator.h:139
 BayesianCalculator.h:140
 BayesianCalculator.h:141
 BayesianCalculator.h:142
 BayesianCalculator.h:143
 BayesianCalculator.h:144
 BayesianCalculator.h:145
 BayesianCalculator.h:146
 BayesianCalculator.h:147
 BayesianCalculator.h:148
 BayesianCalculator.h:149
 BayesianCalculator.h:150
 BayesianCalculator.h:151
 BayesianCalculator.h:152
 BayesianCalculator.h:153
 BayesianCalculator.h:154
 BayesianCalculator.h:155
 BayesianCalculator.h:156
 BayesianCalculator.h:157
 BayesianCalculator.h:158
 BayesianCalculator.h:159
 BayesianCalculator.h:160
 BayesianCalculator.h:161
 BayesianCalculator.h:162
 BayesianCalculator.h:163
 BayesianCalculator.h:164
 BayesianCalculator.h:165
 BayesianCalculator.h:166
 BayesianCalculator.h:167
 BayesianCalculator.h:168
 BayesianCalculator.h:169
 BayesianCalculator.h:170
 BayesianCalculator.h:171
 BayesianCalculator.h:172
 BayesianCalculator.h:173
 BayesianCalculator.h:174
 BayesianCalculator.h:175
 BayesianCalculator.h:176
 BayesianCalculator.h:177
 BayesianCalculator.h:178
 BayesianCalculator.h:179
 BayesianCalculator.h:180
 BayesianCalculator.h:181
 BayesianCalculator.h:182
 BayesianCalculator.h:183
 BayesianCalculator.h:184
 BayesianCalculator.h:185
 BayesianCalculator.h:186
 BayesianCalculator.h:187
 BayesianCalculator.h:188
 BayesianCalculator.h:189
 BayesianCalculator.h:190
 BayesianCalculator.h:191
 BayesianCalculator.h:192
 BayesianCalculator.h:193
 BayesianCalculator.h:194
 BayesianCalculator.h:195
 BayesianCalculator.h:196
 BayesianCalculator.h:197