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