ROOT   Reference Guide
Searching...
No Matches
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
15BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation
16of a credible interval using a Bayesian method.
17The class works only for one single parameter of interest and it integrates the likelihood function with the given prior
18probability 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.
20This calculator works then only for model with a single parameter of interest.
21The model can instead have several nuisance parameters which are integrated (marginalized) in the computation of the posterior function.
22The integration and normalization of the posterior is computed using numerical integration methods provided by ROOT.
23See the MCMCCalculator for model with multiple parameters of interest.
24
25The interface allows one to construct the class by passing the data set, probability density function for the model, the prior
26functions and then the parameter of interest to scan. The nuisance parameters can also be passed to be marginalized when
27computing the posterior. Alternatively, the class can be constructed by passing the data and the ModelConfig containing
28all the needed information (model pdf, prior pdf, parameter of interest, nuisance parameters, etc..)
29
30After configuring the calculator, one only needs to ask GetInterval(), which
31will return an SimpleInterval object. By default the extreme of the integral are obtained by inverting directly the
32cumulative posterior distribution. By using the method SetScanOfPosterior(nbins) the interval is then obtained by
33scanning the posterior function in the given number of points. The first method is in general faster but it requires an
34integration one extra dimension ( in the poi in addition to the nuisance parameters), therefore in some case it can be
35less robust.
36
37The class can also return the posterior function (method GetPosteriorFunction) or if needed the normalized
38posterior function (the posterior pdf) (method GetPosteriorPdf). A posterior plot is also obtained using
39the GetPosteriorPlot method.
40
41The class allows to use different integration methods for integrating in (marginalizing) the nuisances and in the poi. All the numerical
42integration methods of ROOT can be used via the method SetIntegrationType (see more in the documentation of
43this method).
44
45Calculator estimating a credible interval using the Bayesian procedure.
46The calculator computes given the model the posterior distribution and estimates the
47credible interval from the given function.
48*/
49
50
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
68
69#include "Math/IFunction.h"
71#include "Math/Integrator.h"
72#include "Math/RootFinder.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#include <memory>
87
88#include "RConfigure.h"
89
91
92namespace RooStats {
93
94
95// first some utility classes and functions
96
97#ifdef R__HAS_MATHMORE
99#else
101#endif
102
103
104
105
107 LikelihoodFunction(RooFunctor & f, RooFunctor * prior = nullptr, double offset = 0) :
108 fFunc(f), fPrior(prior),
110 {
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(nullptr,Eval) << "Likelihood evaluation ncalls = " << nCalls
125 << " x0 " << x[0] << " nll = " << nll+fOffset;
126 if (fPrior) ooccoutD(nullptr,Eval) << " prior(x) = " << (*fPrior)(x);
127 ooccoutD(nullptr,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(nullptr,Eval) << "LikelihoodFunction::() WARNING - Huge likelihood value found for parameters ";
135 for (int i = 0; i < fFunc.nObs(); ++i)
136 ooccoutW(nullptr,Eval) << " x[" << i << " ] = " << x[i];
137 ooccoutW(nullptr,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 = 0;
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
164
165public:
166
167 PosteriorCdfFunction(RooAbsReal & nll, RooArgList & bindParams, RooAbsReal * prior = nullptr, const char * integType = nullptr, double nllMinimum = 0) :
168 fFunctor(nll, bindParams, RooArgList() ), // functor
169 fPriorFunc(nullptr),
170 fLikelihood(fFunctor, nullptr, nllMinimum), // integral of exp(-nll) function
171 fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ), // integrator
172 fXmin(bindParams.size()), // vector of parameters (min values)
173 fXmax(bindParams.size())
174 {
175 if (prior) {
176 fPriorFunc = std::make_shared<RooFunctor>(*prior, bindParams, RooArgList());
178 }
179
180 fIntegrator.SetFunction(fLikelihood, bindParams.size() );
181
182 ooccoutD(nullptr,NumIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. "
183 << " nllMinimum is " << nllMinimum << std::endl;
184
185 std::vector<double> par(bindParams.size());
186 for (unsigned int i = 0; i < fXmin.size(); ++i) {
187 RooRealVar & var = static_cast<RooRealVar &>( bindParams[i]);
188 fXmin[i] = var.getMin();
189 fXmax[i] = var.getMax();
190 par[i] = var.getVal();
191 ooccoutD(nullptr,NumIntegration) << "PosteriorFunction::Integrate" << var.GetName()
192 << " in interval [ " << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
193 }
194
195 fIntegrator.Options().Print(ooccoutD(nullptr,NumIntegration));
196
197 // store max POI value because it will be changed when evaluating the function
198 fMaxPOI = fXmax[0];
199
200 // compute first the normalization with the poi
201 fNorm = (*this)( fMaxPOI );
202 if (fError) ooccoutE(nullptr,NumIntegration) << "PosteriorFunction::Error computing normalization - norm = " << fNorm << std::endl;
203 fHasNorm = true;
204 fNormCdfValues.insert(std::make_pair(fXmin[0], 0) );
205 fNormCdfValues.insert(std::make_pair(fXmax[0], 1.0) );
206
207 }
208
209 // copy constructor (needed for Cloning the object)
210 // need special treatment because integrator
211 // has no copy constructor
213 ROOT::Math::IGenFunction(),
214 fFunctor(rhs.fFunctor),
215 //fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
218 fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ), // integrator
219 fXmin( rhs.fXmin),
220 fXmax( rhs.fXmax),
221 fNorm( rhs.fNorm),
222 fNormErr( rhs.fNormErr),
223 fOffset(rhs.fOffset),
224 fMaxPOI(rhs.fMaxPOI),
225 fHasNorm(rhs.fHasNorm),
227 fError(rhs.fError),
229 {
231 // need special treatment for the smart pointer
232 // if (rhs.fPriorFunc.get() ) {
233 // fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*(rhs.fPriorFunc) ) );
234 // fLikelihood.SetPrior( fPriorFunc.get() );
235 // }
236 }
237
238
239 bool HasError() const { return fError; }
240
241
242 ROOT::Math::IGenFunction * Clone() const override {
243 ooccoutD(nullptr,NumIntegration) << " cloning function .........." << std::endl;
244 return new PosteriorCdfFunction(*this);
245 }
246
247 // offset value for computing the root
248 void SetOffset(double offset) { fOffset = offset; }
249
250private:
251
252 // make assignment operator private
254 return *this;
255 }
256
257 double DoEval (double x) const override {
258
259 // evaluate cdf at poi value x by integrating poi from [xmin,x] and all the nuisances
260 fXmax[0] = x;
261 if (x <= fXmin[0] ) return -fOffset;
262 // could also avoid a function evaluation at maximum
263 if (x >= fMaxPOI && fHasNorm) return 1. - fOffset; // cdf is bound to these values
264
265 // computes the integral using a previous cdf estimate
266 double normcdf0 = 0;
267 if (fHasNorm && fUseOldValues) {
268 // look in the map of the stored cdf values the closes one
269 std::map<double,double>::iterator itr = fNormCdfValues.upper_bound(x);
270 --itr; // upper bound returns a position 1 up of the value we want
271 if (itr != fNormCdfValues.end() ) {
272 fXmin[0] = itr->first;
273 normcdf0 = itr->second;
274 // ooccoutD(nullptr,NumIntegration) << "PosteriorCdfFunction: computing integral between in poi interval : "
275 // << fXmin[0] << " - " << fXmax[0] << std::endl;
276 }
277 }
278
279 fFunctor.binding().resetNumCall(); // reset number of calls for debug
280
281 double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]);
282 double error = fIntegrator.Error();
283 double normcdf = cdf/fNorm; // normalize the cdf
284
285 ooccoutD(nullptr,NumIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , "
286 << fXmax[0] << "] integral = " << cdf << " +/- " << error
287 << " norm-integ = " << normcdf << " cdf(x) = " << normcdf+normcdf0
288 << " ncalls = " << fFunctor.binding().numCall() << std::endl;
289
290 if (TMath::IsNaN(cdf) || cdf > std::numeric_limits<double>::max()) {
291 ooccoutE(nullptr,NumIntegration) << "PosteriorFunction::Error computing integral - cdf = "
292 << cdf << std::endl;
293 fError = true;
294 }
295
296 if (cdf != 0 && error / cdf > 0.2) {
297 oocoutW(nullptr, NumIntegration)
298 << "PosteriorCdfFunction: integration error is larger than 20 % x0 = " << fXmin[0] << " x = " << x
299 << " cdf(x) = " << cdf << " +/- " << error << std::endl;
300 }
301
302 if (!fHasNorm) {
303 oocoutI(nullptr,NumIntegration) << "PosteriorCdfFunction - integral of posterior = "
304 << cdf << " +/- " << error << std::endl;
305 fNormErr = error;
306 return cdf;
307 }
308
309 normcdf += normcdf0;
310
311 // store values in the map
312 if (fUseOldValues) {
313 fNormCdfValues.insert(std::make_pair(x, normcdf) );
314 }
315
316 double errnorm = sqrt( error*error + normcdf*normcdf * fNormErr * fNormErr )/fNorm;
317 if (normcdf > 1. + 3 * errnorm) {
318 oocoutW(nullptr,NumIntegration) << "PosteriorCdfFunction: normalized cdf values is larger than 1"
319 << " x = " << x << " normcdf(x) = " << normcdf << " +/- " << error/fNorm << std::endl;
320 }
321
322 return normcdf - fOffset; // apply an offset (for finding the roots)
323 }
324
325 mutable RooFunctor fFunctor; // functor binding nll
326 mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
327 LikelihoodFunction fLikelihood; // likelihood function
328 mutable ROOT::Math::IntegratorMultiDim fIntegrator; // integrator (mutable because Integral() is not const
329 mutable std::vector<double> fXmin; // min value of parameters (poi+nuis) -
330 mutable std::vector<double> fXmax; // max value of parameters (poi+nuis) - max poi changes so it is mutable
331 double fNorm = 1.0; // normalization value (computed in constructor)
332 mutable double fNormErr = 0.0; // normalization error value (computed in constructor)
333 double fOffset = 0; // offset for computing the root
334 double fMaxPOI = 0; // maximum value of POI
335 bool fHasNorm = false; // flag to control first call to the function
336 bool fUseOldValues = true; // use old cdf values
337 mutable bool fError = false; // flag to indicate if a numerical evaluation error occurred
338 mutable std::map<double,double> fNormCdfValues;
339};
340
341//__________________________________________________________________
342// Posterior Function class
343// 1-Dim function as function of the poi
344// and it integrated all the nuisance parameters
345
347
348public:
349
350
351 PosteriorFunction(RooAbsReal & nll, RooRealVar & poi, RooArgList & nuisParams, RooAbsReal * prior = nullptr, const char * integType = nullptr, double
352 norm = 1.0, double nllOffset = 0, int niter = 0) :
353 fFunctor(nll, nuisParams, RooArgList() ),
354 fPriorFunc(nullptr),
355 fLikelihood(fFunctor, nullptr, nllOffset),
356 fPoi(&poi),
357 fXmin(nuisParams.size()),
358 fXmax(nuisParams.size()),
359 fNorm(norm)
360 {
361
362 if (prior) {
363 fPriorFunc = std::make_shared<RooFunctor>(*prior, nuisParams, RooArgList());
365 }
366
367 ooccoutD(nullptr,NumIntegration) << "PosteriorFunction::Evaluate the posterior function by integrating the nuisances: " << std::endl;
368 for (unsigned int i = 0; i < fXmin.size(); ++i) {
369 RooRealVar & var = static_cast<RooRealVar &>( nuisParams[i]);
370 fXmin[i] = var.getMin();
371 fXmax[i] = var.getMax();
372 ooccoutD(nullptr,NumIntegration) << "PosteriorFunction::Integrate " << var.GetName()
373 << " in interval [" << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
374 }
375 if (fXmin.size() == 1) { // 1D case
376 fIntegratorOneDim = std::make_unique<ROOT::Math::Integrator>( ROOT::Math::IntegratorOneDim::GetType(integType) );
377
378 fIntegratorOneDim->SetFunction(fLikelihood);
379 // interested only in relative tolerance
380 //fIntegratorOneDim->SetAbsTolerance(1.E-300);
381 fIntegratorOneDim->Options().Print(ooccoutD(nullptr,NumIntegration) );
382 }
383 else if (fXmin.size() > 1) { // multiDim case
384 fIntegratorMultiDim = std::make_unique<ROOT::Math::IntegratorMultiDim>(ROOT::Math::IntegratorMultiDim::GetType(integType));
385 fIntegratorMultiDim->SetFunction(fLikelihood, fXmin.size());
387 if (niter > 0) {
388 opt.SetNCalls(niter);
389 fIntegratorMultiDim->SetOptions(opt);
390 }
391 //fIntegratorMultiDim->SetAbsTolerance(1.E-300);
392 // print the options
393 opt.Print(ooccoutD(nullptr,NumIntegration) );
394 }
395 }
396
397
398 ROOT::Math::IGenFunction * Clone() const override {
399 assert(1);
400 return nullptr; // cannot clone this function for integrator
401 }
402
403 double Error() const { return fError;}
404
405
406private:
407 double DoEval (double x) const override {
408
409 // evaluate posterior function at a poi value x by integrating all nuisance parameters
410
411 fPoi->setVal(x);
412 fFunctor.binding().resetNumCall(); // reset number of calls for debug
413
414 double f = 0;
415 double error = 0;
416 if (fXmin.size() == 1) { // 1D case
417 f = fIntegratorOneDim->Integral(fXmin[0],fXmax[0]);
418 error = fIntegratorOneDim->Error();
419 }
420 else if (fXmin.size() > 1) { // multi-dim case
421 f = fIntegratorMultiDim->Integral(&fXmin[0],&fXmax[0]);
422 error = fIntegratorMultiDim->Error();
423 } else {
424 // no integration to be done
425 f = fLikelihood(x);
426 }
427
428 // debug
429 ooccoutD(nullptr,NumIntegration) << "PosteriorFunction: POI value = "
430 << x << "\tf(x) = " << f << " +/- " << error
431 << " norm-f(x) = " << f/fNorm
432 << " ncalls = " << fFunctor.binding().numCall() << std::endl;
433
434 if (f != 0 && error / f > 0.2) {
435 ooccoutW(nullptr, NumIntegration)
436 << "PosteriorFunction::DoEval - Error from integration in " << fXmin.size() << " Dim is larger than 20 % "
437 << "x = " << x << " p(x) = " << f << " +/- " << error << std::endl;
438 }
439
440 fError = error / fNorm;
441 return f / fNorm;
442 }
443
445 mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
448 std::unique_ptr<ROOT::Math::Integrator> fIntegratorOneDim;
449 std::unique_ptr<ROOT::Math::IntegratorMultiDim> fIntegratorMultiDim;
450 std::vector<double> fXmin;
451 std::vector<double> fXmax;
452 double fNorm;
453 mutable double fError = 0;
454};
455
456////////////////////////////////////////////////////////////////////////////////
457/// Posterior function obtaining sampling toy MC for the nuisance according to their pdf
458
460
461public:
463 RooAbsReal *prior = nullptr, double nllOffset = 0, int niter = 0, bool redoToys = true)
464 : fFunctor(nll, nuisParams, RooArgList()),
465 fPriorFunc(nullptr),
466 fLikelihood(fFunctor, nullptr, nllOffset),
467 fPdf(&pdf),
468 fPoi(&poi),
469 fNuisParams(nuisParams),
470 fNumIterations(niter),
471
472 fRedoToys(redoToys)
473 {
474 if (niter == 0) fNumIterations = 100; // default value
475
476 if (prior) {
477 fPriorFunc = std::make_shared<RooFunctor>(*prior, nuisParams, RooArgList());
479 }
480
481 ooccoutI(nullptr,InputArguments) << "PosteriorFunctionFromToyMC::Evaluate the posterior function by randomizing the nuisances: niter " << fNumIterations << std::endl;
482
483 ooccoutI(nullptr,InputArguments) << "PosteriorFunctionFromToyMC::Pdf used for randomizing the nuisance is " << fPdf->GetName() << std::endl;
484 // check that pdf contains the nuisance
485 std::unique_ptr<RooArgSet> vars{fPdf->getVariables()};
486 for (std::size_t i = 0; i < fNuisParams.size(); ++i) {
487 if (!vars->find( fNuisParams[i].GetName() ) ) {
488 ooccoutW(nullptr,InputArguments) << "Nuisance parameter " << fNuisParams[i].GetName()
489 << " is not part of sampling pdf. "
490 << "they will be treated as constant " << std::endl;
491 }
492 }
493
494 if (!fRedoToys) {
495 ooccoutI(nullptr,InputArguments) << "PosteriorFunctionFromToyMC::Generate nuisance toys only one time (for all POI points)" << std::endl;
496 GenerateToys();
497 }
498 }
499
500 // generate first n-samples of the nuisance parameters
501 void GenerateToys() const {
502 fGenParams = std::unique_ptr<RooDataSet>{fPdf->generate(fNuisParams, fNumIterations)};
503 if(fGenParams==nullptr) {
504 ooccoutE(nullptr,InputArguments) << "PosteriorFunctionFromToyMC - failed to generate nuisance parameters" << std::endl;
505 }
506 }
507
508 double Error() const { return fError;}
509
510 ROOT::Math::IGenFunction * Clone() const override {
511 // use default copy constructor
512 //return new PosteriorFunctionFromToyMC(*this);
513 // clone not implemented
514 assert(1);
515 return nullptr;
516 }
517
518private:
519 // evaluate the posterior at the poi value x
520 double DoEval( double x) const override {
521
522 int npar = fNuisParams.size();
523 assert (npar > 0);
524
525
526 // generate the toys
527 if (fRedoToys) GenerateToys();
528 if (!fGenParams) return 0;
529
530 // evaluate posterior function at a poi value x by integrating all nuisance parameters
531
532 fPoi->setVal(x);
533
534 // loop over all of the generate data
535 double sum = 0;
536 double sum2 = 0;
537
538 for(int iter=0; iter<fNumIterations; ++iter) {
539
540 // get the set of generated parameters and set the nuisance parameters to the generated values
541 std::vector<double> p(npar);
542 for (int i = 0; i < npar; ++i) {
543 const RooArgSet* genset=fGenParams->get(iter);
544 RooAbsArg * arg = genset->find( fNuisParams[i].GetName() );
545 RooRealVar * var = dynamic_cast<RooRealVar*>(arg);
546 assert(var != nullptr);
547 p[i] = var->getVal();
548 (static_cast<RooRealVar &>( fNuisParams[i])).setVal(p[i]);
549 }
550
551 // evaluate now the likelihood function
552 double fval = fLikelihood( &p.front() );
553
554 // likelihood already must contained the pdf we have sampled
555 // so we must divided by it. The value must be normalized on all
556 // other parameters
558 double nuisPdfVal = fPdf->getVal(&arg);
559 fval /= nuisPdfVal;
560
561
562 if( fval > std::numeric_limits<double>::max() ) {
563 ooccoutE(nullptr,Eval) << "BayesianCalculator::EvalPosteriorFunctionFromToy : "
564 << "Likelihood evaluates to infinity " << std::endl;
565 ooccoutE(nullptr,Eval) << "poi value = " << x << std::endl;
566 ooccoutE(nullptr,Eval) << "Nuisance parameter values : ";
567 for (int i = 0; i < npar; ++i)
568 ooccoutE(nullptr,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
569 ooccoutE(nullptr,Eval) << " - return 0 " << std::endl;
570
571 fError = 1.E30;
572 return 0;
573 }
574 if( TMath::IsNaN(fval) ) {
575 ooccoutE(nullptr,Eval) << "BayesianCalculator::EvalPosteriorFunctionFromToy : "
576 << "Likelihood is a NaN " << std::endl;
577 ooccoutE(nullptr,Eval) << "poi value = " << x << std::endl;
578 ooccoutE(nullptr,Eval) << "Nuisance parameter values : ";
579 for (int i = 0; i < npar; ++i)
580 ooccoutE(nullptr,Eval) << fNuisParams[i].GetName() << " = " << p[i] << " ";
581 ooccoutE(nullptr,Eval) << " - return 0 " << std::endl;
582 fError = 1.E30;
583 return 0;
584 }
585
586
587
588 sum += fval;
589 sum2 += fval*fval;
590 }
591
592 // compute the average and variance
593 double val = sum/double(fNumIterations);
594 double dval2 = std::max( sum2/double(fNumIterations) - val*val, 0.0);
595 fError = std::sqrt( dval2 / fNumIterations);
596
597 // debug
598 ooccoutD(nullptr,NumIntegration) << "PosteriorFunctionFromToyMC: POI value = "
599 << x << "\tp(x) = " << val << " +/- " << fError << std::endl;
600
601
602 if (val != 0 && fError/val > 0.2 ) {
603 ooccoutW(nullptr,NumIntegration) << "PosteriorFunctionFromToyMC::DoEval"
604 << " - Error in estimating posterior is larger than 20% ! "
605 << "x = " << x << " p(x) = " << val << " +/- " << fError << std::endl;
606 }
607
608
609 return val;
610 }
611
613 mutable std::shared_ptr<RooFunctor> fPriorFunc; // functor binding the prior
615 mutable RooAbsPdf * fPdf;
618 mutable std::unique_ptr<RooDataSet> fGenParams;
620 mutable double fError = -1;
621 bool fRedoToys; // do toys every iteration
622
623};
624
625////////////////////////////////////////////////////////////////////////////////
626// Implementation of BayesianCalculator
627////////////////////////////////////////////////////////////////////////////////
628
629////////////////////////////////////////////////////////////////////////////////
630/// default constructor
631
633 fData(nullptr),
634 fPdf(nullptr),
635 fPriorPdf(nullptr),
636 fNuisancePdf(nullptr),
637 fProductPdf (nullptr), fLikelihood (nullptr), fIntegratedLikelihood (nullptr), fPosteriorPdf(nullptr),
638 fPosteriorFunction(nullptr), fApproxPosterior(nullptr),
639 fLower(0), fUpper(0),
640 fNLLMin(0),
641 fSize(0.05), fLeftSideFraction(0.5),
642 fBrfPrecision(0.00005),
643 fNScanBins(-1),
644 fNumIterations(0),
645 fValidInterval(false)
646{
647}
648
649////////////////////////////////////////////////////////////////////////////////
650/// Constructor from data set, model pdf, parameter of interests and prior pdf
651/// If nuisance parameters are given they will be integrated according either to the prior or
652/// their constraint term included in the model
653
654BayesianCalculator::BayesianCalculator( /* const char* name, const char* title, */
656 RooAbsPdf& pdf,
657 const RooArgSet& POI,
658 RooAbsPdf& priorPdf,
659 const RooArgSet* nuisanceParameters ) :
660 fData(&data),
661 fPdf(&pdf),
662 fPOI(POI),
663 fPriorPdf(&priorPdf),
664 fNuisancePdf(nullptr),
665 fProductPdf (nullptr), fLikelihood (nullptr), fIntegratedLikelihood (nullptr), fPosteriorPdf(nullptr),
666 fPosteriorFunction(nullptr), fApproxPosterior(nullptr),
667 fLower(0), fUpper(0),
668 fNLLMin(0),
669 fSize(0.05), fLeftSideFraction(0.5),
670 fBrfPrecision(0.00005),
671 fNScanBins(-1),
672 fNumIterations(0),
673 fValidInterval(false)
674{
675
677 // remove constant nuisance parameters
679}
680
681////////////////////////////////////////////////////////////////////////////////
682/// Constructor from a data set and a ModelConfig
683/// model pdf, poi and nuisances will be taken from the ModelConfig
684
686 ModelConfig & model) :
687 fData(&data),
688 fPdf(model.GetPdf()),
689 fPriorPdf( model.GetPriorPdf()),
690 fNuisancePdf(nullptr),
691 fProductPdf (nullptr), fLikelihood (nullptr), fIntegratedLikelihood (nullptr), fPosteriorPdf(nullptr),
692 fPosteriorFunction(nullptr), fApproxPosterior(nullptr),
693 fLower(0), fUpper(0),
694 fNLLMin(0),
695 fSize(0.05), fLeftSideFraction(0.5),
696 fBrfPrecision(0.00005),
697 fNScanBins(-1),
698 fNumIterations(0),
699 fValidInterval(false)
700{
701 SetModel(model);
702}
703
704
706{
707 // destructor
708 ClearAll();
709}
710
711////////////////////////////////////////////////////////////////////////////////
712/// clear all cached pdf objects
713
715 if (fProductPdf) delete fProductPdf;
716 fLogLike.reset();
717 if (fLikelihood) delete fLikelihood;
719 if (fPosteriorPdf) delete fPosteriorPdf;
722 fPosteriorPdf = nullptr;
723 fPosteriorFunction = nullptr;
724 fProductPdf = nullptr;
725 fLikelihood = nullptr;
726 fIntegratedLikelihood = nullptr;
727 fLower = 0;
728 fUpper = 0;
729 fNLLMin = 0;
730 fValidInterval = false;
731}
732
733////////////////////////////////////////////////////////////////////////////////
734/// set the model to use
735/// The model pdf, prior pdf, parameter of interest and nuisances
736/// will be taken according to the model
737
739
740 fPdf = model.GetPdf();
741 fPriorPdf = model.GetPriorPdf();
742 // assignment operator = does not do a real copy the sets (must use add method)
743 fPOI.removeAll();
747 if (model.GetParametersOfInterest()) fPOI.add( *(model.GetParametersOfInterest()) );
750 if (model.GetGlobalObservables()) fGlobalObs.add( *(model.GetGlobalObservables() ) );
751 // remove constant nuisance parameters
753
754 // invalidate the cached pointers
755 ClearAll();
756}
757
758////////////////////////////////////////////////////////////////////////////////
759/// Build and return the posterior function (not normalized) as a RooAbsReal
760/// the posterior is obtained from the product of the likelihood function and the
761/// prior pdf which is then integrated in the nuisance parameters (if existing).
762/// A prior function for the nuisance can be specified either in the prior pdf object
763/// or in the model itself. If no prior nuisance is specified, but prior parameters are then
764/// the integration is performed assuming a flat prior for the nuisance parameters.
765///
766/// NOTE: the return object is managed by the BayesianCalculator class, users do not need to delete it,
767/// but the object will be deleted when the BayesiabCalculator object is deleted
768
770{
771
773 if (fLikelihood) return fLikelihood;
774
775 // run some sanity checks
776 if (!fPdf ) {
777 coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing pdf model" << std::endl;
778 return nullptr;
779 }
780 if (fPOI.empty()) {
781 coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - missing parameter of interest" << std::endl;
782 return nullptr;
783 }
784 if (fPOI.size() > 1) {
785 coutE(InputArguments) << "BayesianCalculator::GetPosteriorPdf - current implementation works only on 1D intervals" << std::endl;
786 return nullptr;
787 }
788
789
790 std::unique_ptr<RooArgSet> constrainedParams{fPdf->getParameters(*fData)};
791 // remove the constant parameters
792 RemoveConstantParameters(&*constrainedParams);
793
794 //constrainedParams->Print("V");
795
796 // use RooFit::Constrain() to be sure constraints terms are taken into account
798
799
800
801 ccoutD(Eval) << "BayesianCalculator::GetPosteriorFunction : "
802 << " pdf value " << fPdf->getVal()
803 << " neglogLikelihood = " << fLogLike->getVal() << std::endl;
804
805 if (fPriorPdf)
806 ccoutD(Eval) << "\t\t\t priorPOI value " << fPriorPdf->getVal() << std::endl;
807
808 // check that likelihood evaluation is not infinity
809 double nllVal = fLogLike->getVal();
810 if ( nllVal > std::numeric_limits<double>::max() ) {
811 coutE(Eval) << "BayesianCalculator::GetPosteriorFunction : "
812 << " Negative log likelihood evaluates to infinity " << std::endl
813 << " Non-const Parameter values : ";
814 RooArgList p(*constrainedParams);
815 for (std::size_t i = 0; i < p.size(); ++i) {
816 RooRealVar * v = dynamic_cast<RooRealVar *>(&p[i] );
817 if (v!=nullptr) ccoutE(Eval) << v->GetName() << " = " << v->getVal() << " ";
818 }
819 ccoutE(Eval) << std::endl;
820 ccoutE(Eval) << "-- Perform a full likelihood fit of the model before or set more reasonable parameter values"
821 << std::endl;
822 coutE(Eval) << "BayesianCalculator::GetPosteriorFunction : " << " cannot compute posterior function " << std::endl;
823 return nullptr;
824 }
825
826
827
828 // need do find minimum of log-likelihood in the range to shift function
829 // to avoid numerical errors when we compute the likelihood (overflows in the exponent)
830 // N.B.: this works for only 1 parameter of interest otherwise Minuit should be used for finding the minimum
831 RooFunctor * nllFunc = fLogLike->functor(fPOI);
832 assert(nllFunc);
833 ROOT::Math::Functor1D wnllFunc(*nllFunc);
834 RooRealVar* poi = dynamic_cast<RooRealVar*>(fPOI.first() );
835 assert(poi);
836
837 // try to reduce some error messages
838 //bool silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
840
841
842 coutI(Eval) << "BayesianCalculator::GetPosteriorFunction : "
843 << " nll value " << nllVal << " poi value = " << poi->getVal() << std::endl;
844
845
847 minim.SetFunction(wnllFunc,poi->getMin(),poi->getMax() );
848 bool ret = minim.Minimize(100,1.E-3,1.E-3);
849 fNLLMin = 0;
850 if (ret) fNLLMin = minim.FValMinimum();
851
852 coutI(Eval) << "BayesianCalculator::GetPosteriorFunction : minimum of NLL vs POI for POI = "
853 << poi->getVal() << " min NLL = " << fNLLMin << std::endl;
854
855 delete nllFunc;
856
857
858 if ( fNuisanceParameters.empty() || fIntegrationType.Contains("ROOFIT") ) {
859
860 ccoutD(Eval) << "BayesianCalculator::GetPosteriorFunction : use ROOFIT integration "
861 << std::endl;
862
863#ifdef DOLATER // (not clear why this does not work)
864 // need to make in this case a likelihood from the nll and make the product with the prior
865 TString likeName = TString("likelihood_times_prior_") + TString(fPriorPdf->GetName());
866 TString formula;
867 formula.Form("exp(-@0+%f+log(@1))",fNLLMin);
868 fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike,*fPriorPdf));
869#else
870 // here use RooProdPdf (not very nice) but working
871
872 if (fLogLike) fLogLike.reset();
873 if (fProductPdf) {
874 delete fProductPdf;
875 fProductPdf = nullptr;
876 }
877
878 // // create a unique name for the product pdf
879 RooAbsPdf * pdfAndPrior = fPdf;
880 if (fPriorPdf) {
881 TString prodName = TString("product_") + TString(fPdf->GetName()) + TString("_") + TString(fPriorPdf->GetName() );
882 // save this as data member since it needs to be deleted afterwards
883 fProductPdf = new RooProdPdf(prodName,"",RooArgList(*fPdf,*fPriorPdf));
884 pdfAndPrior = fProductPdf;
885 }
886
887 std::unique_ptr<RooArgSet> constrParams{fPdf->getParameters(*fData)};
888 // remove the constant parameters
889 RemoveConstantParameters(&*constrParams);
891
892 TString likeName = TString("likelihood_times_prior_") + TString(pdfAndPrior->GetName());
893 TString formula;
894 formula.Form("exp(-@0+%f)",fNLLMin);
895 fLikelihood = new RooFormulaVar(likeName,formula,RooArgList(*fLogLike));
896#endif
897
898
899 // if no nuisance parameter we can just return the likelihood function
902 fLikelihood = nullptr;
903 } else {
904 // case of using RooFit for the integration
906 std::unique_ptr<RooAbsReal>{fLikelihood->createIntegral(fNuisanceParameters)}.release();
907 }
908
909 }
910
911 else if ( fIntegrationType.Contains("TOYMC") ) {
912 // compute the posterior as expectation values of the likelihood function
913 // sampling on the nuisance parameters
914
916
917 bool doToysEveryIteration = true;
918 // if type is 1-TOYMC or TOYMC-1
919 if ( fIntegrationType.Contains("1") || fIntegrationType.Contains("ONE") ) doToysEveryIteration = false;
920
921 RooAbsPdf * samplingPdf = (fNuisancePdf) ? fNuisancePdf : fPdf;
922 if (!fNuisancePdf) {
923 ccoutI(Eval) << "BayesianCalculator::GetPosteriorFunction : no nuisance pdf is provided, try using global pdf (this will be slower)"
924 << std::endl;
925 }
926 fPosteriorFunction = new PosteriorFunctionFromToyMC(*fLogLike, *samplingPdf, *poi, nuisParams, fPriorPdf, fNLLMin,
927 fNumIterations, doToysEveryIteration );
928
929 TString name = "toyposteriorfunction_from_";
930 name += fLogLike->GetName();
932
933 // need to scan likelihood in this case
934 if (fNScanBins <= 0) fNScanBins = 100;
935
936 }
937
938 else {
939
940 // use ROOT integration method if there are nuisance parameters
941
944
945 TString name = "posteriorfunction_from_";
946 name += fLogLike->GetName();
948
949 }
950
951 if (RooAbsReal::numEvalErrors() > 0) {
952 coutW(Eval) << "BayesianCalculator::GetPosteriorFunction : " << RooAbsReal::numEvalErrors()
953 << " errors reported in evaluating log-likelihood function " << std::endl;
954 }
957
959
960}
961
962////////////////////////////////////////////////////////////////////////////////
963/// Build and return the posterior pdf (i.e posterior function normalized to all range of poi)
964/// Note that an extra integration in the POI is required for the normalization
965/// NOTE: user must delete the returned object
966
968{
969
971 if (!plike) return nullptr;
972
973
974 // create a unique name on the posterior from the names of the components
975 TString posteriorName = this->GetName() + TString("_posteriorPdf_") + plike->GetName();
976
977 RooAbsPdf * posteriorPdf = new RooGenericPdf(posteriorName,"@0",*plike);
978
979 return posteriorPdf;
980}
981
982////////////////////////////////////////////////////////////////////////////////
983/// When am approximate posterior is computed binninig the parameter of interest (poi) range
984/// (see SetScanOfPosteriors) an histogram is created and can be returned to the user
985/// A nullptr is instead returned when the posterior is computed without binning the poi.
986///
987/// NOTE: the returned object is managed by the BayesianCalculator class,
988/// if the user wants to take ownership of the returned histogram, he needs to clone
989/// or copy the return object.
990
992{
993 return (fApproxPosterior) ? fApproxPosterior->GetHistogram() : nullptr;
994}
995
996
997////////////////////////////////////////////////////////////////////////////////
998/// return a RooPlot with the posterior and the credibility region
999/// NOTE: User takes ownership of the returned object
1000
1001RooPlot* BayesianCalculator::GetPosteriorPlot(bool norm, double precision ) const
1002{
1003
1005
1006 // if a scan is requested approximate the posterior
1007 if (fNScanBins > 0)
1009
1010 RooAbsReal * posterior = fIntegratedLikelihood;
1011 if (norm) {
1012 // delete and re-do always posterior pdf (could be invalid after approximating it)
1013 if (fPosteriorPdf) delete fPosteriorPdf;
1015 posterior = fPosteriorPdf;
1016 }
1017 if (!posterior) return nullptr;
1018
1020
1021 RooAbsRealLValue* poi = dynamic_cast<RooAbsRealLValue*>(fPOI.first() );
1022 assert(poi);
1023
1024
1025 RooPlot* plot = poi->frame();
1026 if (!plot) return nullptr;
1027
1028 // try to reduce some error messages
1030
1031 plot->SetTitle(TString("Posterior probability of parameter \"")+TString(poi->GetName())+TString("\""));
1033 posterior->plotOn(plot);
1034 plot->GetYaxis()->SetTitle("posterior function");
1035
1036 // reset the counts and default mode
1039
1040 return plot;
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044/// set the integration type (possible type are) :
1045///
1046/// - 1D integration ( used when only one nuisance and when the posterior is scanned):
1048/// - multidim:
1050/// The parameter numIters (settable with SetNumIters) is the max number of function calls.
1051/// It can be reduced to make the integration faster but it will be difficult to reach the required tolerance
1052/// - VEGAS MC integration method based on importance sampling - numIters is number of function calls
1053/// Extra Vegas parameter can be set using IntegratorMultiDimOptions class
1054/// - MISER MC integration method based on stratified sampling
1056/// - PLAIN simple MC integration method, where the max number of calls can be specified using SetNumIters(numIters)
1057///
1058/// Extra integration types are:
1059///
1060/// - TOYMC:
1061/// evaluate posterior by generating toy MC for the nuisance parameters. It is a MC
1062/// integration, where the function is sampled according to the nuisance. It is convenient to use when all
1063/// the nuisance are uncorrelated and it is efficient to generate them
1064/// The toy are generated by default for each poi values
1065/// (this method has been proposed and provided by J.P Chou)
1066/// - 1-TOYMC : same method as before but in this case the toys are generated only one time and then used for
1067/// each poi value. It can be convenient when the generation time is much larger than the evaluation time,
1068/// otherwise it is recommended to re-generate the toy for each poi scanned point of the posterior function
1069/// - ROOFIT:
1070/// use roofit default integration methods which will produce a nested integral (not recommended for more
1071/// than 1 nuisance parameters)
1072
1074 // if type = 0 use default specified via class IntegratorMultiDimOptions::SetDefaultIntegrator
1077}
1078
1079////////////////////////////////////////////////////////////////////////////////
1080/// Compute the interval. By Default a central interval is computed
1081/// and the result is a SimpleInterval object.
1082///
1083/// Using the method (to be called before SetInterval) SetLeftSideTailFraction the user can choose the type of interval.
1084/// By default the returned interval is a central interval with the confidence level specified
1085/// previously in the constructor ( LeftSideTailFraction = 0.5).
1086/// - For lower limit use SetLeftSideTailFraction = 1
1087/// - For upper limit use SetLeftSideTailFraction = 0
1088/// - for shortest intervals use SetLeftSideTailFraction = -1 or call the method SetShortestInterval()
1089///
1090/// NOTE: The BayesianCalculator covers only the case with one
1091/// single parameter of interest
1092///
1093/// NOTE: User takes ownership of the returned object
1094
1096{
1097
1098 if (fValidInterval)
1099 coutW(Eval) << "BayesianCalculator::GetInterval - recomputing interval for the same CL and same model" << std::endl;
1100
1101 RooRealVar* poi = dynamic_cast<RooRealVar*>(fPOI.first() );
1102 if (!poi) {
1103 coutE(Eval) << "BayesianCalculator::GetInterval - no parameter of interest is set " << std::endl;
1104 return nullptr;
1105 }
1106
1107
1108
1109 // get integrated likelihood (posterior function)
1111
1112 //bool silentMode = (RooMsgService::instance().globalKillBelow() >= RooFit::ERROR || RooMsgService::instance().silentMode()) ;
1114
1115 if (fLeftSideFraction < 0 ) {
1116 // compute short intervals
1118 }
1119 else {
1120 // compute the other intervals
1121
1122 double lowerCutOff = fLeftSideFraction * fSize;
1123 double upperCutOff = 1. - (1.- fLeftSideFraction) * fSize;
1124
1125
1126 if (fNScanBins > 0) {
1127 ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
1128 }
1129
1130 else {
1131 // use integration method if there are nuisance parameters
1132 if (!fNuisanceParameters.empty()) {
1133 ComputeIntervalFromCdf(lowerCutOff, upperCutOff);
1134 }
1135 else {
1136 // case of no nuisance - just use createCdf from roofit
1137 ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff);
1138 }
1139 // case cdf failed (scan then the posterior)
1140 if (!fValidInterval) {
1141 fNScanBins = 100;
1142 coutW(Eval) << "BayesianCalculator::GetInterval - computing integral from cdf failed - do a scan in "
1143 << fNScanBins << " nbins " << std::endl;
1144 ComputeIntervalFromApproxPosterior(lowerCutOff, upperCutOff);
1145 }
1146 }
1147 }
1148
1149
1150 // reset the counts and default mode
1151 if (RooAbsReal::numEvalErrors() > 0) {
1152 coutW(Eval) << "BayesianCalculator::GetInterval : " << RooAbsReal::numEvalErrors()
1153 << " errors reported in evaluating log-likelihood function " << std::endl;
1154 }
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
1194void 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
1205 std::unique_ptr<RooAbsReal> cdf{fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf())};
1206 if (!cdf) return;
1207
1208 std::unique_ptr<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) {
1231 coutE(Eval) << "BayesianCalculator::GetInterval "
1232 << "Error returned from Root finder, estimated interval is not fully correct" << std::endl;
1233 } else {
1234 fValidInterval = true;
1235 }
1236
1237 poi->setVal(tmpVal); // patch: restore the original value of poi
1238}
1239
1240////////////////////////////////////////////////////////////////////////////////
1241/// internal function compute the interval using Cdf integration
1242
1243void BayesianCalculator::ComputeIntervalFromCdf(double lowerCutOff, double upperCutOff ) const {
1244
1245 fValidInterval = false;
1246
1247 coutI(InputArguments) << "BayesianCalculator:GetInterval Compute the interval from the posterior cdf " << std::endl;
1248
1249 RooRealVar* poi = dynamic_cast<RooRealVar*>(fPOI.first() );
1250 assert(poi);
1251 if (GetPosteriorFunction() == nullptr) {
1252 coutE(InputArguments) << "BayesianCalculator::GetInterval() cannot make posterior Function " << std::endl;
1253 return;
1254 }
1255
1256 // need to remove the constant parameters
1257 RooArgList bindParams;
1260
1261 // this code could be put inside the PosteriorCdfFunction
1262
1263 //bindParams.Print("V");
1264
1266 if( cdf.HasError() ) {
1267 coutE(Eval) << "BayesianCalculator: Numerical error computing CDF integral - try a different method " << std::endl;
1268 return;
1269 }
1270
1271 //find the roots
1272
1274
1275 ccoutD(Eval) << "BayesianCalculator::GetInterval - finding roots of posterior using RF " << rf.Name()
1276 << " with precision = " << fBrfPrecision;
1277
1278 if (lowerCutOff > 0) {
1279 cdf.SetOffset(lowerCutOff);
1280 ccoutD(NumIntegration) << "Integrating posterior to get cdf and search lower limit at p =" << lowerCutOff << std::endl;
1281 bool ok = rf.Solve(cdf, poi->getMin(),poi->getMax() , 200,fBrfPrecision, fBrfPrecision);
1282 if( cdf.HasError() )
1283 coutW(Eval) << "BayesianCalculator: Numerical error integrating the CDF " << std::endl;
1284 if (!ok) {
1285 coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching lower limit !" << std::endl;
1286 return;
1287 }
1288 fLower = rf.Root();
1289 }
1290 else {
1291 fLower = poi->getMin();
1292 }
1293 if (upperCutOff < 1.0) {
1294 cdf.SetOffset(upperCutOff);
1295 ccoutD(NumIntegration) << "Integrating posterior to get cdf and search upper interval limit at p =" << upperCutOff << std::endl;
1296 bool ok = rf.Solve(cdf, fLower,poi->getMax() , 200, fBrfPrecision, fBrfPrecision);
1297 if( cdf.HasError() )
1298 coutW(Eval) << "BayesianCalculator: Numerical error integrating the CDF " << std::endl;
1299 if (!ok) {
1300 coutE(NumIntegration) << "BayesianCalculator::GetInterval - Error from root finder when searching upper limit !" << std::endl;
1301 return;
1302 }
1303 fUpper = rf.Root();
1304 }
1305 else {
1306 fUpper = poi->getMax();
1307 }
1308
1309 fValidInterval = true;
1310}
1311
1312////////////////////////////////////////////////////////////////////////////////
1313/// approximate posterior in nbins using a TF1
1314/// scan the poi values and evaluate the posterior at each point
1315/// and save the result in a cloned TF1
1316/// For each point the posterior is evaluated by integrating the nuisance
1317/// parameters
1318
1320
1321 if (fApproxPosterior) {
1322 // if number of bins of existing function is >= requested one - no need to redo the scan
1323 if (fApproxPosterior->GetNpx() >= fNScanBins) return;
1324 // otherwise redo the scan
1325 delete fApproxPosterior;
1326 fApproxPosterior = nullptr;
1327 }
1328
1329
1330 RooAbsReal * posterior = GetPosteriorFunction();
1331 if (!posterior) return;
1332
1333
1334 TF1 * tmp = posterior->asTF(fPOI);
1335 assert(tmp != nullptr);
1336 // binned the function in nbins and evaluate at those points
1337 if (fNScanBins > 0) tmp->SetNpx(fNScanBins); // if not use default of TF1 (which is 100)
1338
1339 coutI(Eval) << "BayesianCalculator - scan posterior function in nbins = " << tmp->GetNpx() << std::endl;
1340
1341 fApproxPosterior = static_cast<TF1*>(tmp->Clone());
1342 // save this function for future reuse
1343 // I can delete now original posterior and use this approximated copy
1344 delete tmp;
1345 TString name = posterior->GetName() + TString("_approx");
1346 TString title = posterior->GetTitle() + TString("_approx");
1347 RooAbsReal * posterior2 = new RooTFnBinding(name,title,fApproxPosterior,fPOI);
1348 if (posterior == fIntegratedLikelihood) {
1349 delete fIntegratedLikelihood;
1350 fIntegratedLikelihood = posterior2;
1351 }
1352 else if (posterior == fLikelihood) {
1353 delete fLikelihood;
1354 fLikelihood = posterior2;
1355 }
1356 else {
1357 assert(1); // should never happen this case
1358 }
1359}
1360
1361////////////////////////////////////////////////////////////////////////////////
1362/// compute the interval using the approximate posterior function
1363
1364void BayesianCalculator::ComputeIntervalFromApproxPosterior(double lowerCutOff, double upperCutOff ) const {
1365
1366 ccoutD(Eval) << "BayesianCalculator: Compute interval from the approximate posterior " << std::endl;
1367
1369 if (!fApproxPosterior) return;
1370
1371 double prob[2];
1372 double limits[2] = {0,0};
1373 prob[0] = lowerCutOff;
1374 prob[1] = upperCutOff;
1375 fApproxPosterior->GetQuantiles(2,limits,prob);
1376 fLower = limits[0];
1377 fUpper = limits[1];
1378 fValidInterval = true;
1379}
1380
1381////////////////////////////////////////////////////////////////////////////////
1382/// compute the shortest interval from the histogram representing the posterior
1383
1384
1386 coutI(Eval) << "BayesianCalculator - computing shortest interval with CL = " << 1.-fSize << std::endl;
1387
1388 // compute via the approx posterior function
1390 if (!fApproxPosterior) return;
1391
1392 TH1D * h1 = dynamic_cast<TH1D*>(fApproxPosterior->GetHistogram() );
1393 assert(h1 != nullptr);
1395 // get bins and sort them
1396 double * bins = h1->GetArray();
1397 // exclude under/overflow bins
1398 int n = h1->GetSize()-2;
1399 std::vector<int> index(n);
1400 // exclude bins[0] (the underflow bin content)
1401 TMath::Sort(n, bins+1, &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 // histogram bin content starts from 1
1419 if ( h1->GetBinLowEdge(idx+1) < lower)
1420 lower = h1->GetBinLowEdge(idx+1);
1421 if ( h1->GetXaxis()->GetBinUpEdge(idx+1) > upper)
1422 upper = h1->GetXaxis()->GetBinUpEdge(idx+1);
1423 }
1424
1425 ccoutD(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = "
1426 << actualCL << " difference from requested is " << (actualCL-(1.-fSize))/fSize*100. << "% "
1427 << " limits are [ " << lower << " , " << " upper ] " << std::endl;
1428
1429
1430 if (lower < upper) {
1431 fLower = lower;
1432 fUpper = upper;
1433
1434 if (std::abs(actualCL - (1. - fSize)) > 0.1 * (1. - fSize)) {
1435 coutW(Eval) << "BayesianCalculator::ComputeShortestInterval - actual interval CL = " << actualCL
1436 << " differs more than 10% from desired CL value - must increase nbins " << n
1437 << " to an higher value " << std::endl;
1438 }
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
dim_t fSize
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define ccoutE(a)
#define oocoutW(o, a)
#define coutW(a)
#define oocoutI(o, a)
#define coutE(a)
#define ooccoutW(o, a)
#define ooccoutI(o, a)
#define ooccoutD(o, a)
#define ccoutI(a)
#define ccoutD(a)
#define ooccoutE(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
@ kGray
Definition Rtypes.h:65
winID h TVirtualViewer3D TVirtualGLPainter p
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition TGX11.cxx:110
User class for performing function minimization.
void SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets function to be minimized.
bool Minimize(int maxIter, double absTol=1.E-8, double relTol=1.E-10) override
Find minimum position iterating until convergence specified by the absolute and relative tolerance or...
double FValMinimum() const override
Return function value at current estimate of the minimum.
Functor1D class for one-dimensional functions.
Definition Functor.h:95
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition IFunction.h:112
Numerical multi dimensional integration options.
void Print(std::ostream &os=std::cout) const
print all the options
void SetNCalls(unsigned int calls)
set maximum number of function calls
User class for performing multidimensional integration.
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[]
void SetFunction(Function &f, unsigned int dim)
set integration function using a generic function implementing the operator()(double *x) The dimensio...
static IntegrationMultiDim::Type GetType(const char *name)
static function to get the enumeration from a string
ROOT::Math::IntegratorMultiDimOptions Options() const
retrieve the options
double Error() const
return integration error
static IntegrationOneDim::Type GetType(const char *name)
static function to get the enumeration from a string
User Class to find the Root of one dimensional functions.
Definition RootFinder.h:73
const char * Name() const
Return the current and latest estimate of the lower value of the Root-finding interval (for bracketin...
Definition RootFinder.h:190
bool Solve(Function &f, Derivative &d, double start, int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Solve f(x) = 0, given a derivative d.
Definition RootFinder.h:224
double Root() const
Return the current and latest estimate of the Root.
Definition RootFinder.h:160
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
Storage_t const & get() const
const char * GetName() const override
Returns name of object.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
void resetNumCall() const
Reset function call counter.
Definition RooAbsFunc.h:52
Int_t numCall() const
Return number of function calls since last reset.
Definition RooAbsFunc.h:47
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:163
RooFit::OwningPtr< 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.
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:57
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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...
static Int_t numEvalErrors()
Return the number of logged evaluation errors since the last clearing.
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
static void clearEvalErrorLog()
Clear the stack of evaluation error messages.
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const
Plot (project) PDF on specified frame.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
bool findRoot(double &result, double xlo, double xhi, double value=0) const
Do the root finding using the Brent-Decker method.
void setTol(double tol)
Set convergence tolerance parameter.
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
RooCFunction1Binding is a templated implementation of class RooAbsReal that binds generic C(++) funct...
Lightweight interface adaptor that exports a RooAbsPdf as a functor.
Definition RooFunctor.h:25
Int_t nObs() const
Definition RooFunctor.h:34
RooAbsFunc & binding()
Definition RooFunctor.h:51
Implementation of a probability density function that takes a RooArgList of servers and a C++ express...
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
ROOT::Math::IGenFunction * fPosteriorFunction
function representing the posterior
RooAbsReal * fLikelihood
internal pointer to likelihood function
double fNLLMin
minimum value of Nll
int fNumIterations
number of iterations (when using ToyMC)
RooAbsPdf * GetPosteriorPdf() const
return posterior pdf (object is managed by the user)
RooPlot * GetPosteriorPlot(bool norm=false, double precision=0.01) const
get the plot with option to get it normalized
int fNScanBins
number of bins to scan, if = -1 no scan is done (default)
void ClearAll() const
clear all cached pdf objects
void ComputeShortestInterval() const
compute the shortest interval from the histogram representing the posterior
RooAbsPdf * fNuisancePdf
nuisance pdf (needed when using nuisance sampling technique)
RooArgSet fConditionalObs
conditional observables
double fBrfPrecision
root finder precision
RooAbsReal * fIntegratedLikelihood
integrated likelihood function, i.e - unnormalized posterior function
void ApproximatePosterior() const
approximate posterior in nbins using a TF1 scan the poi values and evaluate the posterior at each poi...
double fSize
size used for getting the interval
RooArgSet fNuisanceParameters
nuisance parameters
double fLeftSideFraction
fraction of probability content on left side of interval
RooArgSet fGlobalObs
global observables
double GetMode() const
return the mode (most probable value of the posterior function)
RooAbsPdf * fPdf
model pdf (could contain the nuisance pdf as constraint term)
SimpleInterval * GetInterval() const override
compute the interval.
RooAbsPdf * fProductPdf
internal pointer to model * prior
TF1 * fApproxPosterior
TF1 representing the scanned posterior function.
RooAbsReal * GetPosteriorFunction() const
return posterior function (object is managed by the BayesianCalculator class)
void SetIntegrationType(const char *type)
set the integration type (possible type are) :
RooAbsPdf * fPosteriorPdf
normalized (on the poi) posterior pdf
double fUpper
upper interval bound
double fLower
computer lower interval bound
TH1 * GetPosteriorHistogram() const
return the approximate posterior as histogram (TH1 object). Note the object is managed by the Bayesia...
void ComputeIntervalFromApproxPosterior(double c1, double c2) const
compute the interval using the approximate posterior function
void SetModel(const ModelConfig &model) override
set the model via the ModelConfig
RooAbsPdf * fPriorPdf
prior pdf (typically for the POI)
std::unique_ptr< RooAbsReal > fLogLike
internal pointer to log likelihood function
double ConfidenceLevel() const override
Get the Confidence level for the test.
~BayesianCalculator() override
destructor
void ComputeIntervalUsingRooFit(double c1, double c2) const
internal function compute the interval using RooFit
void ComputeIntervalFromCdf(double c1, double c2) const
internal function compute the interval using Cdf integration
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return nullptr if not existing)
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return nullptr if not existing)
std::shared_ptr< RooFunctor > fPriorFunc
PosteriorCdfFunction & operator=(const PosteriorCdfFunction &)
PosteriorCdfFunction(const PosteriorCdfFunction &rhs)
ROOT::Math::IntegratorMultiDim fIntegrator
PosteriorCdfFunction(RooAbsReal &nll, RooArgList &bindParams, RooAbsReal *prior=nullptr, const char *integType=nullptr, double nllMinimum=0)
double DoEval(double x) const override
implementation of the evaluation function. Must be implemented by derived classes
std::map< double, double > fNormCdfValues
ROOT::Math::IGenFunction * Clone() const override
Clone a function.
Posterior function obtaining sampling toy MC for the nuisance according to their pdf.
std::shared_ptr< RooFunctor > fPriorFunc
std::unique_ptr< RooDataSet > fGenParams
double DoEval(double x) const override
implementation of the evaluation function. Must be implemented by derived classes
PosteriorFunctionFromToyMC(RooAbsReal &nll, RooAbsPdf &pdf, RooRealVar &poi, RooArgList &nuisParams, RooAbsReal *prior=nullptr, double nllOffset=0, int niter=0, bool redoToys=true)
ROOT::Math::IGenFunction * Clone() const override
Clone a function.
PosteriorFunction(RooAbsReal &nll, RooRealVar &poi, RooArgList &nuisParams, RooAbsReal *prior=nullptr, const char *integType=nullptr, double norm=1.0, double nllOffset=0, int niter=0)
std::shared_ptr< RooFunctor > fPriorFunc
std::unique_ptr< ROOT::Math::IntegratorMultiDim > fIntegratorMultiDim
ROOT::Math::IGenFunction * Clone() const override
Clone a function.
double DoEval(double x) const override
implementation of the evaluation function. Must be implemented by derived classes
std::unique_ptr< ROOT::Math::Integrator > fIntegratorOneDim
SimpleInterval is a concrete implementation of the ConfInterval interface.
Use TF1, TF2, TF3 functions as RooFit objects.
const Float_t * GetArray() const
Definition TArrayF.h:43
Int_t GetSize() const
Definition TArray.h:47
Double_t GetXmax() const
Definition TAxis.h:140
Double_t GetXmin() const
Definition TAxis.h:139
virtual Double_t GetBinUpEdge(Int_t bin) const
Return up edge of bin.
Definition TAxis.cxx:528
1-Dim function class
Definition TF1.h:233
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function Note that this histogram is managed ...
Definition TF1.cxx:1586
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:1994
virtual Int_t GetNpx() const
Definition TF1.h:518
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:670
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:9132
TAxis * GetXaxis()
Definition TH1.h:324
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition TH1.cxx:9143
void SetName(const char *name) override
Change the name of this histogram.
Definition TH1.cxx:8951
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition TH1.cxx:7908
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Basic string class.
Definition TString.h:139
void ToUpper()
Change string to upper case.
Definition TString.cxx:1195
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
Definition TString.cxx:2356
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
RooCmdArg ScanNoCdf()
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
RooCmdArg FillColor(Color_t color)
RooCmdArg Precision(double prec)
RooCmdArg DrawOption(const char *opt)
RooCmdArg Range(const char *rangeName, bool adjustNorm=true)
RooCmdArg MoveToBack()
RooCmdArg VLines()
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Namespace for the RooStats classes.
Definition Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
const ROOT::Math::RootFinder::EType kRootFinderType
Bool_t IsNaN(Double_t x)
Definition TMath.h:892
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 TMathBase.h:431
void SetPrior(RooFunctor *prior)
LikelihoodFunction(RooFunctor &f, RooFunctor *prior=nullptr, double offset=0)
double operator()(const double *x) const
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345