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