ROOT  6.06/09
Reference Guide
BayesianCalculator.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3 /*************************************************************************
4  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5  * All rights reserved. *
6  * *
7  * For the licensing terms see $ROOTSYS/LICENSE. *
8  * For the list of contributors see $ROOTSYS/README/CREDITS. *
9  *************************************************************************/
10 
11 ////////////////////////////////////////////////////////////////////////////////
12 /// BayesianCalculator implementation
13 ///////////////////////////////////////////////////////////////////////////
14 
15 
16 // include other header files
17 
18 #include "RooAbsFunc.h"
19 #include "RooAbsReal.h"
20 #include "RooRealVar.h"
21 #include "RooArgSet.h"
22 #include "RooBrentRootFinder.h"
23 #include "RooFormulaVar.h"
24 #include "RooGenericPdf.h"
25 #include "RooPlot.h"
26 #include "RooProdPdf.h"
27 #include "RooDataSet.h"
28 
29 // include header file of this class
31 #include "RooStats/ModelConfig.h"
32 #include "RooStats/RooStatsUtils.h"
33 
34 #include "Math/IFunction.h"
36 #include "Math/Integrator.h"
37 #include "Math/RootFinder.h"
38 #include "Math/BrentMinimizer1D.h"
39 #include "RooFunctor.h"
40 #include "RooFunctor1DBinding.h"
41 #include "RooTFnBinding.h"
42 #include "RooMsgService.h"
43 
44 #include "TAxis.h"
45 #include "TF1.h"
46 #include "TH1.h"
47 #include "TMath.h"
48 #include "TCanvas.h"
49 
50 #include <map>
51 #include <cmath>
52 
53 //#include "TRandom.h"
54 #include "RConfigure.h"
55 
57 
58 using namespace std;
59 
60 namespace RooStats {
61 
62 
63 // first some utility classes and functions
64 
65 #ifdef R__HAS_MATHMORE
67 #else
69 #endif
70 
71 
72 
73 
74 struct LikelihoodFunction {
75  LikelihoodFunction(RooFunctor & f, RooFunctor * prior = 0, double offset = 0) :
76  fFunc(f), fPrior(prior),
77  fOffset(offset), fMaxL(0) {
78  fFunc.binding().resetNumCall();
79  }
80 
81  void SetPrior(RooFunctor * prior) { fPrior = prior; }
82 
83  double operator() (const double *x ) const {
84  double nll = fFunc(x) - fOffset;
85  double likelihood = std::exp(-nll);
86 
87  if (fPrior) likelihood *= (*fPrior)(x);
88 
89  int nCalls = fFunc.binding().numCall();
90  if (nCalls > 0 && nCalls % 1000 == 0) {
91  ooccoutD((TObject*)0,Eval) << "Likelihood evaluation ncalls = " << nCalls
92  << " x0 " << x[0] << " nll = " << nll+fOffset;
93  if (fPrior) ooccoutD((TObject*)0,Eval) << " prior(x) = " << (*fPrior)(x);
94  ooccoutD((TObject*)0,Eval) << " likelihood " << likelihood
95  << " max Likelihood " << fMaxL << std::endl;
96  }
97 
98  if (likelihood > fMaxL ) {
99  fMaxL = likelihood;
100  if ( likelihood > 1.E10) {
101  ooccoutW((TObject*)0,Eval) << "LikelihoodFunction::() WARNING - Huge likelihood value found for parameters ";
102  for (int i = 0; i < fFunc.nObs(); ++i)
103  ooccoutW((TObject*)0,Eval) << " x[" << i << " ] = " << x[i];
104  ooccoutW((TObject*)0,Eval) << " nll = " << nll << " L = " << likelihood << std::endl;
105  }
106  }
107 
108  return likelihood;
109  }
110 
111  // for the 1D case
112  double operator() (double x) const {
113  // just call the previous method
114  assert(fFunc.nObs() == 1); // check nobs = 1
115  double tmp = x;
116  return (*this)(&tmp);
117  }
118 
119  RooFunctor & fFunc; // functor representing the nll function
120  RooFunctor * fPrior; // functor representing the prior function
121  double fOffset; // offset used to bring the nll in a reasanble range for computing the exponent
122  mutable double fMaxL;
123 };
124 
125 
126 // Posterior CDF Function class
127 // for integral of posterior function in nuisance and POI
128 // 1-Dim function as function of the poi
129 
130 class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
131 
132 public:
133 
134  PosteriorCdfFunction(RooAbsReal & nll, RooArgList & bindParams, RooAbsReal * prior = 0, const char * integType = 0, double nllMinimum = 0) :
135  fFunctor(nll, bindParams, RooArgList() ), // functor
136  fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
137  fLikelihood(fFunctor, 0, nllMinimum), // integral of exp(-nll) function
138  fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ), // integrator
139  fXmin(bindParams.getSize() ), // vector of parameters (min values)
140  fXmax(bindParams.getSize() ), // vector of parameter (max values)
141  fNorm(1.0), fNormErr(0.0), fOffset(0), fMaxPOI(0),
142  fHasNorm(false), fUseOldValues(true), fError(false)
143  {
144 
145  if (prior) {
146  fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, bindParams, RooArgList() ));
147  fLikelihood.SetPrior(fPriorFunc.get() );
148  }
149 
150  fIntegrator.SetFunction(fLikelihood, bindParams.getSize() );
151 
152  ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. "
153  << " nllMinimum is " << nllMinimum << std::endl;
154 
155  std::vector<double> par(bindParams.getSize());
156  for (unsigned int i = 0; i < fXmin.size(); ++i) {
157  RooRealVar & var = (RooRealVar &) bindParams[i];
158  fXmin[i] = var.getMin();
159  fXmax[i] = var.getMax();
160  par[i] = var.getVal();
161  ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Integrate" << var.GetName()
162  << " in interval [ " << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
163  }
164 
165  fIntegrator.Options().Print(ooccoutD((TObject*)0,NumIntegration));
166 
167  // store max POI value because it will be changed when evaluating the function
168  fMaxPOI = fXmax[0];
169 
170  // compute first the normalization with the poi
171  fNorm = (*this)( fMaxPOI );
172  if (fError) ooccoutE((TObject*)0,NumIntegration) << "PosteriorFunction::Error computing normalization - norm = " << fNorm << std::endl;
173  fHasNorm = true;
174  fNormCdfValues.insert(std::make_pair(fXmin[0], 0) );
175  fNormCdfValues.insert(std::make_pair(fXmax[0], 1.0) );
176 
177  }
178 
179  // copy constructor (needed for Cloning the object)
180  // need special treatment because integrator
181  // has no copy constuctor
182  PosteriorCdfFunction(const PosteriorCdfFunction & rhs) :
183  ROOT::Math::IGenFunction(),
184  fFunctor(rhs.fFunctor),
185  //fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
186  fPriorFunc(rhs.fPriorFunc),
187  fLikelihood(fFunctor, fPriorFunc.get(), rhs.fLikelihood.fOffset),
188  fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ), // integrator
189  fXmin( rhs.fXmin),
190  fXmax( rhs.fXmax),
191  fNorm( rhs.fNorm),
192  fNormErr( rhs.fNormErr),
193  fOffset(rhs.fOffset),
194  fMaxPOI(rhs.fMaxPOI),
195  fHasNorm(rhs.fHasNorm),
196  fUseOldValues(rhs.fUseOldValues),
197  fError(rhs.fError),
198  fNormCdfValues(rhs.fNormCdfValues)
199  {
200  fIntegrator.SetFunction(fLikelihood, fXmin.size() );
201  // need special treatment for the smart pointer
202  // if (rhs.fPriorFunc.get() ) {
203  // fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*(rhs.fPriorFunc) ) );
204  // fLikelihood.SetPrior( fPriorFunc.get() );
205  // }
206  }
207 
208 
209  bool HasError() const { return fError; }
210 
211 
212  ROOT::Math::IGenFunction * Clone() const {
213  ooccoutD((TObject*)0,NumIntegration) << " cloning function .........." << std::endl;
214  return new PosteriorCdfFunction(*this);
215  }
216 
217  // offset value for computing the root
218  void SetOffset(double offset) { fOffset = offset; }
219 
220 private:
221 
222  // make assignment operator private
223  PosteriorCdfFunction& operator=(const PosteriorCdfFunction &) {
224  return *this;
225  }
226 
227  double DoEval (double x) const {
228 
229  // evaluate cdf at poi value x by integrating poi from [xmin,x] and all the nuisances
230  fXmax[0] = x;
231  if (x <= fXmin[0] ) return -fOffset;
232  // could also avoid a function evaluation at maximum
233  if (x >= fMaxPOI && fHasNorm) return 1. - fOffset; // cdf is bound to these values
234 
235  // computes the integral using a previous cdf estimate
236  double normcdf0 = 0;
237  if (fHasNorm && fUseOldValues) {
238  // look in the map of the stored cdf values the closes one
239  std::map<double,double>::iterator itr = fNormCdfValues.upper_bound(x);
240  itr--; // upper bound returns a poistion 1 up of the value we want
241  if (itr != fNormCdfValues.end() ) {
242  fXmin[0] = itr->first;
243  normcdf0 = itr->second;
244  // ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction: computing integral between in poi interval : "
245  // << fXmin[0] << " - " << fXmax[0] << std::endl;
246  }
247  }
248 
249  fFunctor.binding().resetNumCall(); // reset number of calls for debug
250 
251  double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]);
252  double error = fIntegrator.Error();
253  double normcdf = cdf/fNorm; // normalize the cdf
254 
255  ooccoutD((TObject*)0,NumIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , "
256  << fXmax[0] << "] integral = " << cdf << " +/- " << error
257  << " norm-integ = " << normcdf << " cdf(x) = " << normcdf+normcdf0
258  << " ncalls = " << fFunctor.binding().numCall() << std::endl;
259 
260  if (TMath::IsNaN(cdf) || cdf > std::numeric_limits<double>::max()) {
261  ooccoutE((TObject*)0,NumIntegration) << "PosteriorFunction::Error computing integral - cdf = "
262  << cdf << std::endl;
263  fError = true;
264  }
265 
266  if (cdf != 0 && error/cdf > 0.2 )
267  oocoutW((TObject*)0,NumIntegration) << "PosteriorCdfFunction: integration error is larger than 20 % x0 = " << fXmin[0]
268  << " x = " << x << " cdf(x) = " << cdf << " +/- " << error << std::endl;
269 
270  if (!fHasNorm) {
271  oocoutI((TObject*)0,NumIntegration) << "PosteriorCdfFunction - integral of posterior = "
272  << cdf << " +/- " << error << std::endl;
273  fNormErr = error;
274  return cdf;
275  }
276 
277  normcdf += normcdf0;
278 
279  // store values in the map
280  if (fUseOldValues) {
281  fNormCdfValues.insert(std::make_pair(x, normcdf) );
282  }
283 
284  double errnorm = sqrt( error*error + normcdf*normcdf * fNormErr * fNormErr )/fNorm;
285  if (normcdf > 1. + 3 * errnorm) {
286  oocoutW((TObject*)0,NumIntegration) << "PosteriorCdfFunction: normalized cdf values is larger than 1"
287  << " x = " << x << " normcdf(x) = " << normcdf << " +/- " << error/fNorm << std::endl;
288  }
289 
290  return normcdf - fOffset; // apply an offset (for finding the roots)
291  }
292 
293  mutable RooFunctor fFunctor; // functor binding nll
294  mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
295  LikelihoodFunction fLikelihood; // likelihood function
296  mutable ROOT::Math::IntegratorMultiDim fIntegrator; // integrator (mutable because Integral() is not const
297  mutable std::vector<double> fXmin; // min value of parameters (poi+nuis) -
298  mutable std::vector<double> fXmax; // max value of parameters (poi+nuis) - max poi changes so it is mutable
299  double fNorm; // normalization value (computed in ctor)
300  mutable double fNormErr; // normalization error value (computed in ctor)
301  double fOffset; // offset for computing the root
302  double fMaxPOI; // maximum value of POI
303  bool fHasNorm; // flag to control first call to the function
304  bool fUseOldValues; // use old cdf values
305  mutable bool fError; // flag to indicate if a numerical evaluation error occurred
306  mutable std::map<double,double> fNormCdfValues;
307 };
308 
309 //__________________________________________________________________
310 // Posterior Function class
311 // 1-Dim function as function of the poi
312 // and it integrated all the nuisance parameters
313 
314 class PosteriorFunction : public ROOT::Math::IGenFunction {
315 
316 public:
317 
318 
319  PosteriorFunction(RooAbsReal & nll, RooRealVar & poi, RooArgList & nuisParams, RooAbsReal * prior = 0, const char * integType = 0, double
320  norm = 1.0, double nllOffset = 0, int niter = 0) :
321  fFunctor(nll, nuisParams, RooArgList() ),
322  fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
323  fLikelihood(fFunctor, 0, nllOffset),
324  fPoi(&poi),
325  fXmin(nuisParams.getSize() ),
326  fXmax(nuisParams.getSize() ),
327  fNorm(norm),
328  fError(0)
329  {
330 
331  if (prior) {
332  fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, nuisParams, RooArgList() ));
333  fLikelihood.SetPrior(fPriorFunc.get() );
334  }
335 
336  ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Evaluate the posterior function by integrating the nuisances: " << std::endl;
337  for (unsigned int i = 0; i < fXmin.size(); ++i) {
338  RooRealVar & var = (RooRealVar &) nuisParams[i];
339  fXmin[i] = var.getMin();
340  fXmax[i] = var.getMax();
341  ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction::Integrate " << var.GetName()
342  << " in interval [" << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
343  }
344  if (fXmin.size() == 1) { // 1D case
345  fIntegratorOneDim.reset( new ROOT::Math::Integrator(ROOT::Math::IntegratorOneDim::GetType(integType) ) );
346 
347  fIntegratorOneDim->SetFunction(fLikelihood);
348  // interested only in relative tolerance
349  //fIntegratorOneDim->SetAbsTolerance(1.E-300);
350  fIntegratorOneDim->Options().Print(ooccoutD((TObject*)0,NumIntegration) );
351  }
352  else if (fXmin.size() > 1) { // multiDim case
353  fIntegratorMultiDim.reset(new ROOT::Math::IntegratorMultiDim(ROOT::Math::IntegratorMultiDim::GetType(integType) ) );
354  fIntegratorMultiDim->SetFunction(fLikelihood, fXmin.size());
355  ROOT::Math::IntegratorMultiDimOptions opt = fIntegratorMultiDim->Options();
356  if (niter > 0) {
357  opt.SetNCalls(niter);
358  fIntegratorMultiDim->SetOptions(opt);
359  }
360  //fIntegratorMultiDim->SetAbsTolerance(1.E-300);
361  // print the options
363  }
364  }
365 
366 
367  ROOT::Math::IGenFunction * Clone() const {
368  assert(1);
369  return 0; // cannot clone this function for integrator
370  }
371 
372  double Error() const { return fError;}
373 
374 
375 private:
376  double DoEval (double x) const {
377 
378  // evaluate posterior function at a poi value x by integrating all nuisance parameters
379 
380  fPoi->setVal(x);
381  fFunctor.binding().resetNumCall(); // reset number of calls for debug
382 
383  double f = 0;
384  double error = 0;
385  if (fXmin.size() == 1) { // 1D case
386  f = fIntegratorOneDim->Integral(fXmin[0],fXmax[0]);
387  error = fIntegratorOneDim->Error();
388  }
389  else if (fXmin.size() > 1) { // multi-dim case
390  f = fIntegratorMultiDim->Integral(&fXmin[0],&fXmax[0]);
391  error = fIntegratorMultiDim->Error();
392  } else {
393  // no integration to be done
394  f = fLikelihood(x);
395  }
396 
397  // debug
398  ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunction: POI value = "
399  << x << "\tf(x) = " << f << " +/- " << error
400  << " norm-f(x) = " << f/fNorm
401  << " ncalls = " << fFunctor.binding().numCall() << std::endl;
402 
403 
404 
405 
406  if (f != 0 && error/f > 0.2 )
407  ooccoutW((TObject*)0,NumIntegration) << "PosteriorFunction::DoEval - Error from integration in "
408  << fXmin.size() << " Dim is larger than 20 % "
409  << "x = " << x << " p(x) = " << f << " +/- " << error << std::endl;
410 
411  fError = error / fNorm;
412  return f / fNorm;
413  }
414 
415  mutable RooFunctor fFunctor;
416  mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
417  LikelihoodFunction fLikelihood;
418  RooRealVar * fPoi;
419  std::unique_ptr<ROOT::Math::Integrator> fIntegratorOneDim;
420  std::unique_ptr<ROOT::Math::IntegratorMultiDim> fIntegratorMultiDim;
421  std::vector<double> fXmin;
422  std::vector<double> fXmax;
423  double fNorm;
424  mutable double fError;
425 };
426 
427 //////////////////////////////////////////////////////////////////////////////////////////////////////////
428 // Posterior function obtaining sampling toy MC for th enuisance according to their pdf
429 class PosteriorFunctionFromToyMC : public ROOT::Math::IGenFunction {
430 
431 public:
432 
433 
434  PosteriorFunctionFromToyMC(RooAbsReal & nll, RooAbsPdf & pdf, RooRealVar & poi, RooArgList & nuisParams, RooAbsReal * prior = 0, double
435  nllOffset = 0, int niter = 0, bool redoToys = true ) :
436  fFunctor(nll, nuisParams, RooArgList() ),
437  fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
438  fLikelihood(fFunctor, 0, nllOffset),
439  fPdf(&pdf),
440  fPoi(&poi),
441  fNuisParams(nuisParams),
442  fGenParams(0),
443  fNumIterations(niter),
444  fError(-1),
445  fRedoToys(redoToys)
446  {
447  if (niter == 0) fNumIterations = 100; // default value
448 
449  if (prior) {
450  fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*prior, nuisParams, RooArgList() ));
451  fLikelihood.SetPrior(fPriorFunc.get() );
452  }
453 
454  ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Evaluate the posterior function by randomizing the nuisances: niter " << fNumIterations << std::endl;
455 
456  ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Pdf used for randomizing the nuisance is " << fPdf->GetName() << std::endl;
457  // check that pdf contains the nuisance
458  RooArgSet * vars = fPdf->getVariables();
459  for (int i = 0; i < fNuisParams.getSize(); ++i) {
460  if (!vars->find( fNuisParams[i].GetName() ) ) {
461  ooccoutW((TObject*)0,InputArguments) << "Nuisance parameter " << fNuisParams[i].GetName()
462  << " is not part of sampling pdf. "
463  << "they will be trated as constant " << std::endl;
464  }
465  }
466  delete vars;
467 
468  if (!fRedoToys) {
469  ooccoutI((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC::Generate nuisance toys only one time (for all POI points)" << std::endl;
470  GenerateToys();
471  }
472  }
473 
474  virtual ~PosteriorFunctionFromToyMC() { if (fGenParams) delete fGenParams; }
475 
476  // generate first n-samples of the nuisance parameters
477  void GenerateToys() const {
478  if (fGenParams) delete fGenParams;
479  fGenParams = fPdf->generate(fNuisParams, fNumIterations);
480  if(fGenParams==0) {
481  ooccoutE((TObject*)0,InputArguments) << "PosteriorFunctionFromToyMC - failed to generate nuisance parameters" << std::endl;
482  }
483  }
484 
485  double Error() const { return fError;}
486 
487  ROOT::Math::IGenFunction * Clone() const {
488  // use defsult copy constructor
489  //return new PosteriorFunctionFromToyMC(*this);
490  // clone not implemented
491  assert(1);
492  return 0;
493  }
494 
495 private:
496  // evaluate the posterior at the poi value x
497  double DoEval( double x) const {
498 
499  int npar = fNuisParams.getSize();
500  assert (npar > 0);
501 
502 
503  // generate the toys
504  if (fRedoToys) GenerateToys();
505  if (!fGenParams) return 0;
506 
507  // evaluate posterior function at a poi value x by integrating all nuisance parameters
508 
509  fPoi->setVal(x);
510 
511  // loop over all of the generate data
512  double sum = 0;
513  double sum2 = 0;
514 
515  for(int iter=0; iter<fNumIterations; ++iter) {
516 
517  // get the set of generated parameters and set the nuisance parameters to the generated values
518  std::vector<double> p(npar);
519  for (int i = 0; i < npar; ++i) {
520  const RooArgSet* genset=fGenParams->get(iter);
521  RooAbsArg * arg = genset->find( fNuisParams[i].GetName() );
522  RooRealVar * var = dynamic_cast<RooRealVar*>(arg);
523  assert(var != 0);
524  p[i] = var->getVal();
525  ((RooRealVar &) fNuisParams[i]).setVal(p[i]);
526  }
527 
528  // evaluate now the likelihood function
529  double fval = fLikelihood( &p.front() );
530 
531  // liklihood already must contained the pdf we have sampled
532  // so we must divided by it. The value must be normalized on all
533  // other parameters
534  RooArgSet arg(fNuisParams);
535  double nuisPdfVal = fPdf->getVal(&arg);
536  fval /= nuisPdfVal;
537 
538 
539  if( fval > std::numeric_limits<double>::max() ) {
540  ooccoutE((TObject*)0,Eval) << "BayesianCalculator::EvalPosteriorFunctionFromToy : "
541  << "Likelihood evaluates to infinity " << std::endl;
542  ooccoutE((TObject*)0,Eval) << "poi value = " << x << std::endl;
543  ooccoutE((TObject*)0,Eval) << "Nuisance parameter values : ";
544  for (int i = 0; i < npar; ++i)
545  ooccoutE((TObject*)0,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
546  ooccoutE((TObject*)0,Eval) << " - return 0 " << std::endl;
547 
548  fError = 1.E30;
549  return 0;
550  }
551  if( TMath::IsNaN(fval) ) {
552  ooccoutE((TObject*)0,Eval) << "BayesianCalculator::EvalPosteriorFunctionFromToy : "
553  << "Likelihood is a NaN " << std::endl;
554  ooccoutE((TObject*)0,Eval) << "poi value = " << x << std::endl;
555  ooccoutE((TObject*)0,Eval) << "Nuisance parameter values : ";
556  for (int i = 0; i < npar; ++i)
557  ooccoutE((TObject*)0,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
558  ooccoutE((TObject*)0,Eval) << " - return 0 " << std::endl;
559  fError = 1.E30;
560  return 0;
561  }
562 
563 
564 
565  sum += fval;
566  sum2 += fval*fval;
567  }
568 
569  // compute the average and variance
570  double val = sum/double(fNumIterations);
571  double dval2 = std::max( sum2/double(fNumIterations) - val*val, 0.0);
572  fError = std::sqrt( dval2 / fNumIterations);
573 
574  // debug
575  ooccoutD((TObject*)0,NumIntegration) << "PosteriorFunctionFromToyMC: POI value = "
576  << x << "\tp(x) = " << val << " +/- " << fError << std::endl;
577 
578 
579  if (val != 0 && fError/val > 0.2 ) {
580  ooccoutW((TObject*)0,NumIntegration) << "PosteriorFunctionFromToyMC::DoEval"
581  << " - Error in estimating posterior is larger than 20% ! "
582  << "x = " << x << " p(x) = " << val << " +/- " << fError << std::endl;
583  }
584 
585 
586  return val;
587  }
588 
589  mutable RooFunctor fFunctor;
590  mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
591  LikelihoodFunction fLikelihood;
592  mutable RooAbsPdf * fPdf;
593  RooRealVar * fPoi;
594  RooArgList fNuisParams;
595  mutable RooDataSet * fGenParams;
596  int fNumIterations;
597  mutable double fError;
598  bool fRedoToys; // do toys every iteration
599 
600 };
601 
602 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
603 // Implementation of BayesianCalculator
604 /////////////////////////////////////////////////////////////////////////////////////////////////////
605 
606 BayesianCalculator::BayesianCalculator() :
607  fData(0),
608  fPdf(0),
609  fPriorPdf(0),
610  fNuisancePdf(0),
611  fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
612  fPosteriorFunction(0), fApproxPosterior(0),
613  fLower(0), fUpper(0),
614  fNLLMin(0),
615  fSize(0.05), fLeftSideFraction(0.5),
616  fBrfPrecision(0.00005),
617  fNScanBins(-1),
618  fNumIterations(0),
619  fValidInterval(false)
620 {
621  // default constructor
622 }
623 
624 BayesianCalculator::BayesianCalculator( /* const char* name, const char* title, */
625  RooAbsData& data,
626  RooAbsPdf& pdf,
627  const RooArgSet& POI,
628  RooAbsPdf& priorPdf,
629  const RooArgSet* nuisanceParameters ) :
630  fData(&data),
631  fPdf(&pdf),
632  fPOI(POI),
633  fPriorPdf(&priorPdf),
634  fNuisancePdf(0),
635  fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
636  fPosteriorFunction(0), fApproxPosterior(0),
637  fLower(0), fUpper(0),
638  fNLLMin(0),
639  fSize(0.05), fLeftSideFraction(0.5),
640  fBrfPrecision(0.00005),
641  fNScanBins(-1),
642  fNumIterations(0),
643  fValidInterval(false)
644 {
645  // Constructor from data set, model pdf, parameter of interests and prior pdf
646  // If nuisance parameters are given they will be integrated according either to the prior or
647  // their constraint term included in the model
648 
649  if (nuisanceParameters) fNuisanceParameters.add(*nuisanceParameters);
650  // remove constant nuisance parameters
652 }
653 
655  ModelConfig & model) :
656  fData(&data),
657  fPdf(model.GetPdf()),
658  fPriorPdf( model.GetPriorPdf()),
659  fNuisancePdf(0),
660  fProductPdf (0), fLogLike(0), fLikelihood (0), fIntegratedLikelihood (0), fPosteriorPdf(0),
661  fPosteriorFunction(0), fApproxPosterior(0),
662  fLower(0), fUpper(0),
663  fNLLMin(0),
664  fSize(0.05), fLeftSideFraction(0.5),
665  fBrfPrecision(0.00005),
666  fNScanBins(-1),
667  fNumIterations(0),
668  fValidInterval(false)
669 {
670  // Constructor from a data set and a ModelConfig
671  // model pdf, poi and nuisances will be taken from the ModelConfig
672  SetModel(model);
673 }
674 
675 
677 {
678  // destructor
679  ClearAll();
680 }
681 
683  // clear all cached pdf objects
684  if (fProductPdf) delete fProductPdf;
685  if (fLogLike) delete fLogLike;
686  if (fLikelihood) delete fLikelihood;
688  if (fPosteriorPdf) delete fPosteriorPdf;
690  if (fApproxPosterior) delete fApproxPosterior;
691  fPosteriorPdf = 0;
692  fPosteriorFunction = 0;
693  fProductPdf = 0;
694  fLogLike = 0;
695  fLikelihood = 0;
697  fLower = 0;
698  fUpper = 0;
699  fNLLMin = 0;
700  fValidInterval = false;
701 }
702 
704  // set the model to use
705  // The model pdf, prior pdf, parameter of interest and nuisances
706  // will be taken according to the model
707 
708  fPdf = model.GetPdf();
709  fPriorPdf = model.GetPriorPdf();
710  // assignment operator = does not do a real copy the sets (must use add method)
711  fPOI.removeAll();
714  if (model.GetParametersOfInterest()) fPOI.add( *(model.GetParametersOfInterest()) );
717  // remove constant nuisance parameters
719 
720  // invalidate the cached pointers
721  ClearAll();
722 }
723 
724 
725 
727 {
728  // Build and return the posterior function (not normalized) as a RooAbsReal
729  // the posterior is obtained from the product of the likelihood function and the
730  // prior pdf which is then intergated in the nuisance parameters (if existing).
731  // A prior function for the nuisance can be specified either in the prior pdf object
732  // or in the model itself. If no prior nuisance is specified, but prior parameters are then
733  // the integration is performed assuming a flat prior for the nuisance parameters.
734  // NOTE: the return object is managed by the class, users do not need to delete it
735 
737  if (fLikelihood) return fLikelihood;
738 
739  // run some sanity checks
740  if (!fPdf ) {
741  coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing pdf model" << std::endl;
742  return 0;
743  }
744  if (fPOI.getSize() == 0) {
745  coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing parameter of interest" << std::endl;
746  return 0;
747  }
748  if (fPOI.getSize() > 1) {
749  coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - current implementation works only on 1D intervals" << std::endl;
750  return 0;
751  }
752 
753 
754  RooArgSet* constrainedParams = fPdf->getParameters(*fData);
755  // remove the constant parameters
756  RemoveConstantParameters(constrainedParams);
757 
758  //constrainedParams->Print("V");
759 
760  // use RooFit::Constrain() to be sure constraints terms are taken into account
762 
763 
764 
765  ccoutD(Eval) << "BayesianCalculator::GetPosteriorFunction : "
766  << " pdf value " << fPdf->getVal()
767  << " neglogLikelihood = " << fLogLike->getVal() << std::endl;
768 
769  if (fPriorPdf)
770  ccoutD(Eval) << "\t\t\t priorPOI value " << fPriorPdf->getVal() << std::endl;
771 
772  // check that likelihood evaluation is not inifinity
773  double nllVal = fLogLike->getVal();
774  if ( nllVal > std::numeric_limits<double>::max() ) {
775  coutE(Eval) << "BayesianCalculator::GetPosteriorFunction : "
776  << " Negative log likelihood evaluates to infinity " << std::endl
777  << " Non-const Parameter values : ";
778  RooArgList p(*constrainedParams);
779  for (int i = 0; i < p.getSize(); ++i) {
780  RooRealVar * v = dynamic_cast<RooRealVar *>(&p[i] );
781  if (v!=0) ccoutE(Eval) << v->GetName() << " = " << v->getVal() << " ";
782  }
783  ccoutE(Eval) << std::endl;
784  ccoutE(Eval) << "-- Perform a full likelihood fit of the model before or set more reasanable parameter values"
785  << std::endl;
786  coutE(Eval) << "BayesianCalculator::GetPosteriorFunction : " << " cannot compute posterior function " << std::endl;
787  return 0;
788  }
789 
790 
791 
792  // need do find minimum of log-likelihood in the range to shift function
793  // to avoid numerical errors when we compute the likelihood (overflows in the exponent)
794  // N.B.: this works for only 1 parameter of interest otherwise Minuit should be used for finding the minimum
795  RooFunctor * nllFunc = fLogLike->functor(fPOI);
796  assert(nllFunc);
797  ROOT::Math::Functor1D wnllFunc(*nllFunc);
798  RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() );
799  assert(poi);
800 
801  // try to reduce some error messages
802  //Bool_t silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
804 
805 
806  coutI(Eval) << "BayesianCalculator::GetPosteriorFunction : "
807  << " nll value " << nllVal << " poi value = " << poi->getVal() << std::endl;
808 
809 
811  minim.SetFunction(wnllFunc,poi->getMin(),poi->getMax() );
812  bool ret = minim.Minimize(100,1.E-3,1.E-3);
813  fNLLMin = 0;
814  if (ret) fNLLMin = minim.FValMinimum();
815 
816  coutI(Eval) << "BayesianCalculator::GetPosteriorFunction : minimum of NLL vs POI for POI = "
817  << poi->getVal() << " min NLL = " << fNLLMin << std::endl;
818 
819  delete nllFunc;
820 
821  delete constrainedParams;
822 
823 
824  if ( fNuisanceParameters.getSize() == 0 || fIntegrationType.Contains("ROOFIT") ) {
825 
826  ccoutD(Eval) << "BayesianCalculator::GetPosteriorFunction : use ROOFIT integration "
827  << std::endl;
828 
829 #ifdef DOLATER // (not clear why this does not work)
830  // need to make in this case a likelihood from the nll and make the product with the prior
831  TString likeName = TString("likelihood_times_prior_") + TString(fPriorPdf->GetName());
832  TString formula;
833  formula.Form("exp(-@0+%f+log(@1))",fNLLMin);
834  fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike,*fPriorPdf));
835 #else
836  // here use RooProdPdf (not very nice) but working
837 
838  if (fLogLike) delete fLogLike;
839  if (fProductPdf) {
840  delete fProductPdf;
841  fProductPdf = 0;
842  }
843 
844  // // create a unique name for the product pdf
845  RooAbsPdf * pdfAndPrior = fPdf;
846  if (fPriorPdf) {
847  TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );
848  // save this as data member since it needs to be deleted afterwards
849  fProductPdf = new RooProdPdf(prodName,"",RooArgList(*fPdf,*fPriorPdf));
850  pdfAndPrior = fProductPdf;
851  }
852 
853  RooArgSet* constrParams = fPdf->getParameters(*fData);
854  // remove the constant parameters
855  RemoveConstantParameters(constrParams);
857  delete constrParams;
858 
859  TString likeName = TString("likelihood_times_prior_") + TString(pdfAndPrior->GetName());
860  TString formula;
861  formula.Form("exp(-@0+%f)",fNLLMin);
862  fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike));
863 #endif
864 
865 
866  // if no nuisance parameter we can just return the likelihood function
867  if (fNuisanceParameters.getSize() == 0) {
869  fLikelihood = 0;
870  }
871  else
872  // case of using RooFit for the integration
874 
875 
876  }
877 
878  else if ( fIntegrationType.Contains("TOYMC") ) {
879  // compute the posterior as expectation values of the likelihood function
880  // sampling on the nuisance parameters
881 
882  RooArgList nuisParams(fNuisanceParameters);
883 
884  bool doToysEveryIteration = true;
885  // if type is 1-TOYMC or TOYMC-1
886  if ( fIntegrationType.Contains("1") || fIntegrationType.Contains("ONE") ) doToysEveryIteration = false;
887 
888  RooAbsPdf * samplingPdf = (fNuisancePdf) ? fNuisancePdf : fPdf;
889  if (!fNuisancePdf) {
890  ccoutI(Eval) << "BayesianCalculator::GetPosteriorFunction : no nuisance pdf is provided, try using global pdf (this will be slower)"
891  << std::endl;
892  }
893  fPosteriorFunction = new PosteriorFunctionFromToyMC(*fLogLike, *samplingPdf, *poi, nuisParams, fPriorPdf, fNLLMin,
894  fNumIterations, doToysEveryIteration );
895 
896  TString name = "toyposteriorfunction_from_";
897  name += fLogLike->GetName();
899 
900  // need to scan likelihood in this case
901  if (fNScanBins <= 0) fNScanBins = 100;
902 
903  }
904 
905  else {
906 
907  // use ROOT integration method if there are nuisance parameters
908 
909  RooArgList nuisParams(fNuisanceParameters);
910  fPosteriorFunction = new PosteriorFunction(*fLogLike, *poi, nuisParams, fPriorPdf, fIntegrationType, 1., fNLLMin, fNumIterations );
911 
912  TString name = "posteriorfunction_from_";
913  name += fLogLike->GetName();
915 
916  }
917 
918 
919  if (RooAbsReal::numEvalErrors() > 0)
920  coutW(Eval) << "BayesianCalculator::GetPosteriorFunction : " << RooAbsReal::numEvalErrors() << " errors reported in evaluating log-likelihood function "
921  << std::endl;
924 
925  return fIntegratedLikelihood;
926 
927 }
928 
930 {
931  // Build and return the posterior pdf (i.e posterior function normalized to all range of poi)
932  // Note that an extra integration in the POI is required for the normalization
933  // NOTE: user must delete the returned object
934 
935  RooAbsReal * plike = GetPosteriorFunction();
936  if (!plike) return 0;
937 
938 
939  // create a unique name on the posterior from the names of the components
940  TString posteriorName = this->GetName() + TString("_posteriorPdf_") + plike->GetName();
941 
942  RooAbsPdf * posteriorPdf = new RooGenericPdf(posteriorName,"@0",*plike);
943 
944  return posteriorPdf;
945 }
946 
947 
948 RooPlot* BayesianCalculator::GetPosteriorPlot(bool norm, double precision ) const
949 {
950  // return a RooPlot with the posterior and the credibility region
951  // NOTE: User takes ownership of the returned object
952 
954 
955  // if a scan is requested approximate the posterior
956  if (fNScanBins > 0)
958 
959  RooAbsReal * posterior = fIntegratedLikelihood;
960  if (norm) {
961  // delete and re-do always posterior pdf (could be invalid after approximating it)
962  if (fPosteriorPdf) delete fPosteriorPdf;
964  posterior = fPosteriorPdf;
965  }
966  if (!posterior) return 0;
967 
968  if (!fValidInterval) GetInterval();
969 
970  RooAbsRealLValue* poi = dynamic_cast<RooAbsRealLValue*>( fPOI.first() );
971  assert(poi);
972 
973 
974  RooPlot* plot = poi->frame();
975  if (!plot) return 0;
976 
977  // try to reduce some error messages
979 
980  plot->SetTitle(TString("Posterior probability of parameter \"")+TString(poi->GetName())+TString("\""));
982  posterior->plotOn(plot);
983  plot->GetYaxis()->SetTitle("posterior function");
984 
985  // reset the counts and default mode
988 
989  return plot;
990 }
991 
993  // set the integration type (possible type are) :
994  // 1D integration ( used when only one nuisance and when the posterior is scanned):
995  // adaptive , gauss, nonadaptive
996  // multidim:
997  // ADAPTIVE, adaptive numerical integration
998  // The parameter numIters (settable with SetNumIters) is the max number of function calls.
999  // It can be reduced to make teh integration faster but it will be difficult to reach the required tolerance
1000  // VEGAS MC integration method based on importance sampling - numIters is number of function calls
1001  // Extra Vegas parameter can be set using IntegratorMultiDimOptions class
1002  // MISER MC integration method based on stratified sampling
1003  // See also http://en.wikipedia.org/wiki/Monte_Carlo_integration for VEGAS and MISER description
1004  // PLAIN simple MC integration method, where the max number of calls can be specified using SetNumIters(numIters)
1005  //
1006  // Extra integration types are:
1007  //
1008  // TOYMC:
1009  // evaluate posterior by generating toy MC for the nuisance parameters. It is a MC
1010  // integration, where the function is sampled according to the nuisance. It is convenient to use when all
1011  // the nuisance are uncorrelated and it is efficient to generate them
1012  // The toy are generated by default for each poi values
1013  // (this method has been proposed and provided by J.P Chou)
1014  // 1-TOYMC : same method as before but in this case the toys are generated only one time and then used for
1015  // each poi value. It can be convenient when the generation time is much larger than the evaluation time,
1016  // otherwise it is recoomended to re-generate the toy for each poi scanned point of the posterior function
1017  //
1018  //
1019  // ROOFIT:
1020  // use roofit default integration methods which will produce a nested integral (not reccomended for more
1021  // than 1 nuisance parameters)
1022  //
1023  // if type = 0 use default specified via class IntegratorMultiDimOptions::SetDefaultIntegrator
1024  fIntegrationType = TString(type);
1026 }
1027 
1028 
1029 
1031 {
1032  // Compute the interval. By Default a central interval is computed
1033  // and the result is a SimpleInterval object.
1034  // Using the method (to be called before SetInterval) SetLeftSideTailFraction the user can choose the type of interval.
1035  // By default the returned interval is a central interval with the confidence level specified
1036  // previously in the constructor ( LeftSideTailFraction = 0.5).
1037  // For lower limit use SetLeftSideTailFraction = 1
1038  // For upper limit use SetLeftSideTailFraction = 0
1039  // for shortest intervals use SetLeftSideTailFraction = -1 or call the method SetShortestInterval()
1040  // NOTE: The BayesianCalculator covers only the case with one
1041  // single parameter of interest
1042  // NOTE: User takes ownership of the returned object
1043 
1044  if (fValidInterval)
1045  coutW(Eval) << "BayesianCalculator::GetInterval - recomputing interval for the same CL and same model" << std::endl;
1046 
1047  RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() );
1048  if (!poi) {
1049  coutE(Eval) << "BayesianCalculator::GetInterval - no parameter of interest is set " << std::endl;
1050  return 0;
1051  }
1052 
1053 
1054 
1055  // get integrated likelihood (posterior function)
1057 
1058  //Bool_t silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
1060 
1061  if (fLeftSideFraction < 0 ) {
1062  // compute short intervals
1064  }
1065  else {
1066  // compute the other intervals
1067 
1068  double lowerCutOff = fLeftSideFraction * fSize;
1069  double upperCutOff = 1. - (1.- fLeftSideFraction) * fSize;
1070 
1071 
1072  if (fNScanBins > 0) {
1073  ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
1074  }
1075 
1076  else {
1077  // use integration method if there are nuisance parameters
1078  if (fNuisanceParameters.getSize() > 0) {
1079  ComputeIntervalFromCdf(lowerCutOff, upperCutOff);
1080  }
1081  else {
1082  // case of no nuisance - just use createCdf from roofit
1083  ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff);
1084  }
1085  // case cdf failed (scan then the posterior)
1086  if (!fValidInterval) {
1087  fNScanBins = 100;
1088  coutW(Eval) << "BayesianCalculator::GetInterval - computing integral from cdf failed - do a scan in "
1089  << fNScanBins << " nbins " << std::endl;
1090  ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
1091  }
1092  }
1093  }
1094 
1095 
1096  // reset the counts and default mode
1097  if (RooAbsReal::numEvalErrors() > 0)
1098  coutW(Eval) << "BayesianCalculator::GetInterval : " << RooAbsReal::numEvalErrors() << " errors reported in evaluating log-likelihood function "
1099  << std::endl;
1100 
1103 
1104  if (!fValidInterval) {
1105  fLower = 1; fUpper = 0;
1106  coutE(Eval) << "BayesianCalculator::GetInterval - cannot compute a valid interval - return a dummy [1,0] interval"
1107  << std::endl;
1108  }
1109  else {
1110  coutI(Eval) << "BayesianCalculator::GetInterval - found a valid interval : [" << fLower << " , "
1111  << fUpper << " ]" << std::endl;
1112  }
1113 
1114  TString interval_name = TString("BayesianInterval_a") + TString(this->GetName());
1115  SimpleInterval * interval = new SimpleInterval(interval_name,*poi,fLower,fUpper,ConfidenceLevel());
1116  interval->SetTitle("SimpleInterval from BayesianCalculator");
1117 
1118  return interval;
1119 }
1120 
1122  /// Returns the value of the parameter for the point in
1123  /// parameter-space that is the most likely.
1124  /// How do we do if there are points that are equi-probable?
1125  /// use approximate posterior
1126  /// t.b.d use real function to find the mode
1127 
1130  return h->GetBinCenter(h->GetMaximumBin() );
1131  //return fApproxPosterior->GetMaximumX();
1132 }
1133 
1134 void BayesianCalculator::ComputeIntervalUsingRooFit(double lowerCutOff, double upperCutOff ) const {
1135  // internal function compute the interval using RooFit
1136 
1137  coutI(Eval) << "BayesianCalculator: Compute interval using RooFit: posteriorPdf + createCdf + RooBrentRootFinder " << std::endl;
1138 
1139  RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() );
1140  assert(poi);
1141 
1142  fValidInterval = false;
1144  if (!fPosteriorPdf) return;
1145 
1147  if (!cdf) return;
1148 
1149  RooAbsFunc* cdf_bind = cdf->bindVars(fPOI,&fPOI);
1150  if (!cdf_bind) return;
1151 
1152  RooBrentRootFinder brf(*cdf_bind);
1153  brf.setTol(fBrfPrecision); // set the brf precision
1154 
1155  double tmpVal = poi->getVal(); // patch used because findRoot changes the value of poi
1156 
1157  bool ret = true;
1158  if (lowerCutOff > 0) {
1159  double y = lowerCutOff;
1160  ret &= brf.findRoot(fLower,poi->getMin(),poi->getMax(),y);
1161  }
1162  else
1163  fLower = poi->getMin();
1164 
1165  if (upperCutOff < 1.0) {
1166  double y=upperCutOff;
1167  ret &= brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y);
1168  }
1169  else
1170  fUpper = poi->getMax();
1171  if (!ret) coutE(Eval) << "BayesianCalculator::GetInterval "
1172  << "Error returned from Root finder, estimated interval is not fully correct"
1173  << std::endl;
1174  else
1175  fValidInterval = true;
1176 
1177 
1178  poi->setVal(tmpVal); // patch: restore the original value of poi
1179 
1180  delete cdf_bind;
1181  delete cdf;
1182 }
1183 
1184 void BayesianCalculator::ComputeIntervalFromCdf(double lowerCutOff, double upperCutOff ) const {
1185  // internal function compute the interval using Cdf integration
1186 
1187  fValidInterval = false;
1188 
1189  coutI(InputArguments) << "BayesianCalculator:GetInterval Compute the interval from the posterior cdf " << std::endl;
1190 
1191  RooRealVar* poi = dynamic_cast<RooRealVar*>( fPOI.first() );
1192  assert(poi);
1193  if (GetPosteriorFunction() == 0) {
1194  coutE(InputArguments) << "BayesianCalculator::GetInterval() cannot make posterior Function " << std::endl;
1195  return;
1196  }
1197 
1198  // need to remove the constant parameters
1199  RooArgList bindParams;
1200  bindParams.add(fPOI);
1201  bindParams.add(fNuisanceParameters);
1202 
1203  // this code could be put inside the PosteriorCdfFunction
1204 
1205  //bindParams.Print("V");
1206 
1207  PosteriorCdfFunction cdf(*fLogLike, bindParams, fPriorPdf, fIntegrationType, fNLLMin );
1208  if( cdf.HasError() ) {
1209  coutE(Eval) << "BayesianCalculator: Numerical error computing CDF integral - try a different method " << std::endl;
1210  return;
1211  }
1212 
1213  //find the roots
1214 
1215  ROOT::Math::RootFinder rf(kRootFinderType);
1216 
1217  ccoutD(Eval) << "BayesianCalculator::GetInterval - finding roots of posterior using RF " << rf.Name()
1218  << " with precision = " << fBrfPrecision;
1219 
1220  if (lowerCutOff > 0) {
1221  cdf.SetOffset(lowerCutOff);
1222  ccoutD(NumIntegration) << "Integrating posterior to get cdf and search lower limit at p =" << lowerCutOff << std::endl;
1223  bool ok = rf.Solve(cdf, poi->getMin(),poi->getMax() , 200,fBrfPrecision, fBrfPrecision);
1224  if( cdf.HasError() )
1225  coutW(Eval) << "BayesianCalculator: Numerical error integrating the CDF " << std::endl;
1226  if (!ok) {
1227  coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching lower limit !" << std::endl;
1228  return;
1229  }
1230  fLower = rf.Root();
1231  }
1232  else {
1233  fLower = poi->getMin();
1234  }
1235  if (upperCutOff < 1.0) {
1236  cdf.SetOffset(upperCutOff);
1237  ccoutD(NumIntegration) << "Integrating posterior to get cdf and search upper interval limit at p =" << upperCutOff << std::endl;
1238  bool ok = rf.Solve(cdf, fLower,poi->getMax() , 200, fBrfPrecision, fBrfPrecision);
1239  if( cdf.HasError() )
1240  coutW(Eval) << "BayesianCalculator: Numerical error integrating the CDF " << std::endl;
1241  if (!ok) {
1242  coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching upper limit !" << std::endl;
1243  return;
1244  }
1245  fUpper = rf.Root();
1246  }
1247  else {
1248  fUpper = poi->getMax();
1249  }
1250 
1251  fValidInterval = true;
1252 }
1253 
1255  // approximate posterior in nbins using a TF1
1256  // scan the poi values and evaluate the posterior at each point
1257  // and save the result in a cloned TF1
1258  // For each point the posterior is evaluated by integrating the nuisance
1259  // parameters
1260 
1261  if (fApproxPosterior) {
1262  // if number of bins of existing function is >= requested one - no need to redo the scan
1263  if (fApproxPosterior->GetNpx() >= fNScanBins) return;
1264  // otherwise redo the scan
1265  delete fApproxPosterior;
1266  fApproxPosterior = 0;
1267  }
1268 
1269 
1270  RooAbsReal * posterior = GetPosteriorFunction();
1271  if (!posterior) return;
1272 
1273 
1274  TF1 * tmp = posterior->asTF(fPOI);
1275  assert(tmp != 0);
1276  // binned the function in nbins and evaluate at thos points
1277  if (fNScanBins > 0) tmp->SetNpx(fNScanBins); // if not use default of TF1 (which is 100)
1278 
1279  coutI(Eval) << "BayesianCalculator - scan posterior function in nbins = " << tmp->GetNpx() << std::endl;
1280 
1281  fApproxPosterior = (TF1*) tmp->Clone();
1282  // save this function for future reuse
1283  // I can delete now original posterior and use this approximated copy
1284  delete tmp;
1285  TString name = posterior->GetName() + TString("_approx");
1286  TString title = posterior->GetTitle() + TString("_approx");
1287  RooAbsReal * posterior2 = new RooTFnBinding(name,title,fApproxPosterior,fPOI);
1288  if (posterior == fIntegratedLikelihood) {
1289  delete fIntegratedLikelihood;
1290  fIntegratedLikelihood = posterior2;
1291  }
1292  else if (posterior == fLikelihood) {
1293  delete fLikelihood;
1294  fLikelihood = posterior2;
1295  }
1296  else {
1297  assert(1); // should never happen this case
1298  }
1299 }
1300 
1301 void BayesianCalculator::ComputeIntervalFromApproxPosterior(double lowerCutOff, double upperCutOff ) const {
1302  // compute the interval using the approximate posterior function
1303 
1304  ccoutD(Eval) << "BayesianCalculator: Compute interval from the approximate posterior " << std::endl;
1305 
1307  if (!fApproxPosterior) return;
1308 
1309  double prob[2];
1310  double limits[2] = {0,0};
1311  prob[0] = lowerCutOff;
1312  prob[1] = upperCutOff;
1313  fApproxPosterior->GetQuantiles(2,limits,prob);
1314  fLower = limits[0];
1315  fUpper = limits[1];
1316  fValidInterval = true;
1317 }
1318 
1320  // compute the shortest interval
1321  coutI(Eval) << "BayesianCalculator - computing shortest interval with CL = " << 1.-fSize << std::endl;
1322 
1323  // compute via the approx posterior function
1325  if (!fApproxPosterior) return;
1326 
1327  TH1D * h1 = dynamic_cast<TH1D*>(fApproxPosterior->GetHistogram() );
1328  assert(h1 != 0);
1330  // get bins and sort them
1331  double * bins = h1->GetArray();
1332  int n = h1->GetSize()-2; // exclude under/overflow bins
1333  std::vector<int> index(n);
1334  TMath::Sort(n, bins, &index[0]);
1335  // find cut off as test size
1336  double sum = 0;
1337  double actualCL = 0;
1338  double upper = h1->GetXaxis()->GetXmin();
1339  double lower = h1->GetXaxis()->GetXmax();
1340  double norm = h1->GetSumOfWeights();
1341 
1342  for (int i = 0; i < n; ++i) {
1343  int idx = index[i];
1344  double p = bins[ idx] / norm;
1345  sum += p;
1346  if (sum > 1.-fSize ) {
1347  actualCL = sum - p;
1348  break;
1349  }
1350 
1351  if ( h1->GetBinLowEdge(idx) < lower)
1352  lower = h1->GetBinLowEdge(idx);
1353  if ( h1->GetXaxis()->GetBinUpEdge(idx) > upper)
1354  upper = h1->GetXaxis()->GetBinUpEdge(idx);
1355  }
1356 
1357  ccoutD(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = "
1358  << actualCL << " difference from requested is " << (actualCL-(1.-fSize))/fSize*100. << "% "
1359  << " limits are [ " << lower << " , " << " upper ] " << std::endl;
1360 
1361 
1362  if (lower < upper) {
1363  fLower = lower;
1364  fUpper = upper;
1365 
1366 
1367 
1368  if ( std::abs(actualCL-(1.-fSize)) > 0.1*(1.-fSize) )
1369  coutW(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = "
1370  << actualCL << " differs more than 10% from desired CL value - must increase nbins "
1371  << n << " to an higher value " << std::endl;
1372  }
1373  else
1374  coutE(Eval) << "BayesianCalculator::ComputeShortestInterval " << n << " bins are not sufficient " << std::endl;
1375 
1376  fValidInterval = true;
1377 
1378 }
1379 
1380 
1381 
1382 
1383 } // end namespace RooStats
1384 
User Class to find the Root of one dimensional functions.
Definition: RootFinder.h:88
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:777
RooPlot * GetPosteriorPlot(bool norm=false, double precision=0.01) const
get the plot with option to get it normalized
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
virtual Int_t GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
Compute Quantiles for density distribution of this function.
Definition: TF1.cxx:1692
double par[1]
Definition: unuranDistr.cxx:38
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg(), const RooCmdArg &arg10=RooCmdArg()) const
Plot (project) PDF on specified frame.
static IntegrationMultiDim::Type GetType(const char *name)
static function to get the enumeration from a string
Definition: Integrator.cxx:79
#define coutE(a)
Definition: RooMsgService.h:35
RooAbsPdf * GetPosteriorPdf() const
return posterior pdf (object is managed by the BayesianCalculator class)
ROOT::Math::IGenFunction * fPosteriorFunction
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
RooCmdArg DrawOption(const char *opt)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:52
#define ooccoutI(o, a)
Definition: RooMsgService.h:52
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:244
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3139
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition: TH1.cxx:7952
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return NULL if not existing)
Definition: ModelConfig.h:262
#define coutI(a)
Definition: RooMsgService.h:32
#define assert(cond)
Definition: unittest.h:542
TH1 * h
Definition: legend2.C:5
void SetIntegrationType(const char *type)
set the integration type (possible type are) :
#define oocoutI(o, a)
Definition: RooMsgService.h:45
Definition: Rtypes.h:60
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1101
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1118
static void clearEvalErrorLog()
Clear the stack of evaluation error messages.
Basic string class.
Definition: TString.h:137
virtual Double_t getMin(const char *name=0) const
void SetNCalls(unsigned int calls)
set maximum number of function calls
const Bool_t kFALSE
Definition: Rtypes.h:92
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition: RooPlot.cxx:1099
Type GetType(const std::string &Name)
Definition: Systematics.cxx:34
#define coutW(a)
Definition: RooMsgService.h:34
const TKDTreeBinning * bins
double Root() const
Return the current and latest estimate of the Root.
Definition: RootFinder.h:175
static const float upper
Definition: main.cpp:49
ClassImp(TIterator) Bool_t TIterator return false
Compare two iterator objects.
Definition: TIterator.cxx:20
#define ooccoutW(o, a)
Definition: RooMsgService.h:54
#define ccoutI(a)
Definition: RooMsgService.h:39
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:250
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to vusualize the function.
Definition: TF1.cxx:1297
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
const char * Name
Definition: TXMLSetup.cxx:67
#define ooccoutE(o, a)
Definition: RooMsgService.h:55
double GetMode() const
return the mode (most probable value of the posterior function)
virtual SimpleInterval * GetInterval() const
compute the interval.
RooAbsReal * createCdf(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Create a cumulative distribution function of this p.d.f in terms of the observables listed in iset...
Definition: RooAbsPdf.cxx:3027
RooAbsArg * first() const
RooAbsReal * GetPosteriorFunction() const
return posterior function (object is managed by the BayesianCalculator class)
double sqrt(double)
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
void ComputeIntervalUsingRooFit(double c1, double c2) const
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
Double_t x[n]
Definition: legend1.C:17
virtual Double_t GetBinLowEdge(Int_t bin) const
return bin lower edge for 1D historam Better to use h1.GetXaxis().GetBinLowEdge(bin) ...
Definition: TH1.cxx:8481
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
double cdf(double *x, double *p)
Definition: unuranDistr.cxx:44
void Print(std::ostream &os=std::cout) const
print all the options
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
TF1 * asTF(const RooArgList &obs, const RooArgList &pars=RooArgList(), const RooArgSet &nset=RooArgSet()) const
Return a ROOT TF1,2,3 object bound to this RooAbsReal with given definition of observables and parame...
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1002
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:63
RooCmdArg MoveToBack()
ClassImp(RooStats::BayesianCalculator) using namespace std
BayesianCalculator implementation.
TH1F * h1
Definition: legend1.C:5
User class for performing function minimization.
Double_t GetXmin() const
Definition: TAxis.h:137
void ComputeIntervalFromApproxPosterior(double c1, double c2) const
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:559
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:202
void ComputeIntervalFromCdf(double c1, double c2) const
virtual Double_t GetBinCenter(Int_t bin) const
return bin center for 1D historam Better to use h1.GetXaxis().GetBinCenter(bin)
Definition: TH1.cxx:8470
Numerical multi dimensional integration options.
#define ccoutE(a)
Definition: RooMsgService.h:42
RooAbsArg * find(const char *name) const
Find object with given name in list.
static Int_t numEvalErrors()
Return the number of logged evaluation errors since the last clearing.
SVector< double, 2 > v
Definition: Dict.h:5
IBaseFunctionOneDim IGenFunction
Definition: IFunctionfwd.h:24
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition: TAxis.cxx:499
static IntegrationOneDim::Type GetType(const char *name)
static function to get the enumeration from a string
Definition: Integrator.cxx:54
void SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets function to be minimized.
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
Definition: TString.cxx:2321
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:102
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7354
Double_t E()
Definition: TMath.h:54
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
void setTol(Double_t tol)
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
const Double_t * GetArray() const
Definition: TArrayD.h:45
virtual void SetModel(const ModelConfig &model)
set the model via the ModelConfig
Int_t GetSize() const
Definition: TArray.h:49
virtual void SetName(const char *name)
Change the name of this histogram.
Definition: TH1.cxx:8288
virtual Int_t GetNpx() const
Definition: TF1.h:352
TRObject operator()(const T1 &t1) const
const char * Name() const
Return the current and latest estimate of the lower value of the Root-finding interval (for bracketin...
Definition: RootFinder.h:205
RooCmdArg Precision(Double_t prec)
RooCmdArg FillColor(Color_t color)
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
int type
Definition: TGX11.cxx:120
Double_t GetXmax() const
Definition: TAxis.h:138
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:80
virtual double DoEval(double x) const =0
implementation of the evaluation function. Must be implemented by derived classes ...
#define ooccoutD(o, a)
Definition: RooMsgService.h:51
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
Namespace for new Math classes and functions.
#define name(a, b)
Definition: linkTestLib0.cpp:5
Binding & operator=(OUT(*fun)(void))
#define oocoutW(o, a)
Definition: RooMsgService.h:47
Mother of all ROOT objects.
Definition: TObject.h:58
virtual Double_t getMax(const char *name=0) const
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:73
virtual Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const
Do the root finding using the Brent-Decker method.
virtual double FValMinimum() const
Return function value at current estimate of the minimum.
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
User class for performing multidimensional integration.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
1-Dim function class
Definition: TF1.h:149
RooCmdArg VLines()
Int_t IsNaN(Double_t x)
Definition: TMath.h:617
RooCmdArg ScanNoCdf()
virtual IBaseFunctionOneDim * Clone() const =0
Clone a function.
virtual ~BayesianCalculator()
destructor
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:247
virtual bool Minimize(int maxIter, double absTol=1.E-8, double relTol=1.E-10)
Find minimum position iterating until convergence specified by the absolute and relative tolerance or...
Functor1D class for one-dimensional functions.
Definition: Functor.h:492
Int_t getSize() const
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
double exp(double)
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
Definition: ModelConfig.h:256
RooFunctor * functor(const RooArgList &obs, const RooArgList &pars=RooArgList(), const RooArgSet &nset=RooArgSet()) const
Return a RooFunctor object bound to this RooAbsReal with given definition of observables and paramete...
RooCmdArg ConditionalObservables(const RooArgSet &set)
#define ccoutD(a)
Definition: RooMsgService.h:38
virtual void SetTitle(const char *title="")
Change (i.e. set) the title of the TNamed.
Definition: TNamed.cxx:152
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
const Int_t n
Definition: legend1.C:16
bool Solve(Function &f, Derivative &d, double start, int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Definition: RootFinder.h:233
RooAbsFunc * bindVars(const RooArgSet &vars, const RooArgSet *nset=0, Bool_t clipInvalid=kFALSE) const
Create an interface adaptor f(vars) that binds us to the specified variables (in arbitrary order)...
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add element to non-owning set.
Definition: RooArgSet.cxx:448
RooCmdArg Constrain(const RooArgSet &params)
static const float lower
Definition: main.cpp:48
TAxis * GetXaxis()
Definition: TH1.h:319
void Error(ErrorHandler_t func, int code, const char *va_(fmt),...)
Write error message and call a handler, if required.
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:503