ROOT logo
// @(#)root/roostats:$Id: BayesianCalculator.cxx 31876 2009-12-14 11:11:27Z brun $
// 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.             *
*************************************************************************/

//_________________________________________________
/**
   BayesianCalculator is a concrete implementation of IntervalCalculator. 
   It computes the posterior probability density functions using the  
   numerical (or analytical integration) for integrating the product of the 
   likelihood and prior functions (Bayes theorem). 
   The class works only for problems with only one parameter of interest,  
   the posterior is a one-dimensional function
   The class computes via  GetInterval() the central Bayesian credible intervals

   Note: when nuisance parameters are present a multi-dimensional integration is 
   needed. In some cases, when the integration must be performed numerically, evaluating the posterior or 
   getting the interval (calling GetInterval) can result in long execution time. 
   In these case using the MCMCCalculator could be more convenient
**/

// include other header files

#include "RooAbsFunc.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgSet.h"
#include "RooBrentRootFinder.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPlot.h"
#include "RooProdPdf.h"

// include header file of this class 
#include "RooStats/BayesianCalculator.h"
#include "RooStats/ModelConfig.h"

#include "TAxis.h"

ClassImp(RooStats::BayesianCalculator)

namespace RooStats { 


BayesianCalculator::BayesianCalculator() :
  fData(0),
  fPdf(0),
  fPriorPOI(0),
  fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0), 
  fInterval(0),
  fSize(0.05)
{
   // default constructor. Need to call the Setter methods afterwards
}

BayesianCalculator::BayesianCalculator( /* const char* name,  const char* title, */						   
						    RooAbsData& data,
                                                    RooAbsPdf& pdf,
						    const RooArgSet& POI,
						    RooAbsPdf& priorPOI,
						    const RooArgSet* nuisanceParameters ) :
   //TNamed( TString(name), TString(title) ),
  fData(&data),
  fPdf(&pdf),
  fPOI(POI),
  fPriorPOI(&priorPOI),
  fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
  fInterval(0),
  fSize(0.05)
{
   // constructor from data set, model pdf, set with the parameter of interest 
   // (must contain only one parameter for the moment) and prior pdf
   // Optionally an additional set of parameters can be specified (nuisance parameters) 
   // which will be integrated (marginalized) when creating the posterior pdf.
   // A default size of 0.05 is used (for 95% CL interval)
   if (nuisanceParameters) fNuisanceParameters.add(*nuisanceParameters); 
}

BayesianCalculator::BayesianCalculator( RooAbsData& data,
                       ModelConfig & model) : 
   fData(&data), 
   fPdf(model.GetPdf()),
   fPriorPOI( model.GetPriorPdf()),
   fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
   fInterval(0),
   fSize(0.05)
{
   // Same constructor but from data and a ModelConfig describing the model pdf and the prior, the parameter
   // of interest and the nuisance parameters
   SetModel(model);
}


BayesianCalculator::~BayesianCalculator()
{
   // destructor cleaning all managed objects
   ClearAll(); 
}

void BayesianCalculator::ClearAll() const { 
   // clear cached pdf objects (posterior pdf, Likelihood, NLL, etc.) 
   if (fProductPdf) delete fProductPdf; 
   if (fLogLike) delete fLogLike; 
   if (fLikelihood) delete fLikelihood; 
   if (fIntegratedLikelihood) delete fIntegratedLikelihood; 
   if (fPosteriorPdf) delete fPosteriorPdf;      
   if (fInterval) delete fInterval; 
   fPosteriorPdf = 0; 
   fProductPdf = 0;
   fLogLike = 0; 
   fLikelihood = 0; 
   fIntegratedLikelihood = 0; 
   fInterval = 0; 
}

void BayesianCalculator::SetModel(const ModelConfig & model) {
   // set the model configuration 
   fPdf = model.GetPdf();
   fPriorPOI =  model.GetPriorPdf(); 
   // assignment operator = does not do a real copy the sets (must use add method) 
   fPOI.removeAll();
   fNuisanceParameters.removeAll();
   if (model.GetParametersOfInterest()) fPOI.add( *(model.GetParametersOfInterest()) );
   if (model.GetNuisanceParameters())  fNuisanceParameters.add( *(model.GetNuisanceParameters() ) );

   // invalidate the cached pointers
   ClearAll(); 
}

   RooArgSet* BayesianCalculator::GetMode(RooArgSet* /* parameters */) const
{
   // return the mode (not yet implemented) but can be easly obtained from 
   //  GetPosteriorPdf()->asTF(poi)->GetMaximumX();
   return 0;
}

RooAbsPdf* BayesianCalculator::GetPosteriorPdf() const
{
   // get the posterior pdf as a RooAbsPdf 
   // the posterior is obtained from the product of the likelihood function and the 
   // prior pdf which is then intergated in the nuisance parameters (if existing). 
   // A prior function for the nuisance can be specified either in the prior pdf object 
   // or in the model itself. If no prior nuisance is specified, but prior parameters are then 
   // the integration is performed assuming a flat prior for the nuisance parameters.

   // run some checks
   if (!fPdf ) return 0; 
   if (!fPriorPOI) { 
      std::cerr << "BayesianCalculator::GetPosteriorPdf - missing prior pdf" << std::endl;
   }
   if (fPOI.getSize() == 0) return 0; 
   if (fPOI.getSize() > 1) { 
      std::cerr << "BayesianCalculator::GetPosteriorPdf - current implementation works only on 1D intervals" << std::endl;
      return 0; 
   }


   // create a unique name for the product pdf 
   TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPOI->GetName() );   
   fProductPdf = new RooProdPdf(prodName,"",RooArgList(*fPdf,*fPriorPOI));

   RooArgSet* constrainedParams = fProductPdf->getParameters(*fData);

   // use RooFit::Constrain() to make product of likelihood with prior pdf
   fLogLike = fProductPdf->createNLL(*fData, RooFit::Constrain(*constrainedParams) );

   TString likeName = TString("likelihood_") + TString(fProductPdf->GetName());   
   fLikelihood = new RooFormulaVar(likeName,"exp(-@0)",RooArgList(*fLogLike));
   RooAbsReal * plike = fLikelihood; 
   if (fNuisanceParameters.getSize() > 0) { 
      fIntegratedLikelihood = fLikelihood->createIntegral(fNuisanceParameters);
      plike = fIntegratedLikelihood; 
   }

   // create a unique name on the posterior from the names of the components
   TString posteriorName = this->GetName() + TString("_posteriorPdf_") + plike->GetName(); 
   fPosteriorPdf = new RooGenericPdf(posteriorName,"@0",*plike);

   delete constrainedParams;

   return fPosteriorPdf;
}


RooPlot* BayesianCalculator::GetPosteriorPlot() const
{
  /// return a RooPlot with the posterior PDF and the credibility region

  if (!fPosteriorPdf) GetPosteriorPdf();
  if (!fInterval) GetInterval();

  RooAbsRealLValue* poi = dynamic_cast<RooAbsRealLValue*>( fPOI.first() );
  assert(poi);

   RooPlot* plot = poi->frame();

   plot->SetTitle(TString("Posterior probability of parameter \"")+TString(poi->GetName())+TString("\""));  
   fPosteriorPdf->plotOn(plot,RooFit::Range(fInterval->LowerLimit(),fInterval->UpperLimit(),kFALSE),RooFit::VLines(),RooFit::DrawOption("F"),RooFit::MoveToBack(),RooFit::FillColor(kGray));
   fPosteriorPdf->plotOn(plot);
   plot->GetYaxis()->SetTitle("posterior probability");
   
   return plot; 
}


SimpleInterval* BayesianCalculator::GetInterval() const
{
  /// returns a SimpleInterval with the lower/upper limit on 
  /// the scanned variable (the parameter of interest specified in the constructor).
  /// The returned interval is a central interval with the confidence level specified  
  /// previously in SetConfidenceLevel (default is 0.95).
  /// NOTE1: for finding only an upper/lower limit of 95 % the CL must be set to 0.90
  /// NOTE2: The method can result very slow when nuisance parameters are present due to 
  ///        the time needed for performing multi-dimensional numerical integration. 
  ///        In these case using the MCMCCalculator could be more convenient.

  if (fInterval) return fInterval; 

  RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() ); 
  assert(poi);

  if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf();

  RooAbsReal* cdf = fPosteriorPdf->createCdf(fPOI);

  RooAbsFunc* cdf_bind = cdf->bindVars(fPOI,&fPOI);
  RooBrentRootFinder brf(*cdf_bind);
  brf.setTol(0.00005); // precision

  double lowerLimit = 0; 
  double upperLimit = 0; 

  double tmpVal = poi->getVal(); // patch because findRoot changes the value of poi

  double y = fSize/2;
  brf.findRoot(lowerLimit,poi->getMin(),poi->getMax(),y);

  y=1-fSize/2;
  brf.findRoot(upperLimit,poi->getMin(),poi->getMax(),y);

  poi->setVal(tmpVal); // patch: restore the original value of poi

  delete cdf_bind;
  delete cdf;

  TString interval_name = TString("BayesianInterval_a") + TString(this->GetName());
  fInterval = new SimpleInterval(interval_name,*poi,lowerLimit,upperLimit,ConfidenceLevel());
  fInterval->SetTitle("SimpleInterval from BayesianCalculator");

  return fInterval;
}

} // end namespace RooStats

 BayesianCalculator.cxx:1
 BayesianCalculator.cxx:2
 BayesianCalculator.cxx:3
 BayesianCalculator.cxx:4
 BayesianCalculator.cxx:5
 BayesianCalculator.cxx:6
 BayesianCalculator.cxx:7
 BayesianCalculator.cxx:8
 BayesianCalculator.cxx:9
 BayesianCalculator.cxx:10
 BayesianCalculator.cxx:11
 BayesianCalculator.cxx:12
 BayesianCalculator.cxx:13
 BayesianCalculator.cxx:14
 BayesianCalculator.cxx:15
 BayesianCalculator.cxx:16
 BayesianCalculator.cxx:17
 BayesianCalculator.cxx:18
 BayesianCalculator.cxx:19
 BayesianCalculator.cxx:20
 BayesianCalculator.cxx:21
 BayesianCalculator.cxx:22
 BayesianCalculator.cxx:23
 BayesianCalculator.cxx:24
 BayesianCalculator.cxx:25
 BayesianCalculator.cxx:26
 BayesianCalculator.cxx:27
 BayesianCalculator.cxx:28
 BayesianCalculator.cxx:29
 BayesianCalculator.cxx:30
 BayesianCalculator.cxx:31
 BayesianCalculator.cxx:32
 BayesianCalculator.cxx:33
 BayesianCalculator.cxx:34
 BayesianCalculator.cxx:35
 BayesianCalculator.cxx:36
 BayesianCalculator.cxx:37
 BayesianCalculator.cxx:38
 BayesianCalculator.cxx:39
 BayesianCalculator.cxx:40
 BayesianCalculator.cxx:41
 BayesianCalculator.cxx:42
 BayesianCalculator.cxx:43
 BayesianCalculator.cxx:44
 BayesianCalculator.cxx:45
 BayesianCalculator.cxx:46
 BayesianCalculator.cxx:47
 BayesianCalculator.cxx:48
 BayesianCalculator.cxx:49
 BayesianCalculator.cxx:50
 BayesianCalculator.cxx:51
 BayesianCalculator.cxx:52
 BayesianCalculator.cxx:53
 BayesianCalculator.cxx:54
 BayesianCalculator.cxx:55
 BayesianCalculator.cxx:56
 BayesianCalculator.cxx:57
 BayesianCalculator.cxx:58
 BayesianCalculator.cxx:59
 BayesianCalculator.cxx:60
 BayesianCalculator.cxx:61
 BayesianCalculator.cxx:62
 BayesianCalculator.cxx:63
 BayesianCalculator.cxx:64
 BayesianCalculator.cxx:65
 BayesianCalculator.cxx:66
 BayesianCalculator.cxx:67
 BayesianCalculator.cxx:68
 BayesianCalculator.cxx:69
 BayesianCalculator.cxx:70
 BayesianCalculator.cxx:71
 BayesianCalculator.cxx:72
 BayesianCalculator.cxx:73
 BayesianCalculator.cxx:74
 BayesianCalculator.cxx:75
 BayesianCalculator.cxx:76
 BayesianCalculator.cxx:77
 BayesianCalculator.cxx:78
 BayesianCalculator.cxx:79
 BayesianCalculator.cxx:80
 BayesianCalculator.cxx:81
 BayesianCalculator.cxx:82
 BayesianCalculator.cxx:83
 BayesianCalculator.cxx:84
 BayesianCalculator.cxx:85
 BayesianCalculator.cxx:86
 BayesianCalculator.cxx:87
 BayesianCalculator.cxx:88
 BayesianCalculator.cxx:89
 BayesianCalculator.cxx:90
 BayesianCalculator.cxx:91
 BayesianCalculator.cxx:92
 BayesianCalculator.cxx:93
 BayesianCalculator.cxx:94
 BayesianCalculator.cxx:95
 BayesianCalculator.cxx:96
 BayesianCalculator.cxx:97
 BayesianCalculator.cxx:98
 BayesianCalculator.cxx:99
 BayesianCalculator.cxx:100
 BayesianCalculator.cxx:101
 BayesianCalculator.cxx:102
 BayesianCalculator.cxx:103
 BayesianCalculator.cxx:104
 BayesianCalculator.cxx:105
 BayesianCalculator.cxx:106
 BayesianCalculator.cxx:107
 BayesianCalculator.cxx:108
 BayesianCalculator.cxx:109
 BayesianCalculator.cxx:110
 BayesianCalculator.cxx:111
 BayesianCalculator.cxx:112
 BayesianCalculator.cxx:113
 BayesianCalculator.cxx:114
 BayesianCalculator.cxx:115
 BayesianCalculator.cxx:116
 BayesianCalculator.cxx:117
 BayesianCalculator.cxx:118
 BayesianCalculator.cxx:119
 BayesianCalculator.cxx:120
 BayesianCalculator.cxx:121
 BayesianCalculator.cxx:122
 BayesianCalculator.cxx:123
 BayesianCalculator.cxx:124
 BayesianCalculator.cxx:125
 BayesianCalculator.cxx:126
 BayesianCalculator.cxx:127
 BayesianCalculator.cxx:128
 BayesianCalculator.cxx:129
 BayesianCalculator.cxx:130
 BayesianCalculator.cxx:131
 BayesianCalculator.cxx:132
 BayesianCalculator.cxx:133
 BayesianCalculator.cxx:134
 BayesianCalculator.cxx:135
 BayesianCalculator.cxx:136
 BayesianCalculator.cxx:137
 BayesianCalculator.cxx:138
 BayesianCalculator.cxx:139
 BayesianCalculator.cxx:140
 BayesianCalculator.cxx:141
 BayesianCalculator.cxx:142
 BayesianCalculator.cxx:143
 BayesianCalculator.cxx:144
 BayesianCalculator.cxx:145
 BayesianCalculator.cxx:146
 BayesianCalculator.cxx:147
 BayesianCalculator.cxx:148
 BayesianCalculator.cxx:149
 BayesianCalculator.cxx:150
 BayesianCalculator.cxx:151
 BayesianCalculator.cxx:152
 BayesianCalculator.cxx:153
 BayesianCalculator.cxx:154
 BayesianCalculator.cxx:155
 BayesianCalculator.cxx:156
 BayesianCalculator.cxx:157
 BayesianCalculator.cxx:158
 BayesianCalculator.cxx:159
 BayesianCalculator.cxx:160
 BayesianCalculator.cxx:161
 BayesianCalculator.cxx:162
 BayesianCalculator.cxx:163
 BayesianCalculator.cxx:164
 BayesianCalculator.cxx:165
 BayesianCalculator.cxx:166
 BayesianCalculator.cxx:167
 BayesianCalculator.cxx:168
 BayesianCalculator.cxx:169
 BayesianCalculator.cxx:170
 BayesianCalculator.cxx:171
 BayesianCalculator.cxx:172
 BayesianCalculator.cxx:173
 BayesianCalculator.cxx:174
 BayesianCalculator.cxx:175
 BayesianCalculator.cxx:176
 BayesianCalculator.cxx:177
 BayesianCalculator.cxx:178
 BayesianCalculator.cxx:179
 BayesianCalculator.cxx:180
 BayesianCalculator.cxx:181
 BayesianCalculator.cxx:182
 BayesianCalculator.cxx:183
 BayesianCalculator.cxx:184
 BayesianCalculator.cxx:185
 BayesianCalculator.cxx:186
 BayesianCalculator.cxx:187
 BayesianCalculator.cxx:188
 BayesianCalculator.cxx:189
 BayesianCalculator.cxx:190
 BayesianCalculator.cxx:191
 BayesianCalculator.cxx:192
 BayesianCalculator.cxx:193
 BayesianCalculator.cxx:194
 BayesianCalculator.cxx:195
 BayesianCalculator.cxx:196
 BayesianCalculator.cxx:197
 BayesianCalculator.cxx:198
 BayesianCalculator.cxx:199
 BayesianCalculator.cxx:200
 BayesianCalculator.cxx:201
 BayesianCalculator.cxx:202
 BayesianCalculator.cxx:203
 BayesianCalculator.cxx:204
 BayesianCalculator.cxx:205
 BayesianCalculator.cxx:206
 BayesianCalculator.cxx:207
 BayesianCalculator.cxx:208
 BayesianCalculator.cxx:209
 BayesianCalculator.cxx:210
 BayesianCalculator.cxx:211
 BayesianCalculator.cxx:212
 BayesianCalculator.cxx:213
 BayesianCalculator.cxx:214
 BayesianCalculator.cxx:215
 BayesianCalculator.cxx:216
 BayesianCalculator.cxx:217
 BayesianCalculator.cxx:218
 BayesianCalculator.cxx:219
 BayesianCalculator.cxx:220
 BayesianCalculator.cxx:221
 BayesianCalculator.cxx:222
 BayesianCalculator.cxx:223
 BayesianCalculator.cxx:224
 BayesianCalculator.cxx:225
 BayesianCalculator.cxx:226
 BayesianCalculator.cxx:227
 BayesianCalculator.cxx:228
 BayesianCalculator.cxx:229
 BayesianCalculator.cxx:230
 BayesianCalculator.cxx:231
 BayesianCalculator.cxx:232
 BayesianCalculator.cxx:233
 BayesianCalculator.cxx:234
 BayesianCalculator.cxx:235
 BayesianCalculator.cxx:236
 BayesianCalculator.cxx:237
 BayesianCalculator.cxx:238
 BayesianCalculator.cxx:239
 BayesianCalculator.cxx:240
 BayesianCalculator.cxx:241
 BayesianCalculator.cxx:242
 BayesianCalculator.cxx:243
 BayesianCalculator.cxx:244
 BayesianCalculator.cxx:245
 BayesianCalculator.cxx:246
 BayesianCalculator.cxx:247
 BayesianCalculator.cxx:248
 BayesianCalculator.cxx:249
 BayesianCalculator.cxx:250
 BayesianCalculator.cxx:251
 BayesianCalculator.cxx:252
 BayesianCalculator.cxx:253
 BayesianCalculator.cxx:254
 BayesianCalculator.cxx:255
 BayesianCalculator.cxx:256
 BayesianCalculator.cxx:257
 BayesianCalculator.cxx:258
 BayesianCalculator.cxx:259