Logo ROOT  
Reference Guide
MCMCInterval.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Authors: Kevin Belasco 17/06/2009
3// Authors: Kyle Cranmer 17/06/2009
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class RooStats::MCMCInterval
13 \ingroup Roostats
14
15 MCMCInterval is a concrete implementation of the RooStats::ConfInterval
16 interface. It takes as input Markov Chain of data points in the parameter
17 space generated by Monte Carlo using the Metropolis algorithm. From the Markov
18 Chain, the confidence interval can be determined in two ways:
19
20#### Using a Kernel-Estimated PDF: (not the default method)
21
22 A RooNDKeysPdf is constructed from the data set using adaptive kernel width.
23 With this RooNDKeysPdf F, we then integrate over the most likely domain in the
24 parameter space (tallest points in the posterior RooNDKeysPdf) until the target
25 confidence level is reached within an acceptable neighborhood as defined by
26 SetEpsilon(). More specifically: we calculate the following for different
27 cutoff values C until we reach the target confidence level: \f$\int_{ F >= C } F
28 d{normset} \f$.
29 Important note: this is not the default method because of a bug in constructing
30 the RooNDKeysPdf from a weighted data set. Configure to use this method by
31 calling SetUseKeys(true), and the data set will be interpreted without weights.
32
33
34#### Using a binned data set: (the default method)
35
36 This is the binned analog of the continuous integrative method that uses the
37 kernel-estimated PDF. The points in the Markov Chain are put into a binned
38 data set and the interval is then calculated by adding the heights of the bins
39 in decreasing order until the desired level of confidence has been reached.
40 Note that this means the actual confidence level is >= the confidence level
41 prescribed by the client (unless the user calls SetHistStrict(false)). This
42 method is the default but may not remain as such in future releases, so you may
43 wish to explicitly configure to use this method by calling SetUseKeys(false)
44
45
46 These are not the only ways for the confidence interval to be determined, and
47 other possibilities are being considered being added, especially for the
48 1-dimensional case.
49
50
51 One can ask an MCMCInterval for the lower and upper limits on a specific
52 parameter of interest in the interval. Note that this works better for some
53 distributions (ones with exactly one local maximum) than others, and sometimes
54 has little value.
55*/
56
57
58#include "Rtypes.h"
59
60#include "TMath.h"
61
64#include "RooStats/Heaviside.h"
65#include "RooDataHist.h"
66#include "RooNDKeysPdf.h"
67#include "RooProduct.h"
69#include "RooRealVar.h"
70#include "RooArgList.h"
71#include "TH1.h"
72#include "TH1F.h"
73#include "TH2F.h"
74#include "TH3F.h"
75#include "RooMsgService.h"
76#include "RooGlobalFunc.h"
77#include "TObject.h"
78#include "THnSparse.h"
79#include "RooNumber.h"
80
81#include <cstdlib>
82#include <string>
83#include <algorithm>
84
86
87using namespace RooFit;
88using namespace RooStats;
89using namespace std;
90
91static const double DEFAULT_EPSILON = 0.01;
92static const double DEFAULT_DELTA = 10e-6;
93
94////////////////////////////////////////////////////////////////////////////////
95
96MCMCInterval::MCMCInterval(const char* name)
98{
99 fConfidenceLevel = 0.0;
100 fHistConfLevel = 0.0;
101 fKeysConfLevel = 0.0;
102 fTFConfLevel = 0.0;
103 fFull = 0.0;
104 fChain = nullptr;
105 fAxes = nullptr;
106 fDataHist = nullptr;
107 fSparseHist = nullptr;
108 fVector.clear();
109 fKeysPdf = nullptr;
110 fProduct = nullptr;
111 fHeaviside = nullptr;
112 fKeysDataHist = nullptr;
113 fCutoffVar = nullptr;
114 fHist = nullptr;
115 fNumBurnInSteps = 0;
116 fHistCutoff = -1;
117 fKeysCutoff = -1;
118 fDimension = 1;
119 fUseKeys = false;
120 fUseSparseHist = false;
121 fIsHistStrict = true;
125 fTFLower = -1.0 * RooNumber::infinity();
127 fVecWeight = 0;
128 fLeftSideTF = -1;
129}
130
131////////////////////////////////////////////////////////////////////////////////
132
134 const RooArgSet& parameters, MarkovChain& chain) : ConfInterval(name)
135{
136 fNumBurnInSteps = 0;
137 fConfidenceLevel = 0.0;
138 fHistConfLevel = 0.0;
139 fKeysConfLevel = 0.0;
140 fTFConfLevel = 0.0;
141 fFull = 0.0;
142 fAxes = nullptr;
143 fChain = &chain;
144 fDataHist = nullptr;
145 fSparseHist = nullptr;
146 fVector.clear();
147 fKeysPdf = nullptr;
148 fProduct = nullptr;
149 fHeaviside = nullptr;
150 fKeysDataHist = nullptr;
151 fCutoffVar = nullptr;
152 fHist = nullptr;
153 fHistCutoff = -1;
154 fKeysCutoff = -1;
155 fUseKeys = false;
156 fUseSparseHist = false;
157 fIsHistStrict = true;
159 SetParameters(parameters);
162 fTFLower = -1.0 * RooNumber::infinity();
164 fVecWeight = 0;
165 fLeftSideTF = -1;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169
171{
172 // destructor
173 delete[] fAxes;
174 delete fHist;
175 delete fChain;
176 // kbelasco: check here for memory management errors
177 delete fDataHist;
178 delete fSparseHist;
179 delete fKeysPdf;
180 delete fProduct;
181 delete fHeaviside;
182 delete fKeysDataHist;
183 delete fCutoffVar;
184}
185
187 CompareDataHistBins(RooDataHist* hist) : fDataHist(hist) {}
188 bool operator() (Int_t bin1 , Int_t bin2) {
189 fDataHist->get(bin1);
190 double n1 = fDataHist->weight();
191 fDataHist->get(bin2);
192 double n2 = fDataHist->weight();
193 return (n1 < n2);
194 }
196};
197
199 CompareSparseHistBins(THnSparse* hist) : fSparseHist(hist) {}
200 bool operator() (Long_t bin1, Long_t bin2) {
201 double n1 = fSparseHist->GetBinContent(bin1);
202 double n2 = fSparseHist->GetBinContent(bin2);
203 return (n1 < n2);
204 }
206};
207
210 fChain(chain), fParam(param) {}
212 double xi = fChain->Get(i)->getRealValue(fParam->GetName());
213 double xj = fChain->Get(j)->getRealValue(fParam->GetName());
214 return (xi < xj);
215 }
218};
219
220////////////////////////////////////////////////////////////////////////////////
221/// kbelasco: for this method, consider running DetermineInterval() if
222/// fKeysPdf==nullptr, fSparseHist==nullptr, fDataHist==nullptr, or fVector.empty()
223/// rather than just returning false. Though this should not be an issue
224/// because nobody should be able to get an MCMCInterval that has their interval
225/// or posterior representation nullptr/empty since they should only get this
226/// through the MCMCCalculator
227
228bool MCMCInterval::IsInInterval(const RooArgSet& point) const
229{
230 if (fIntervalType == kShortest) {
231 if (fUseKeys) {
232 if (fKeysPdf == nullptr)
233 return false;
234
235 // evaluate keyspdf at point and return whether >= cutoff
236 RooStats::SetParameters(&point, const_cast<RooArgSet *>(&fParameters));
238 } else {
239 if (fUseSparseHist) {
240 if (fSparseHist == nullptr)
241 return false;
242
243 // evaluate sparse hist at bin where point lies and return
244 // whether >= cutoff
246 const_cast<RooArgSet*>(&fParameters));
247 Long_t bin;
248 // kbelasco: consider making x static
249 std::vector<double> x(fDimension);
250 for (Int_t i = 0; i < fDimension; i++)
251 x[i] = fAxes[i]->getVal();
252 bin = fSparseHist->GetBin(x.data(), false);
253 double weight = fSparseHist->GetBinContent((Long64_t)bin);
254 return (weight >= (double)fHistCutoff);
255 } else {
256 if (fDataHist == nullptr)
257 return false;
258
259 // evaluate data hist at bin where point lies and return whether
260 // >= cutoff
261 Int_t bin;
262 bin = fDataHist->getIndex(point);
263 fDataHist->get(bin);
264 return (fDataHist->weight() >= (double)fHistCutoff);
265 }
266 }
267 } else if (fIntervalType == kTailFraction) {
268 if (fVector.empty())
269 return false;
270
271 // return whether value of point is within the range
272 double x = point.getRealValue(fAxes[0]->GetName());
273 if (fTFLower <= x && x <= fTFUpper)
274 return true;
275
276 return false;
277 }
278
279 coutE(InputArguments) << "Error in MCMCInterval::IsInInterval: "
280 << "Interval type not set. Returning false." << endl;
281 return false;
282}
283
284////////////////////////////////////////////////////////////////////////////////
285
287{
288 fConfidenceLevel = cl;
290}
291
292// kbelasco: update this or just take it out
293// kbelasco: consider keeping this around but changing the implementation
294// to set the number of bins for each RooRealVar and then recreating the
295// histograms
296//void MCMCInterval::SetNumBins(Int_t numBins)
297//{
298// if (numBins > 0) {
299// fPreferredNumBins = numBins;
300// for (Int_t d = 0; d < fDimension; d++)
301// fNumBins[d] = numBins;
302// }
303// else {
304// coutE(Eval) << "* Error in MCMCInterval::SetNumBins: " <<
305// "Negative number of bins given: " << numBins << endl;
306// return;
307// }
308//
309// // If the histogram already exists, recreate it with the new bin numbers
310// if (fHist != nullptr)
311// CreateHist();
312//}
313
314////////////////////////////////////////////////////////////////////////////////
315
317{
318 Int_t size = axes.getSize();
319 if (size != fDimension) {
320 coutE(InputArguments) << "* Error in MCMCInterval::SetAxes: " <<
321 "number of variables in axes (" << size <<
322 ") doesn't match number of parameters ("
323 << fDimension << ")" << endl;
324 return;
325 }
326 for (Int_t i = 0; i < size; i++)
327 fAxes[i] = (RooRealVar*)axes.at(i);
328}
329
330////////////////////////////////////////////////////////////////////////////////
331
333{
334 // kbelasco: check here for memory leak. does RooNDKeysPdf use
335 // the RooArgList passed to it or does it make a clone?
336 // also check for memory leak from chain, does RooNDKeysPdf clone that?
337 if (fAxes == nullptr || fParameters.empty()) {
338 coutE(InputArguments) << "Error in MCMCInterval::CreateKeysPdf: "
339 << "parameters have not been set." << endl;
340 return;
341 }
342
343 if (fNumBurnInSteps >= fChain->Size()) {
345 "MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: " <<
346 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
347 "in Markov chain." << endl;
348 delete fKeysPdf;
349 delete fCutoffVar;
350 delete fHeaviside;
351 delete fProduct;
352 fKeysPdf = nullptr;
353 fCutoffVar = nullptr;
354 fHeaviside = nullptr;
355 fProduct = nullptr;
356 return;
357 }
358
361 RooArgList* paramsList = new RooArgList();
362 for (Int_t i = 0; i < fDimension; i++)
363 paramsList->add(*fAxes[i]);
364
365 // kbelasco: check for memory leaks with chain. who owns it? does
366 // RooNDKeysPdf take ownership?
367 fKeysPdf = new RooNDKeysPdf("keysPDF", "Keys PDF", *paramsList, *chain, "a");
368 fCutoffVar = new RooRealVar("cutoff", "cutoff", 0);
369 fHeaviside = new Heaviside("heaviside", "Heaviside", *fKeysPdf, *fCutoffVar);
370 fProduct = new RooProduct("product", "Keys PDF & Heaviside Product",
372}
373
374////////////////////////////////////////////////////////////////////////////////
375
377{
378 if (fAxes == nullptr || fChain == nullptr) {
379 coutE(Eval) << "* Error in MCMCInterval::CreateHist(): " <<
380 "Crucial data member was nullptr." << endl;
381 coutE(Eval) << "Make sure to fully construct/initialize." << endl;
382 return;
383 }
384 if (fHist != nullptr) {
385 delete fHist;
386 fHist = nullptr;
387 }
388
389 if (fNumBurnInSteps >= fChain->Size()) {
391 "MCMCInterval::CreateHist: creation of histogram failed: " <<
392 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
393 "in Markov chain." << endl;
394 return;
395 }
396
397 if (fDimension == 1)
398 fHist = new TH1F("posterior", "MCMC Posterior Histogram",
399 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax());
400
401 else if (fDimension == 2)
402 fHist = new TH2F("posterior", "MCMC Posterior Histogram",
403 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax(),
404 fAxes[1]->numBins(), fAxes[1]->getMin(), fAxes[1]->getMax());
405
406 else if (fDimension == 3)
407 fHist = new TH3F("posterior", "MCMC Posterior Histogram",
408 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax(),
409 fAxes[1]->numBins(), fAxes[1]->getMin(), fAxes[1]->getMax(),
410 fAxes[2]->numBins(), fAxes[2]->getMin(), fAxes[2]->getMax());
411
412 else {
413 coutE(Eval) << "* Error in MCMCInterval::CreateHist() : " <<
414 "TH1* couldn't handle dimension: " << fDimension << endl;
415 return;
416 }
417
418 // Fill histogram
419 Int_t size = fChain->Size();
420 const RooArgSet* entry;
421 for (Int_t i = fNumBurnInSteps; i < size; i++) {
422 entry = fChain->Get(i);
423 if (fDimension == 1)
424 ((TH1F*)fHist)->Fill(entry->getRealValue(fAxes[0]->GetName()),
425 fChain->Weight());
426 else if (fDimension == 2)
427 ((TH2F*)fHist)->Fill(entry->getRealValue(fAxes[0]->GetName()),
428 entry->getRealValue(fAxes[1]->GetName()),
429 fChain->Weight());
430 else
431 ((TH3F*)fHist)->Fill(entry->getRealValue(fAxes[0]->GetName()),
432 entry->getRealValue(fAxes[1]->GetName()),
433 entry->getRealValue(fAxes[2]->GetName()),
434 fChain->Weight());
435 }
436
437 if (fDimension >= 1)
439 if (fDimension >= 2)
441 if (fDimension >= 3)
443}
444
445////////////////////////////////////////////////////////////////////////////////
446
448{
449 if (fAxes == nullptr || fChain == nullptr) {
450 coutE(InputArguments) << "* Error in MCMCInterval::CreateSparseHist(): "
451 << "Crucial data member was nullptr." << endl;
452 coutE(InputArguments) << "Make sure to fully construct/initialize."
453 << endl;
454 return;
455 }
456 if (fSparseHist != nullptr)
457 delete fSparseHist;
458
459 std::vector<double> min(fDimension);
460 std::vector<double> max(fDimension);
461 std::vector<Int_t> bins(fDimension);
462 for (Int_t i = 0; i < fDimension; i++) {
463 min[i] = fAxes[i]->getMin();
464 max[i] = fAxes[i]->getMax();
465 bins[i] = fAxes[i]->numBins();
466 }
467 fSparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
468 fDimension, bins.data(), min.data(), max.data());
469
470 // kbelasco: it appears we need to call Sumw2() just to get the
471 // histogram to keep a running total of the weight so that Getsumw doesn't
472 // just return 0
474
475 if (fNumBurnInSteps >= fChain->Size()) {
477 "MCMCInterval::CreateSparseHist: creation of histogram failed: " <<
478 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
479 "in Markov chain." << endl;
480 }
481
482 // Fill histogram
483 Int_t size = fChain->Size();
484 const RooArgSet* entry;
485 std::vector<double> x(fDimension);
486 for (Int_t i = fNumBurnInSteps; i < size; i++) {
487 entry = fChain->Get(i);
488 for (Int_t ii = 0; ii < fDimension; ii++)
489 x[ii] = entry->getRealValue(fAxes[ii]->GetName());
490 fSparseHist->Fill(x.data(), fChain->Weight());
491 }
492}
493
494////////////////////////////////////////////////////////////////////////////////
495
497{
498 if (fParameters.empty() || fChain == nullptr) {
499 coutE(Eval) << "* Error in MCMCInterval::CreateDataHist(): " <<
500 "Crucial data member was nullptr or empty." << endl;
501 coutE(Eval) << "Make sure to fully construct/initialize." << endl;
502 return;
503 }
504
505 if (fNumBurnInSteps >= fChain->Size()) {
507 "MCMCInterval::CreateDataHist: creation of histogram failed: " <<
508 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
509 "in Markov chain." << endl;
510 fDataHist = nullptr;
511 return;
512 }
513
516}
517
518////////////////////////////////////////////////////////////////////////////////
519
521{
522 fVector.clear();
523 fVecWeight = 0;
524
525 if (fChain == nullptr) {
526 coutE(InputArguments) << "* Error in MCMCInterval::CreateVector(): " <<
527 "Crucial data member (Markov chain) was nullptr." << endl;
528 coutE(InputArguments) << "Make sure to fully construct/initialize."
529 << endl;
530 return;
531 }
532
533 if (fNumBurnInSteps >= fChain->Size()) {
535 "MCMCInterval::CreateVector: creation of vector failed: " <<
536 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
537 "in Markov chain." << endl;
538 }
539
540 // Fill vector
542 fVector.resize(size);
543 Int_t i;
544 Int_t chainIndex;
545 for (i = 0; i < size; i++) {
546 chainIndex = i + fNumBurnInSteps;
547 fVector[i] = chainIndex;
548 fVecWeight += fChain->Weight(chainIndex);
549 }
550
551 stable_sort(fVector.begin(), fVector.end(),
553}
554
555////////////////////////////////////////////////////////////////////////////////
556
558{
560 fParameters.add(parameters);
562 if (fAxes != nullptr)
563 delete[] fAxes;
565 Int_t n = 0;
566 for (auto *obj : fParameters) {
567 if (dynamic_cast<RooRealVar*>(obj) != nullptr)
568 fAxes[n] = static_cast<RooRealVar*>(obj);
569 else
570 coutE(Eval) << "* Error in MCMCInterval::SetParameters: " <<
571 obj->GetName() << " not a RooRealVar*" << std::endl;
572 n++;
573 }
574}
575
576////////////////////////////////////////////////////////////////////////////////
577
579{
580 switch (fIntervalType) {
581 case kShortest:
583 break;
584 case kTailFraction:
586 break;
587 default:
588 coutE(InputArguments) << "MCMCInterval::DetermineInterval(): " <<
589 "Error: Interval type not set" << endl;
590 break;
591 }
592}
593
594////////////////////////////////////////////////////////////////////////////////
595
597{
598 if (fUseKeys)
600 else
602}
603
604////////////////////////////////////////////////////////////////////////////////
605
607{
608 if (fLeftSideTF < 0 || fLeftSideTF > 1) {
609 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval: "
610 << "Fraction must be in the range [0, 1]. "
611 << fLeftSideTF << "is not allowed." << endl;
612 return;
613 }
614
615 if (fDimension != 1) {
616 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval(): "
617 << "Error: Can only find a tail-fraction interval for 1-D intervals"
618 << endl;
619 return;
620 }
621
622 if (fAxes == nullptr) {
623 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval(): "
624 << "Crucial data member was nullptr." << endl;
625 coutE(InputArguments) << "Make sure to fully construct/initialize."
626 << endl;
627 return;
628 }
629
630 // kbelasco: fill in code here to find interval
631 //
632 // also make changes so that calling GetPosterior...() returns nullptr
633 // when fIntervalType == kTailFraction, since there really
634 // is no posterior for this type of interval determination
635 if (fVector.empty())
637
638 if (fVector.empty() || fVecWeight == 0) {
639 // if size is still 0, then creation failed.
640 // if fVecWeight == 0, then there are no entries (indicates the same
641 // error as fVector.empty() because that only happens when
642 // fNumBurnInSteps >= fChain->Size())
643 // either way, reset and return
644 fVector.clear();
645 fTFLower = -1.0 * RooNumber::infinity();
647 fTFConfLevel = 0.0;
648 fVecWeight = 0;
649 return;
650 }
651
652 RooRealVar* param = fAxes[0];
653
654 double c = fConfidenceLevel;
655 double leftTailCutoff = fVecWeight * (1 - c) * fLeftSideTF;
656 double rightTailCutoff = fVecWeight * (1 - c) * (1 - fLeftSideTF);
657 double leftTailSum = 0;
658 double rightTailSum = 0;
659
660 // kbelasco: consider changing these values to +infinity and -infinity
661 double ll = param->getMin();
662 double ul = param->getMax();
663
664 double x;
665 double w;
666
667 // save a lot of GetName() calls if compiler does not already optimize this
668 const char* name = param->GetName();
669
670 // find lower limit
671 Int_t i;
672 for (i = 0; i < (Int_t)fVector.size(); i++) {
674 w = fChain->Weight();
675
676 if (TMath::Abs(leftTailSum + w - leftTailCutoff) <
677 TMath::Abs(leftTailSum - leftTailCutoff)) {
678 // moving the lower limit to x would bring us closer to the desired
679 // left tail size
680 ll = x;
681 leftTailSum += w;
682 } else
683 break;
684 }
685
686 // find upper limit
687 for (i = (Int_t)fVector.size() - 1; i >= 0; i--) {
689 w = fChain->Weight();
690
691 if (TMath::Abs(rightTailSum + w - rightTailCutoff) <
692 TMath::Abs(rightTailSum - rightTailCutoff)) {
693 // moving the lower limit to x would bring us closer to the desired
694 // left tail size
695 ul = x;
696 rightTailSum += w;
697 } else
698 break;
699 }
700
701 fTFLower = ll;
702 fTFUpper = ul;
703 fTFConfLevel = 1 - (leftTailSum + rightTailSum) / fVecWeight;
704}
705
706////////////////////////////////////////////////////////////////////////////////
707
709{
710 if (fKeysPdf == nullptr)
712
713 if (fKeysPdf == nullptr) {
714 // if fKeysPdf is still nullptr, then it means CreateKeysPdf failed
715 // so clear all the data members this function would normally determine
716 // and return
717 fFull = 0.0;
718 fKeysCutoff = -1;
719 fKeysConfLevel = 0.0;
720 return;
721 }
722
723 // now we have a keys pdf of the posterior
724
725 double cutoff = 0.0;
726 fCutoffVar->setVal(cutoff);
728 double full = integral->getVal(fParameters);
729 fFull = full;
730 delete integral;
731 if (full < 0.98) {
732 coutW(Eval) << "Warning: Integral of Keys PDF came out to " << full
733 << " instead of expected value 1. Will continue using this "
734 << "factor to normalize further integrals of this PDF." << endl;
735 }
736
737 // kbelasco: Is there a better way to set the search range?
738 // from 0 to max value of Keys
739 // kbelasco: how to get max value?
740 //double max = product.maxVal(product.getMaxVal(fParameters));
741
742 double volume = 1.0;
743 for (auto *var : static_range_cast<RooRealVar*>(fParameters))
744 volume *= (var->getMax() - var->getMin());
745
746 double topCutoff = full / volume;
747 double bottomCutoff = topCutoff;
748 double confLevel = CalcConfLevel(topCutoff, full);
749 if (AcceptableConfLevel(confLevel)) {
750 fKeysConfLevel = confLevel;
751 fKeysCutoff = topCutoff;
752 return;
753 }
754 bool changed = false;
755 // find high end of range
756 while (confLevel > fConfidenceLevel) {
757 topCutoff *= 2.0;
758 confLevel = CalcConfLevel(topCutoff, full);
759 if (AcceptableConfLevel(confLevel)) {
760 fKeysConfLevel = confLevel;
761 fKeysCutoff = topCutoff;
762 return;
763 }
764 changed = true;
765 }
766 if (changed) {
767 bottomCutoff = topCutoff / 2.0;
768 } else {
769 changed = false;
770 bottomCutoff /= 2.0;
771 confLevel = CalcConfLevel(bottomCutoff, full);
772 if (AcceptableConfLevel(confLevel)) {
773 fKeysConfLevel = confLevel;
774 fKeysCutoff = bottomCutoff;
775 return;
776 }
777 while (confLevel < fConfidenceLevel) {
778 bottomCutoff /= 2.0;
779 confLevel = CalcConfLevel(bottomCutoff, full);
780 if (AcceptableConfLevel(confLevel)) {
781 fKeysConfLevel = confLevel;
782 fKeysCutoff = bottomCutoff;
783 return;
784 }
785 changed = true;
786 }
787 if (changed) {
788 topCutoff = bottomCutoff * 2.0;
789 }
790 }
791
792 coutI(Eval) << "range set: [" << bottomCutoff << ", " << topCutoff << "]"
793 << endl;
794
795 cutoff = (topCutoff + bottomCutoff) / 2.0;
796 confLevel = CalcConfLevel(cutoff, full);
797
798 // need to use WithinDeltaFraction() because sometimes the integrating the
799 // posterior in this binary search seems to not have enough granularity to
800 // find an acceptable conf level (small no. of strange cases).
801 // WithinDeltaFraction causes the search to terminate when
802 // topCutoff is essentially equal to bottomCutoff (compared to the magnitude
803 // of their mean).
804 while (!AcceptableConfLevel(confLevel) &&
805 !WithinDeltaFraction(topCutoff, bottomCutoff)) {
806 if (confLevel > fConfidenceLevel)
807 bottomCutoff = cutoff;
808 else if (confLevel < fConfidenceLevel)
809 topCutoff = cutoff;
810 cutoff = (topCutoff + bottomCutoff) / 2.0;
811 coutI(Eval) << "cutoff range: [" << bottomCutoff << ", "
812 << topCutoff << "]" << endl;
813 confLevel = CalcConfLevel(cutoff, full);
814 }
815
816 fKeysCutoff = cutoff;
817 fKeysConfLevel = confLevel;
818}
819
820////////////////////////////////////////////////////////////////////////////////
821
823{
824 if (fUseSparseHist)
826 else
828}
829
830////////////////////////////////////////////////////////////////////////////////
831
833{
834 Long_t numBins;
835 if (fSparseHist == nullptr)
837
838 if (fSparseHist == nullptr) {
839 // if fSparseHist is still nullptr, then CreateSparseHist failed
840 fHistCutoff = -1;
841 fHistConfLevel = 0.0;
842 return;
843 }
844
845 numBins = (Long_t)fSparseHist->GetNbins();
846
847 std::vector<Long_t> bins(numBins);
848 for (Int_t ibin = 0; ibin < numBins; ibin++)
849 bins[ibin] = (Long_t)ibin;
850 std::stable_sort(bins.begin(), bins.end(), CompareSparseHistBins(fSparseHist));
851
852 double nEntries = fSparseHist->GetSumw();
853 double sum = 0;
854 double content;
855 Int_t i;
856 // see above note on indexing to understand numBins - 3
857 for (i = numBins - 1; i >= 0; i--) {
858 content = fSparseHist->GetBinContent(bins[i]);
859 if ((sum + content) / nEntries >= fConfidenceLevel) {
860 fHistCutoff = content;
861 if (fIsHistStrict) {
862 sum += content;
863 i--;
864 break;
865 } else {
866 i++;
867 break;
868 }
869 }
870 sum += content;
871 }
872
873 if (fIsHistStrict) {
874 // keep going to find the sum
875 for ( ; i >= 0; i--) {
876 content = fSparseHist->GetBinContent(bins[i]);
877 if (content == fHistCutoff)
878 sum += content;
879 else
880 break; // content must be < fHistCutoff
881 }
882 } else {
883 // backtrack to find the cutoff and sum
884 for ( ; i < numBins; i++) {
885 content = fSparseHist->GetBinContent(bins[i]);
886 if (content > fHistCutoff) {
887 fHistCutoff = content;
888 break;
889 } else // content == fHistCutoff
890 sum -= content;
891 if (i == numBins - 1)
892 // still haven't set fHistCutoff correctly yet, and we have no bins
893 // left, so set fHistCutoff to something higher than the tallest bin
894 fHistCutoff = content + 1.0;
895 }
896 }
897
898 fHistConfLevel = sum / nEntries;
899}
900
901////////////////////////////////////////////////////////////////////////////////
902
904{
905 Int_t numBins;
906 if (fDataHist == nullptr)
908 if (fDataHist == nullptr) {
909 // if fDataHist is still nullptr, then CreateDataHist failed
910 fHistCutoff = -1;
911 fHistConfLevel = 0.0;
912 return;
913 }
914
915 numBins = fDataHist->numEntries();
916
917 std::vector<Int_t> bins(numBins);
918 for (Int_t ibin = 0; ibin < numBins; ibin++)
919 bins[ibin] = ibin;
920 std::stable_sort(bins.begin(), bins.end(), CompareDataHistBins(fDataHist));
921
922 double nEntries = fDataHist->sum(false);
923 double sum = 0;
924 double content;
925 Int_t i;
926 for (i = numBins - 1; i >= 0; i--) {
927 fDataHist->get(bins[i]);
928 content = fDataHist->weight();
929 if ((sum + content) / nEntries >= fConfidenceLevel) {
930 fHistCutoff = content;
931 if (fIsHistStrict) {
932 sum += content;
933 i--;
934 break;
935 } else {
936 i++;
937 break;
938 }
939 }
940 sum += content;
941 }
942
943 if (fIsHistStrict) {
944 // keep going to find the sum
945 for ( ; i >= 0; i--) {
946 fDataHist->get(bins[i]);
947 content = fDataHist->weight();
948 if (content == fHistCutoff)
949 sum += content;
950 else
951 break; // content must be < fHistCutoff
952 }
953 } else {
954 // backtrack to find the cutoff and sum
955 for ( ; i < numBins; i++) {
956 fDataHist->get(bins[i]);
957 content = fDataHist->weight();
958 if (content > fHistCutoff) {
959 fHistCutoff = content;
960 break;
961 } else // content == fHistCutoff
962 sum -= content;
963 if (i == numBins - 1)
964 // still haven't set fHistCutoff correctly yet, and we have no bins
965 // left, so set fHistCutoff to something higher than the tallest bin
966 fHistCutoff = content + 1.0;
967 }
968 }
969
970 fHistConfLevel = sum / nEntries;
971}
972
973////////////////////////////////////////////////////////////////////////////////
974
976{
977 if (fIntervalType == kShortest) {
978 if (fUseKeys)
979 return fKeysConfLevel;
980 else
981 return fHistConfLevel;
982 } else if (fIntervalType == kTailFraction) {
983 return fTFConfLevel;
984 } else {
985 coutE(InputArguments) << "MCMCInterval::GetActualConfidenceLevel: "
986 << "not implemented for this type of interval. Returning 0." << endl;
987 return 0;
988 }
989}
990
991////////////////////////////////////////////////////////////////////////////////
992
994{
995 switch (fIntervalType) {
996 case kShortest:
997 return LowerLimitShortest(param);
998 case kTailFraction:
999 return LowerLimitTailFraction(param);
1000 default:
1001 coutE(InputArguments) << "MCMCInterval::LowerLimit(): " <<
1002 "Error: Interval type not set" << endl;
1003 return RooNumber::infinity();
1004 }
1005}
1006
1007////////////////////////////////////////////////////////////////////////////////
1008
1010{
1011 switch (fIntervalType) {
1012 case kShortest:
1013 return UpperLimitShortest(param);
1014 case kTailFraction:
1015 return UpperLimitTailFraction(param);
1016 default:
1017 coutE(InputArguments) << "MCMCInterval::UpperLimit(): " <<
1018 "Error: Interval type not set" << endl;
1019 return RooNumber::infinity();
1020 }
1021}
1022
1023////////////////////////////////////////////////////////////////////////////////
1024
1026{
1027 if (fTFLower == -1.0 * RooNumber::infinity())
1029
1030 return fTFLower;
1031}
1032
1033////////////////////////////////////////////////////////////////////////////////
1034
1036{
1039
1040 return fTFUpper;
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044
1046{
1047 if (fUseKeys)
1048 return LowerLimitByKeys(param);
1049 else
1050 return LowerLimitByHist(param);
1051}
1052
1053////////////////////////////////////////////////////////////////////////////////
1054
1056{
1057 if (fUseKeys)
1058 return UpperLimitByKeys(param);
1059 else
1060 return UpperLimitByHist(param);
1061}
1062
1063////////////////////////////////////////////////////////////////////////////////
1064/// Determine the lower limit for param on this interval
1065/// using the binned data set
1066
1068{
1069 if (fUseSparseHist)
1070 return LowerLimitBySparseHist(param);
1071 else
1072 return LowerLimitByDataHist(param);
1073}
1074
1075////////////////////////////////////////////////////////////////////////////////
1076/// Determine the upper limit for param on this interval
1077/// using the binned data set
1078
1080{
1081 if (fUseSparseHist)
1082 return UpperLimitBySparseHist(param);
1083 else
1084 return UpperLimitByDataHist(param);
1085}
1086
1087////////////////////////////////////////////////////////////////////////////////
1088/// Determine the lower limit for param on this interval
1089/// using the binned data set
1090
1092{
1093 if (fDimension != 1) {
1094 coutE(InputArguments) << "In MCMCInterval::LowerLimitBySparseHist: "
1095 << "Sorry, will not compute lower limit unless dimension == 1" << endl;
1096 return param.getMin();
1097 }
1098 if (fHistCutoff < 0)
1099 DetermineBySparseHist(); // this initializes fSparseHist
1100
1101 if (fHistCutoff < 0) {
1102 // if fHistCutoff < 0 still, then determination of interval failed
1103 coutE(Eval) << "In MCMCInterval::LowerLimitBySparseHist: "
1104 << "couldn't determine cutoff. Check that num burn in steps < num "
1105 << "steps in the Markov chain. Returning param.getMin()." << endl;
1106 return param.getMin();
1107 }
1108
1109 std::vector<Int_t> coord(fDimension);
1110 for (Int_t d = 0; d < fDimension; d++) {
1111 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1112 Long_t numBins = (Long_t)fSparseHist->GetNbins();
1113 double lowerLimit = param.getMax();
1114 double val;
1115 for (Long_t i = 0; i < numBins; i++) {
1116 if (fSparseHist->GetBinContent(i, &coord[0]) >= fHistCutoff) {
1117 val = fSparseHist->GetAxis(d)->GetBinCenter(coord[d]);
1118 if (val < lowerLimit)
1119 lowerLimit = val;
1120 }
1121 }
1122 return lowerLimit;
1123 }
1124 }
1125 return param.getMin();
1126}
1127
1128////////////////////////////////////////////////////////////////////////////////
1129/// Determine the lower limit for param on this interval
1130/// using the binned data set
1131
1133{
1134 if (fHistCutoff < 0)
1135 DetermineByDataHist(); // this initializes fDataHist
1136
1137 if (fHistCutoff < 0) {
1138 // if fHistCutoff < 0 still, then determination of interval failed
1139 coutE(Eval) << "In MCMCInterval::LowerLimitByDataHist: "
1140 << "couldn't determine cutoff. Check that num burn in steps < num "
1141 << "steps in the Markov chain. Returning param.getMin()." << endl;
1142 return param.getMin();
1143 }
1144
1145 for (Int_t d = 0; d < fDimension; d++) {
1146 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1147 Int_t numBins = fDataHist->numEntries();
1148 double lowerLimit = param.getMax();
1149 double val;
1150 for (Int_t i = 0; i < numBins; i++) {
1151 fDataHist->get(i);
1152 if (fDataHist->weight() >= fHistCutoff) {
1153 val = fDataHist->get()->getRealValue(param.GetName());
1154 if (val < lowerLimit)
1155 lowerLimit = val;
1156 }
1157 }
1158 return lowerLimit;
1159 }
1160 }
1161 return param.getMin();
1162}
1163
1164////////////////////////////////////////////////////////////////////////////////
1165/// Determine the upper limit for param on this interval
1166/// using the binned data set
1167
1169{
1170 if (fDimension != 1) {
1171 coutE(InputArguments) << "In MCMCInterval::UpperLimitBySparseHist: "
1172 << "Sorry, will not compute upper limit unless dimension == 1" << endl;
1173 return param.getMax();
1174 }
1175 if (fHistCutoff < 0)
1176 DetermineBySparseHist(); // this initializes fSparseHist
1177
1178 if (fHistCutoff < 0) {
1179 // if fHistCutoff < 0 still, then determination of interval failed
1180 coutE(Eval) << "In MCMCInterval::UpperLimitBySparseHist: "
1181 << "couldn't determine cutoff. Check that num burn in steps < num "
1182 << "steps in the Markov chain. Returning param.getMax()." << endl;
1183 return param.getMax();
1184 }
1185
1186 std::vector<Int_t> coord(fDimension);
1187 for (Int_t d = 0; d < fDimension; d++) {
1188 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1189 Long_t numBins = (Long_t)fSparseHist->GetNbins();
1190 double upperLimit = param.getMin();
1191 double val;
1192 for (Long_t i = 0; i < numBins; i++) {
1193 if (fSparseHist->GetBinContent(i, &coord[0]) >= fHistCutoff) {
1194 val = fSparseHist->GetAxis(d)->GetBinCenter(coord[d]);
1195 if (val > upperLimit)
1196 upperLimit = val;
1197 }
1198 }
1199 return upperLimit;
1200 }
1201 }
1202 return param.getMax();
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Determine the upper limit for param on this interval
1207/// using the binned data set
1208
1210{
1211 if (fHistCutoff < 0)
1212 DetermineByDataHist(); // this initializes fDataHist
1213
1214 if (fHistCutoff < 0) {
1215 // if fHistCutoff < 0 still, then determination of interval failed
1216 coutE(Eval) << "In MCMCInterval::UpperLimitByDataHist: "
1217 << "couldn't determine cutoff. Check that num burn in steps < num "
1218 << "steps in the Markov chain. Returning param.getMax()." << endl;
1219 return param.getMax();
1220 }
1221
1222 for (Int_t d = 0; d < fDimension; d++) {
1223 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1224 Int_t numBins = fDataHist->numEntries();
1225 double upperLimit = param.getMin();
1226 double val;
1227 for (Int_t i = 0; i < numBins; i++) {
1228 fDataHist->get(i);
1229 if (fDataHist->weight() >= fHistCutoff) {
1230 val = fDataHist->get()->getRealValue(param.GetName());
1231 if (val > upperLimit)
1232 upperLimit = val;
1233 }
1234 }
1235 return upperLimit;
1236 }
1237 }
1238 return param.getMax();
1239}
1240
1241////////////////////////////////////////////////////////////////////////////////
1242/// Determine the lower limit for param on this interval
1243/// using the keys pdf
1244
1246{
1247 if (fKeysCutoff < 0)
1249
1250 if (fKeysDataHist == nullptr)
1252
1253 if (fKeysCutoff < 0 || fKeysDataHist == nullptr) {
1254 // failure in determination of cutoff and/or creation of histogram
1255 coutE(Eval) << "in MCMCInterval::LowerLimitByKeys(): "
1256 << "couldn't find lower limit, check that the number of burn in "
1257 << "steps < number of total steps in the Markov chain. Returning "
1258 << "param.getMin()" << endl;
1259 return param.getMin();
1260 }
1261
1262 for (Int_t d = 0; d < fDimension; d++) {
1263 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1264 Int_t numBins = fKeysDataHist->numEntries();
1265 double lowerLimit = param.getMax();
1266 double val;
1267 for (Int_t i = 0; i < numBins; i++) {
1268 fKeysDataHist->get(i);
1269 if (fKeysDataHist->weight() >= fKeysCutoff) {
1270 val = fKeysDataHist->get()->getRealValue(param.GetName());
1271 if (val < lowerLimit)
1272 lowerLimit = val;
1273 }
1274 }
1275 return lowerLimit;
1276 }
1277 }
1278 return param.getMin();
1279}
1280
1281////////////////////////////////////////////////////////////////////////////////
1282/// Determine the upper limit for param on this interval
1283/// using the keys pdf
1284
1286{
1287 if (fKeysCutoff < 0)
1289
1290 if (fKeysDataHist == nullptr)
1292
1293 if (fKeysCutoff < 0 || fKeysDataHist == nullptr) {
1294 // failure in determination of cutoff and/or creation of histogram
1295 coutE(Eval) << "in MCMCInterval::UpperLimitByKeys(): "
1296 << "couldn't find upper limit, check that the number of burn in "
1297 << "steps < number of total steps in the Markov chain. Returning "
1298 << "param.getMax()" << endl;
1299 return param.getMax();
1300 }
1301
1302 for (Int_t d = 0; d < fDimension; d++) {
1303 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1304 Int_t numBins = fKeysDataHist->numEntries();
1305 double upperLimit = param.getMin();
1306 double val;
1307 for (Int_t i = 0; i < numBins; i++) {
1308 fKeysDataHist->get(i);
1309 if (fKeysDataHist->weight() >= fKeysCutoff) {
1310 val = fKeysDataHist->get()->getRealValue(param.GetName());
1311 if (val > upperLimit)
1312 upperLimit = val;
1313 }
1314 }
1315 return upperLimit;
1316 }
1317 }
1318 return param.getMax();
1319}
1320
1321////////////////////////////////////////////////////////////////////////////////
1322/// Determine the approximate maximum value of the Keys PDF
1323
1325{
1326 if (fKeysCutoff < 0)
1328
1329 if (fKeysDataHist == nullptr)
1331
1332 if (fKeysDataHist == nullptr) {
1333 // failure in determination of cutoff and/or creation of histogram
1334 coutE(Eval) << "in MCMCInterval::KeysMax(): "
1335 << "couldn't find Keys max value, check that the number of burn in "
1336 << "steps < number of total steps in the Markov chain. Returning 0"
1337 << endl;
1338 return 0;
1339 }
1340
1341 Int_t numBins = fKeysDataHist->numEntries();
1342 double max = 0;
1343 double w;
1344 for (Int_t i = 0; i < numBins; i++) {
1345 fKeysDataHist->get(i);
1346 w = fKeysDataHist->weight();
1347 if (w > max)
1348 max = w;
1349 }
1350
1351 return max;
1352}
1353
1354////////////////////////////////////////////////////////////////////////////////
1355
1357{
1358 if (fHistCutoff < 0)
1360
1361 return fHistCutoff;
1362}
1363
1364////////////////////////////////////////////////////////////////////////////////
1365
1367{
1368 if (fKeysCutoff < 0)
1370
1371 // kbelasco: if fFull hasn't been set (because Keys creation failed because
1372 // fNumBurnInSteps >= fChain->Size()) then this will return infinity, which
1373 // seems ok to me since it will indicate error
1374
1375 return fKeysCutoff / fFull;
1376}
1377
1378////////////////////////////////////////////////////////////////////////////////
1379
1380double MCMCInterval::CalcConfLevel(double cutoff, double full)
1381{
1382 RooAbsReal* integral;
1383 double confLevel;
1384 fCutoffVar->setVal(cutoff);
1386 confLevel = integral->getVal(fParameters) / full;
1387 coutI(Eval) << "cutoff = " << cutoff << ", conf = " << confLevel << endl;
1388 //cout << "tmp: cutoff = " << cutoff << ", conf = " << confLevel << endl;
1389 delete integral;
1390 return confLevel;
1391}
1392
1393////////////////////////////////////////////////////////////////////////////////
1394
1396{
1397 if(fConfidenceLevel == 0)
1398 coutE(InputArguments) << "Error in MCMCInterval::GetPosteriorHist: "
1399 << "confidence level not set " << endl;
1400 if (fHist == nullptr)
1401 CreateHist();
1402
1403 if (fHist == nullptr)
1404 // if fHist is still nullptr, then CreateHist failed
1405 return nullptr;
1406
1407 return (TH1*) fHist->Clone("MCMCposterior_hist");
1408}
1409
1410////////////////////////////////////////////////////////////////////////////////
1411
1413{
1414 if (fConfidenceLevel == 0)
1415 coutE(InputArguments) << "Error in MCMCInterval::GetPosteriorKeysPdf: "
1416 << "confidence level not set " << endl;
1417 if (fKeysPdf == nullptr)
1418 CreateKeysPdf();
1419
1420 if (fKeysPdf == nullptr)
1421 // if fKeysPdf is still nullptr, then it means CreateKeysPdf failed
1422 return nullptr;
1423
1424 return (RooNDKeysPdf*) fKeysPdf->Clone("MCMCPosterior_keys");
1425}
1426
1427////////////////////////////////////////////////////////////////////////////////
1428
1430{
1431 if (fConfidenceLevel == 0)
1432 coutE(InputArguments) << "MCMCInterval::GetPosteriorKeysProduct: "
1433 << "confidence level not set " << endl;
1434 if (fProduct == nullptr) {
1435 CreateKeysPdf();
1437 }
1438
1439 if (fProduct == nullptr)
1440 // if fProduct is still nullptr, then it means CreateKeysPdf failed
1441 return nullptr;
1442
1443 return (RooProduct*) fProduct->Clone("MCMCPosterior_keysproduct");
1444}
1445
1446////////////////////////////////////////////////////////////////////////////////
1447
1449{
1450 // returns list of parameters
1451 return new RooArgSet(fParameters);
1452}
1453
1454////////////////////////////////////////////////////////////////////////////////
1455
1457{
1458 return (TMath::Abs(confLevel - fConfidenceLevel) < fEpsilon);
1459}
1460
1461////////////////////////////////////////////////////////////////////////////////
1462
1464{
1465 return (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2));
1466}
1467
1468////////////////////////////////////////////////////////////////////////////////
1469
1471{
1472 if (fAxes == nullptr)
1473 return;
1474 if (fProduct == nullptr)
1476 if (fProduct == nullptr)
1477 // if fProduct still nullptr, then creation failed
1478 return;
1479
1480 //RooAbsBinning** savedBinning = new RooAbsBinning*[fDimension];
1481 std::vector<Int_t> savedBins(fDimension);
1482 Int_t i;
1483 double numBins;
1484 RooRealVar* var;
1485
1486 // kbelasco: Note - the accuracy is only increased here if the binning for
1487 // each RooRealVar is uniform
1488
1489 // kbelasco: look into why saving the binnings and replacing them doesn't
1490 // work (replaces with 1 bin always).
1491 // Note: this code modifies the binning for the parameters (if they are
1492 // uniform) and sets them back to what they were. If the binnings are not
1493 // uniform, this code does nothing.
1494
1495 // first scan through fAxes to make sure all binnings are uniform, or else
1496 // we can't change the number of bins because there seems to be an error
1497 // when setting the binning itself rather than just the number of bins
1498 bool tempChangeBinning = true;
1499 for (i = 0; i < fDimension; i++) {
1500 if (!fAxes[i]->getBinning(nullptr, false, false).isUniform()) {
1501 tempChangeBinning = false;
1502 break;
1503 }
1504 }
1505
1506 // kbelasco: for 1 dimension this should be fine, but for more dimensions
1507 // the total number of bins in the histogram increases exponentially with
1508 // the dimension, so don't do this above 1-D for now.
1509 if (fDimension >= 2)
1510 tempChangeBinning = false;
1511
1512 if (tempChangeBinning) {
1513 // set high number of bins for high accuracy on lower/upper limit by keys
1514 for (i = 0; i < fDimension; i++) {
1515 var = fAxes[i];
1516 //savedBinning[i] = &var->getBinning("__binning_clone", false, true);
1517 savedBins[i] = var->getBinning(nullptr, false, false).numBins();
1518 numBins = (var->getMax() - var->getMin()) / fEpsilon;
1519 var->setBins((Int_t)numBins);
1520 }
1521 }
1522
1523 fKeysDataHist = new RooDataHist("_productDataHist",
1524 "Keys PDF & Heaviside Product Data Hist", fParameters);
1526
1527 if (tempChangeBinning) {
1528 // set the binning back to normal
1529 for (i = 0; i < fDimension; i++)
1530 //fAxes[i]->setBinning(*savedBinning[i], nullptr);
1531 //fAxes[i]->setBins(savedBinning[i]->numBins(), nullptr);
1532 fAxes[i]->setBins(savedBins[i], nullptr);
1533 }
1534}
1535
1536////////////////////////////////////////////////////////////////////////////////
1537
1538bool MCMCInterval::CheckParameters(const RooArgSet& parameterPoint) const
1539{
1540 // check that the parameters are correct
1541
1542 if (parameterPoint.getSize() != fParameters.getSize() ) {
1543 coutE(Eval) << "MCMCInterval: size is wrong, parameters don't match" << std::endl;
1544 return false;
1545 }
1546 if ( ! parameterPoint.equals( fParameters ) ) {
1547 coutE(Eval) << "MCMCInterval: size is ok, but parameters don't match" << std::endl;
1548 return false;
1549 }
1550 return true;
1551}
static const double DEFAULT_EPSILON
static const double DEFAULT_DELTA
#define d(i)
Definition: RSha256.hxx:102
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
Definition: RooMsgService.h:34
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
int Int_t
Definition: RtypesCore.h:45
long Long_t
Definition: RtypesCore.h:54
long long Long64_t
Definition: RtypesCore.h:80
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
char name[80]
Definition: TGX11.cxx:110
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:220
TRObject operator()(const T1 &t1) const
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:84
Int_t numBins() const
Return number of bins.
Definition: RooAbsBinning.h:36
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
bool empty() const
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
Int_t numBins(const char *rangeName=nullptr) const override
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, double scaleFactor, bool correctForBinVolume=false, bool showProgress=false) const
Fill a RooDataHist with values sampled from this function at the bin centers.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:104
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 std::liste...
Definition: RooAbsReal.cxx:523
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
double sum(bool correctForBinSize, bool inverseCorr=false) const
Return the sum of the weights of all bins in the histogram.
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:104
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:76
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
Generic N-dimensional implementation of a kernel estimation p.d.f.
Definition: RooNDKeysPdf.h:48
static double infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:48
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:29
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
Definition: RooRealVar.cxx:317
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:407
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Represents the Heaviside function.
Definition: Heaviside.h:18
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:33
virtual void CreateDataHist()
virtual void DetermineByDataHist()
virtual void DetermineShortestInterval()
virtual double GetActualConfidenceLevel()
virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
double fKeysConfLevel
the actual conf level determined by keys
Definition: MCMCInterval.h:292
virtual void CreateVector(RooRealVar *param)
RooDataHist * fDataHist
the binned Markov Chain data
Definition: MCMCInterval.h:282
virtual double LowerLimitByDataHist(RooRealVar &param)
determine lower limit using histogram
double fDelta
topCutoff (a) considered == bottomCutoff (b) iff
Definition: MCMCInterval.h:320
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
Definition: MCMCInterval.h:280
TH1 * fHist
the binned Markov Chain data
Definition: MCMCInterval.h:303
virtual double UpperLimitBySparseHist(RooRealVar &param)
determine upper limit using histogram
enum IntervalType fIntervalType
Definition: MCMCInterval.h:325
virtual double UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
RooRealVar * fCutoffVar
cutoff variable to use for integrating keys pdf
Definition: MCMCInterval.h:291
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void DetermineByKeys()
double fTFLower
lower limit of the tail-fraction interval
Definition: MCMCInterval.h:300
bool fUseSparseHist
whether to use sparse hist (vs. RooDataHist)
Definition: MCMCInterval.h:306
double fVecWeight
sum of weights of all entries in fVector
Definition: MCMCInterval.h:299
double fLeftSideTF
left side tail-fraction for interval
Definition: MCMCInterval.h:296
bool AcceptableConfLevel(double confLevel)
virtual void DetermineBySparseHist()
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
RooProduct * fProduct
the (keysPdf * heaviside) product
Definition: MCMCInterval.h:288
virtual double LowerLimitBySparseHist(RooRealVar &param)
determine lower limit using histogram
bool WithinDeltaFraction(double a, double b)
double fKeysCutoff
cutoff keys pdf value to be in interval
Definition: MCMCInterval.h:293
virtual double LowerLimitShortest(RooRealVar &param)
get the lower limit of param in the shortest confidence interval Note that this works better for some...
virtual double LowerLimitTailFraction(RooRealVar &param)
determine lower limit of the lower confidence interval
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
Int_t fDimension
number of variables
Definition: MCMCInterval.h:311
void SetConfidenceLevel(double cl) override
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
MCMCInterval(const char *name=nullptr)
default constructor
RooRealVar ** fAxes
array of pointers to RooRealVars representing the axes of the histogram fAxes[0] represents x-axis,...
Definition: MCMCInterval.h:314
virtual void DetermineTailFractionInterval()
double fTFConfLevel
the actual conf level of tail-fraction interval
Definition: MCMCInterval.h:297
double fHistConfLevel
the actual conf level determined by hist
Definition: MCMCInterval.h:284
virtual void CreateSparseHist()
Heaviside * fHeaviside
the Heaviside function
Definition: MCMCInterval.h:289
double fFull
Value of intergral of fProduct.
Definition: MCMCInterval.h:294
virtual double UpperLimitByDataHist(RooRealVar &param)
determine upper limit using histogram
virtual void DetermineInterval()
virtual double CalcConfLevel(double cutoff, double full)
virtual double LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
virtual double LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
virtual double UpperLimitShortest(RooRealVar &param)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
THnSparse * fSparseHist
the binned Markov Chain data
Definition: MCMCInterval.h:283
bool CheckParameters(const RooArgSet &point) const override
check if parameters are correct. (dummy implementation to start)
std::vector< Int_t > fVector
vector containing the Markov chain data
Definition: MCMCInterval.h:298
virtual double LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
double fTFUpper
upper limit of the tail-fraction interval
Definition: MCMCInterval.h:301
double fHistCutoff
cutoff bin size to be in interval
Definition: MCMCInterval.h:285
MarkovChain * fChain
the markov chain
Definition: MCMCInterval.h:279
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
virtual double UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
bool fIsHistStrict
whether the specified confidence level is a
Definition: MCMCInterval.h:307
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual double UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual void CreateKeysPdf()
virtual double UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
Int_t fNumBurnInSteps
number of steps to discard as burn in, starting from the first
Definition: MCMCInterval.h:312
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
Definition: MCMCInterval.h:305
RooDataHist * fKeysDataHist
data hist representing product
Definition: MCMCInterval.h:290
RooArgSet fParameters
parameters of interest for this interval
Definition: MCMCInterval.h:278
bool IsInInterval(const RooArgSet &point) const override
determine whether this point is in the confidence interval
virtual void DetermineByHist()
virtual void SetParameters(const RooArgSet &parameters)
Set the parameters of interest for this interval and change other internal data members accordingly.
virtual void CreateKeysDataHist()
RooNDKeysPdf * fKeysPdf
the kernel estimation pdf
Definition: MCMCInterval.h:287
double fEpsilon
acceptable error for Keys interval determination
Definition: MCMCInterval.h:318
virtual double GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual double GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
virtual double Weight() const
get the weight of the current (last indexed) entry
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:49
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:478
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:577
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetZaxis()
Definition: TH1.h:324
TAxis * GetXaxis()
Definition: TH1.h:322
TAxis * GetYaxis()
Definition: TH1.h:323
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition: TH1.cxx:2727
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:257
3-D histogram with a float per channel (see TH1 documentation)}
Definition: TH3.h:268
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:149
Double_t GetSumw() const
Definition: THnBase.h:211
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:130
Efficient multidimensional histogram.
Definition: THnSparse.h:36
Double_t GetBinContent(const Int_t *idx) const
Forwards to THnBase::GetBinContent() overload.
Definition: THnSparse.h:120
Long64_t GetBin(const Int_t *idx) const override
Definition: THnSparse.h:95
void Sumw2() override
Enable calculation of errors.
Definition: THnSparse.cxx:948
Long64_t GetNbins() const override
Definition: THnSparse.h:92
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
RooCmdArg SelectVars(const RooArgSet &vars)
RooCmdArg EventRange(Int_t nStart, Int_t nStop)
RooCmdArg NormSet(Args_t &&... argsOrArgSet)
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ InputArguments
Definition: RooGlobalFunc.h:63
Namespace for the RooStats classes.
Definition: Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:63
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
RooDataHist * fDataHist
CompareDataHistBins(RooDataHist *hist)
CompareSparseHistBins(THnSparse *hist)
CompareVectorIndices(MarkovChain *chain, RooRealVar *param)
MarkovChain * fChain
TArc a
Definition: textangle.C:12