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