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