Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 std::endl;
90
91////////////////////////////////////////////////////////////////////////////////
92
94 : ConfInterval(name), fTFLower(-RooNumber::infinity()), fTFUpper(RooNumber::infinity())
95{
96 fVector.clear();
97}
98
99////////////////////////////////////////////////////////////////////////////////
100
101MCMCInterval::MCMCInterval(const char *name, const RooArgSet &parameters, MarkovChain &chain)
102 : ConfInterval(name), fChain(&chain), fTFLower(-RooNumber::infinity()), fTFUpper(RooNumber::infinity())
103{
104 fVector.clear();
105 SetParameters(parameters);
106}
107
108////////////////////////////////////////////////////////////////////////////////
109
111{
112 // destructor
113 delete[] fAxes;
114 delete fHist;
115 delete fChain;
116 // kbelasco: check here for memory management errors
117 delete fDataHist;
118 delete fSparseHist;
119 delete fKeysPdf;
120 delete fProduct;
121 delete fHeaviside;
122 delete fKeysDataHist;
123 delete fCutoffVar;
124}
125
127 CompareDataHistBins(RooDataHist* hist) : fDataHist(hist) {}
128 bool operator() (Int_t bin1 , Int_t bin2) {
129 fDataHist->get(bin1);
130 double n1 = fDataHist->weight();
131 fDataHist->get(bin2);
132 double n2 = fDataHist->weight();
133 return (n1 < n2);
134 }
136};
137
139 CompareSparseHistBins(THnSparse* hist) : fSparseHist(hist) {}
140 bool operator() (Long_t bin1, Long_t bin2) {
141 double n1 = fSparseHist->GetBinContent(bin1);
142 double n2 = fSparseHist->GetBinContent(bin2);
143 return (n1 < n2);
144 }
146};
147
150 fChain(chain), fParam(param) {}
152 double xi = fChain->Get(i)->getRealValue(fParam->GetName());
153 double xj = fChain->Get(j)->getRealValue(fParam->GetName());
154 return (xi < xj);
155 }
158};
159
160////////////////////////////////////////////////////////////////////////////////
161/// kbelasco: for this method, consider running DetermineInterval() if
162/// fKeysPdf==nullptr, fSparseHist==nullptr, fDataHist==nullptr, or fVector.empty()
163/// rather than just returning false. Though this should not be an issue
164/// because nobody should be able to get an MCMCInterval that has their interval
165/// or posterior representation nullptr/empty since they should only get this
166/// through the MCMCCalculator
167
168bool MCMCInterval::IsInInterval(const RooArgSet& point) const
169{
170 if (fIntervalType == kShortest) {
171 if (fUseKeys) {
172 if (fKeysPdf == nullptr)
173 return false;
174
175 // evaluate keyspdf at point and return whether >= cutoff
176 RooStats::SetParameters(&point, const_cast<RooArgSet *>(&fParameters));
178 } else {
179 if (fUseSparseHist) {
180 if (fSparseHist == nullptr)
181 return false;
182
183 // evaluate sparse hist at bin where point lies and return
184 // whether >= cutoff
186 const_cast<RooArgSet*>(&fParameters));
187 Long_t bin;
188 // kbelasco: consider making x static
189 std::vector<double> x(fDimension);
190 for (Int_t i = 0; i < fDimension; i++)
191 x[i] = fAxes[i]->getVal();
192 bin = fSparseHist->GetBin(x.data(), false);
193 double weight = fSparseHist->GetBinContent((Long64_t)bin);
194 return (weight >= (double)fHistCutoff);
195 } else {
196 if (fDataHist == nullptr)
197 return false;
198
199 // evaluate data hist at bin where point lies and return whether
200 // >= cutoff
201 Int_t bin;
202 bin = fDataHist->getIndex(point);
203 fDataHist->get(bin);
204 return (fDataHist->weight() >= (double)fHistCutoff);
205 }
206 }
207 } else if (fIntervalType == kTailFraction) {
208 if (fVector.empty())
209 return false;
210
211 // return whether value of point is within the range
212 double x = point.getRealValue(fAxes[0]->GetName());
213 if (fTFLower <= x && x <= fTFUpper)
214 return true;
215
216 return false;
217 }
218
219 coutE(InputArguments) << "Error in MCMCInterval::IsInInterval: "
220 << "Interval type not set. Returning false." << endl;
221 return false;
222}
223
224////////////////////////////////////////////////////////////////////////////////
225
227{
228 fConfidenceLevel = cl;
230}
231
232// kbelasco: update this or just take it out
233// kbelasco: consider keeping this around but changing the implementation
234// to set the number of bins for each RooRealVar and then recreating the
235// histograms
236//void MCMCInterval::SetNumBins(Int_t numBins)
237//{
238// if (numBins > 0) {
239// fPreferredNumBins = numBins;
240// for (Int_t d = 0; d < fDimension; d++)
241// fNumBins[d] = numBins;
242// }
243// else {
244// coutE(Eval) << "* Error in MCMCInterval::SetNumBins: " <<
245// "Negative number of bins given: " << numBins << endl;
246// return;
247// }
248//
249// // If the histogram already exists, recreate it with the new bin numbers
250// if (fHist != nullptr)
251// CreateHist();
252//}
253
254////////////////////////////////////////////////////////////////////////////////
255
257{
258 Int_t size = axes.size();
259 if (size != fDimension) {
260 coutE(InputArguments) << "* Error in MCMCInterval::SetAxes: " <<
261 "number of variables in axes (" << size <<
262 ") doesn't match number of parameters ("
263 << fDimension << ")" << endl;
264 return;
265 }
266 for (Int_t i = 0; i < size; i++)
267 fAxes[i] = static_cast<RooRealVar*>(axes.at(i));
268}
269
270////////////////////////////////////////////////////////////////////////////////
271
273{
274 // kbelasco: check here for memory leak. does RooNDKeysPdf use
275 // the RooArgList passed to it or does it make a clone?
276 // also check for memory leak from chain, does RooNDKeysPdf clone that?
277 if (fAxes == nullptr || fParameters.empty()) {
278 coutE(InputArguments) << "Error in MCMCInterval::CreateKeysPdf: "
279 << "parameters have not been set." << endl;
280 return;
281 }
282
283 if (fNumBurnInSteps >= fChain->Size()) {
285 "MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: " <<
286 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
287 "in Markov chain." << endl;
288 delete fKeysPdf;
289 delete fCutoffVar;
290 delete fHeaviside;
291 delete fProduct;
292 fKeysPdf = nullptr;
293 fCutoffVar = nullptr;
294 fHeaviside = nullptr;
295 fProduct = nullptr;
296 return;
297 }
298
299 std::unique_ptr<RooDataSet> chain{fChain->GetAsDataSet(SelectVars(fParameters),
301 RooArgList* paramsList = new RooArgList();
302 for (Int_t i = 0; i < fDimension; i++)
303 paramsList->add(*fAxes[i]);
304
305 fKeysPdf = new RooNDKeysPdf("keysPDF", "Keys PDF", *paramsList, *chain, "a");
306 fCutoffVar = new RooRealVar("cutoff", "cutoff", 0);
307 fHeaviside = new Heaviside("heaviside", "Heaviside", *fKeysPdf, *fCutoffVar);
308 fProduct = new RooProduct("product", "Keys PDF & Heaviside Product",
310}
311
312////////////////////////////////////////////////////////////////////////////////
313
315{
316 if (fAxes == nullptr || fChain == nullptr) {
317 coutE(Eval) << "* Error in MCMCInterval::CreateHist(): " <<
318 "Crucial data member was nullptr." << endl;
319 coutE(Eval) << "Make sure to fully construct/initialize." << endl;
320 return;
321 }
322 if (fHist != nullptr) {
323 delete fHist;
324 fHist = nullptr;
325 }
326
327 if (fNumBurnInSteps >= fChain->Size()) {
329 "MCMCInterval::CreateHist: creation of histogram failed: " <<
330 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
331 "in Markov chain." << endl;
332 return;
333 }
334
335 if (fDimension == 1) {
336 fHist = new TH1F("posterior", "MCMC Posterior Histogram",
337 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax());
338
339 } else if (fDimension == 2) {
340 fHist = new TH2F("posterior", "MCMC Posterior Histogram",
341 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax(),
342 fAxes[1]->numBins(), fAxes[1]->getMin(), fAxes[1]->getMax());
343
344 } else if (fDimension == 3) {
345 fHist = new TH3F("posterior", "MCMC Posterior Histogram",
346 fAxes[0]->numBins(), fAxes[0]->getMin(), fAxes[0]->getMax(),
347 fAxes[1]->numBins(), fAxes[1]->getMin(), fAxes[1]->getMax(),
348 fAxes[2]->numBins(), fAxes[2]->getMin(), fAxes[2]->getMax());
349
350 } else {
351 coutE(Eval) << "* Error in MCMCInterval::CreateHist() : " <<
352 "TH1* couldn't handle dimension: " << fDimension << endl;
353 return;
354 }
355
356 // Fill histogram
357 Int_t size = fChain->Size();
358 const RooArgSet* entry;
359 for (Int_t i = fNumBurnInSteps; i < size; i++) {
360 entry = fChain->Get(i);
361 if (fDimension == 1) {
362 (static_cast<TH1F*>(fHist))->Fill(entry->getRealValue(fAxes[0]->GetName()),
363 fChain->Weight());
364 } else if (fDimension == 2) {
365 (static_cast<TH2F*>(fHist))->Fill(entry->getRealValue(fAxes[0]->GetName()),
366 entry->getRealValue(fAxes[1]->GetName()),
367 fChain->Weight());
368 } else {
369 (static_cast<TH3F *>(fHist))
370 ->Fill(entry->getRealValue(fAxes[0]->GetName()), entry->getRealValue(fAxes[1]->GetName()),
371 entry->getRealValue(fAxes[2]->GetName()), fChain->Weight());
372 }
373 }
374
375 if (fDimension >= 1)
377 if (fDimension >= 2)
379 if (fDimension >= 3)
381}
382
383////////////////////////////////////////////////////////////////////////////////
384
386{
387 if (fAxes == nullptr || fChain == nullptr) {
388 coutE(InputArguments) << "* Error in MCMCInterval::CreateSparseHist(): "
389 << "Crucial data member was nullptr." << endl;
390 coutE(InputArguments) << "Make sure to fully construct/initialize."
391 << endl;
392 return;
393 }
394 if (fSparseHist != nullptr)
395 delete fSparseHist;
396
397 std::vector<double> min(fDimension);
398 std::vector<double> max(fDimension);
399 std::vector<Int_t> bins(fDimension);
400 for (Int_t i = 0; i < fDimension; i++) {
401 min[i] = fAxes[i]->getMin();
402 max[i] = fAxes[i]->getMax();
403 bins[i] = fAxes[i]->numBins();
404 }
405 fSparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
406 fDimension, bins.data(), min.data(), max.data());
407
408 // kbelasco: it appears we need to call Sumw2() just to get the
409 // histogram to keep a running total of the weight so that Getsumw doesn't
410 // just return 0
412
413 if (fNumBurnInSteps >= fChain->Size()) {
415 "MCMCInterval::CreateSparseHist: creation of histogram failed: " <<
416 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
417 "in Markov chain." << endl;
418 }
419
420 // Fill histogram
421 Int_t size = fChain->Size();
422 const RooArgSet* entry;
423 std::vector<double> x(fDimension);
424 for (Int_t i = fNumBurnInSteps; i < size; i++) {
425 entry = fChain->Get(i);
426 for (Int_t ii = 0; ii < fDimension; ii++)
427 x[ii] = entry->getRealValue(fAxes[ii]->GetName());
428 fSparseHist->Fill(x.data(), fChain->Weight());
429 }
430}
431
432////////////////////////////////////////////////////////////////////////////////
433
435{
436 if (fParameters.empty() || fChain == nullptr) {
437 coutE(Eval) << "* Error in MCMCInterval::CreateDataHist(): " <<
438 "Crucial data member was nullptr or empty." << endl;
439 coutE(Eval) << "Make sure to fully construct/initialize." << endl;
440 return;
441 }
442
443 if (fNumBurnInSteps >= fChain->Size()) {
445 "MCMCInterval::CreateDataHist: creation of histogram failed: " <<
446 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
447 "in Markov chain." << endl;
448 fDataHist = nullptr;
449 return;
450 }
451
452 fDataHist = std::unique_ptr<RooDataHist>{fChain->GetAsDataHist(SelectVars(fParameters),
453 EventRange(fNumBurnInSteps, fChain->Size()))}.release();
454}
455
456////////////////////////////////////////////////////////////////////////////////
457
459{
460 fVector.clear();
461 fVecWeight = 0;
462
463 if (fChain == nullptr) {
464 coutE(InputArguments) << "* Error in MCMCInterval::CreateVector(): " <<
465 "Crucial data member (Markov chain) was nullptr." << endl;
466 coutE(InputArguments) << "Make sure to fully construct/initialize."
467 << endl;
468 return;
469 }
470
471 if (fNumBurnInSteps >= fChain->Size()) {
473 "MCMCInterval::CreateVector: creation of vector failed: " <<
474 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
475 "in Markov chain." << endl;
476 }
477
478 // Fill vector
480 fVector.resize(size);
481 Int_t i;
482 Int_t chainIndex;
483 for (i = 0; i < size; i++) {
484 chainIndex = i + fNumBurnInSteps;
485 fVector[i] = chainIndex;
486 fVecWeight += fChain->Weight(chainIndex);
487 }
488
489 stable_sort(fVector.begin(), fVector.end(),
491}
492
493////////////////////////////////////////////////////////////////////////////////
494
496{
498 fParameters.add(parameters);
500 if (fAxes != nullptr)
501 delete[] fAxes;
503 Int_t n = 0;
504 for (auto *obj : fParameters) {
505 if (dynamic_cast<RooRealVar *>(obj) != nullptr) {
506 fAxes[n] = static_cast<RooRealVar*>(obj);
507 } else {
508 coutE(Eval) << "* Error in MCMCInterval::SetParameters: " << obj->GetName() << " not a RooRealVar*"
509 << std::endl;
510 }
511 n++;
512 }
513}
514
515////////////////////////////////////////////////////////////////////////////////
516
518{
519 switch (fIntervalType) {
520 case kShortest:
522 break;
523 case kTailFraction:
525 break;
526 default:
527 coutE(InputArguments) << "MCMCInterval::DetermineInterval(): " <<
528 "Error: Interval type not set" << endl;
529 break;
530 }
531}
532
533////////////////////////////////////////////////////////////////////////////////
534
536{
537 if (fUseKeys) {
539 } else {
541 }
542}
543
544////////////////////////////////////////////////////////////////////////////////
545
547{
548 if (fLeftSideTF < 0 || fLeftSideTF > 1) {
549 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval: "
550 << "Fraction must be in the range [0, 1]. "
551 << fLeftSideTF << "is not allowed." << endl;
552 return;
553 }
554
555 if (fDimension != 1) {
556 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval(): "
557 << "Error: Can only find a tail-fraction interval for 1-D intervals"
558 << endl;
559 return;
560 }
561
562 if (fAxes == nullptr) {
563 coutE(InputArguments) << "MCMCInterval::DetermineTailFractionInterval(): "
564 << "Crucial data member was nullptr." << endl;
565 coutE(InputArguments) << "Make sure to fully construct/initialize."
566 << endl;
567 return;
568 }
569
570 // kbelasco: fill in code here to find interval
571 //
572 // also make changes so that calling GetPosterior...() returns nullptr
573 // when fIntervalType == kTailFraction, since there really
574 // is no posterior for this type of interval determination
575 if (fVector.empty())
577
578 if (fVector.empty() || fVecWeight == 0) {
579 // if size is still 0, then creation failed.
580 // if fVecWeight == 0, then there are no entries (indicates the same
581 // error as fVector.empty() because that only happens when
582 // fNumBurnInSteps >= fChain->Size())
583 // either way, reset and return
584 fVector.clear();
585 fTFLower = -1.0 * RooNumber::infinity();
587 fTFConfLevel = 0.0;
588 fVecWeight = 0;
589 return;
590 }
591
592 RooRealVar* param = fAxes[0];
593
594 double c = fConfidenceLevel;
595 double leftTailCutoff = fVecWeight * (1 - c) * fLeftSideTF;
596 double rightTailCutoff = fVecWeight * (1 - c) * (1 - fLeftSideTF);
597 double leftTailSum = 0;
598 double rightTailSum = 0;
599
600 // kbelasco: consider changing these values to +infinity and -infinity
601 double ll = param->getMin();
602 double ul = param->getMax();
603
604 double x;
605 double w;
606
607 // save a lot of GetName() calls if compiler does not already optimize this
608 const char* name = param->GetName();
609
610 // find lower limit
611 Int_t i;
612 for (i = 0; i < (Int_t)fVector.size(); i++) {
614 w = fChain->Weight();
615
616 if (std::abs(leftTailSum + w - leftTailCutoff) <
617 std::abs(leftTailSum - leftTailCutoff)) {
618 // moving the lower limit to x would bring us closer to the desired
619 // left tail size
620 ll = x;
621 leftTailSum += w;
622 } else
623 break;
624 }
625
626 // find upper limit
627 for (i = (Int_t)fVector.size() - 1; i >= 0; i--) {
629 w = fChain->Weight();
630
631 if (std::abs(rightTailSum + w - rightTailCutoff) <
632 std::abs(rightTailSum - rightTailCutoff)) {
633 // moving the lower limit to x would bring us closer to the desired
634 // left tail size
635 ul = x;
636 rightTailSum += w;
637 } else
638 break;
639 }
640
641 fTFLower = ll;
642 fTFUpper = ul;
643 fTFConfLevel = 1 - (leftTailSum + rightTailSum) / fVecWeight;
644}
645
646////////////////////////////////////////////////////////////////////////////////
647
649{
650 if (fKeysPdf == nullptr)
652
653 if (fKeysPdf == nullptr) {
654 // if fKeysPdf is still nullptr, then it means CreateKeysPdf failed
655 // so clear all the data members this function would normally determine
656 // and return
657 fFull = 0.0;
658 fKeysCutoff = -1;
659 fKeysConfLevel = 0.0;
660 return;
661 }
662
663 // now we have a keys pdf of the posterior
664
665 double cutoff = 0.0;
666 fCutoffVar->setVal(cutoff);
667 double full = std::unique_ptr<RooAbsReal>{fProduct->createIntegral(fParameters, NormSet(fParameters))}->getVal(fParameters);
668 fFull = full;
669 if (full < 0.98) {
670 coutW(Eval) << "Warning: Integral of Keys PDF came out to " << full
671 << " instead of expected value 1. Will continue using this "
672 << "factor to normalize further integrals of this PDF." << endl;
673 }
674
675 // kbelasco: Is there a better way to set the search range?
676 // from 0 to max value of Keys
677 // kbelasco: how to get max value?
678 //double max = product.maxVal(product.getMaxVal(fParameters));
679
680 double volume = 1.0;
681 for (auto *var : static_range_cast<RooRealVar*>(fParameters))
682 volume *= (var->getMax() - var->getMin());
683
684 double topCutoff = full / volume;
685 double bottomCutoff = topCutoff;
686 double confLevel = CalcConfLevel(topCutoff, full);
687 if (AcceptableConfLevel(confLevel)) {
688 fKeysConfLevel = confLevel;
689 fKeysCutoff = topCutoff;
690 return;
691 }
692 bool changed = false;
693 // find high end of range
694 while (confLevel > fConfidenceLevel) {
695 topCutoff *= 2.0;
696 confLevel = CalcConfLevel(topCutoff, full);
697 if (AcceptableConfLevel(confLevel)) {
698 fKeysConfLevel = confLevel;
699 fKeysCutoff = topCutoff;
700 return;
701 }
702 changed = true;
703 }
704 if (changed) {
705 bottomCutoff = topCutoff / 2.0;
706 } else {
707 changed = false;
708 bottomCutoff /= 2.0;
709 confLevel = CalcConfLevel(bottomCutoff, full);
710 if (AcceptableConfLevel(confLevel)) {
711 fKeysConfLevel = confLevel;
712 fKeysCutoff = bottomCutoff;
713 return;
714 }
715 while (confLevel < fConfidenceLevel) {
716 bottomCutoff /= 2.0;
717 confLevel = CalcConfLevel(bottomCutoff, full);
718 if (AcceptableConfLevel(confLevel)) {
719 fKeysConfLevel = confLevel;
720 fKeysCutoff = bottomCutoff;
721 return;
722 }
723 changed = true;
724 }
725 if (changed) {
726 topCutoff = bottomCutoff * 2.0;
727 }
728 }
729
730 coutI(Eval) << "range set: [" << bottomCutoff << ", " << topCutoff << "]"
731 << endl;
732
733 cutoff = (topCutoff + bottomCutoff) / 2.0;
734 confLevel = CalcConfLevel(cutoff, full);
735
736 // need to use WithinDeltaFraction() because sometimes the integrating the
737 // posterior in this binary search seems to not have enough granularity to
738 // find an acceptable conf level (small no. of strange cases).
739 // WithinDeltaFraction causes the search to terminate when
740 // topCutoff is essentially equal to bottomCutoff (compared to the magnitude
741 // of their mean).
742 while (!AcceptableConfLevel(confLevel) &&
743 !WithinDeltaFraction(topCutoff, bottomCutoff)) {
744 if (confLevel > fConfidenceLevel) {
745 bottomCutoff = cutoff;
746 } else if (confLevel < fConfidenceLevel) {
747 topCutoff = cutoff;
748 }
749 cutoff = (topCutoff + bottomCutoff) / 2.0;
750 coutI(Eval) << "cutoff range: [" << bottomCutoff << ", "
751 << topCutoff << "]" << endl;
752 confLevel = CalcConfLevel(cutoff, full);
753 }
754
755 fKeysCutoff = cutoff;
756 fKeysConfLevel = confLevel;
757}
758
759////////////////////////////////////////////////////////////////////////////////
760
762{
763 if (fUseSparseHist) {
765 } else {
767 }
768}
769
770////////////////////////////////////////////////////////////////////////////////
771
773{
774 Long_t numBins;
775 if (fSparseHist == nullptr)
777
778 if (fSparseHist == nullptr) {
779 // if fSparseHist is still nullptr, then CreateSparseHist failed
780 fHistCutoff = -1;
781 fHistConfLevel = 0.0;
782 return;
783 }
784
785 numBins = (Long_t)fSparseHist->GetNbins();
786
787 std::vector<Long_t> bins(numBins);
788 for (Int_t ibin = 0; ibin < numBins; ibin++)
789 bins[ibin] = (Long_t)ibin;
790 std::stable_sort(bins.begin(), bins.end(), CompareSparseHistBins(fSparseHist));
791
792 double nEntries = fSparseHist->GetSumw();
793 double sum = 0;
794 double content;
795 Int_t i;
796 // see above note on indexing to understand numBins - 3
797 for (i = numBins - 1; i >= 0; i--) {
798 content = fSparseHist->GetBinContent(bins[i]);
799 if ((sum + content) / nEntries >= fConfidenceLevel) {
800 fHistCutoff = content;
801 if (fIsHistStrict) {
802 sum += content;
803 i--;
804 break;
805 } else {
806 i++;
807 break;
808 }
809 }
810 sum += content;
811 }
812
813 if (fIsHistStrict) {
814 // keep going to find the sum
815 for ( ; i >= 0; i--) {
816 content = fSparseHist->GetBinContent(bins[i]);
817 if (content == fHistCutoff) {
818 sum += content;
819 } else {
820 break; // content must be < fHistCutoff
821 }
822 }
823 } else {
824 // backtrack to find the cutoff and sum
825 for ( ; i < numBins; i++) {
826 content = fSparseHist->GetBinContent(bins[i]);
827 if (content > fHistCutoff) {
828 fHistCutoff = content;
829 break;
830 } else // content == fHistCutoff
831 sum -= content;
832 if (i == numBins - 1) {
833 // still haven't set fHistCutoff correctly yet, and we have no bins
834 // left, so set fHistCutoff to something higher than the tallest bin
835 fHistCutoff = content + 1.0;
836 }
837 }
838 }
839
840 fHistConfLevel = sum / nEntries;
841}
842
843////////////////////////////////////////////////////////////////////////////////
844
846{
847 Int_t numBins;
848 if (fDataHist == nullptr)
850 if (fDataHist == nullptr) {
851 // if fDataHist is still nullptr, then CreateDataHist failed
852 fHistCutoff = -1;
853 fHistConfLevel = 0.0;
854 return;
855 }
856
857 numBins = fDataHist->numEntries();
858
859 std::vector<Int_t> bins(numBins);
860 for (Int_t ibin = 0; ibin < numBins; ibin++)
861 bins[ibin] = ibin;
862 std::stable_sort(bins.begin(), bins.end(), CompareDataHistBins(fDataHist));
863
864 double nEntries = fDataHist->sum(false);
865 double sum = 0;
866 double content;
867 Int_t i;
868 for (i = numBins - 1; i >= 0; i--) {
869 fDataHist->get(bins[i]);
870 content = fDataHist->weight();
871 if ((sum + content) / nEntries >= fConfidenceLevel) {
872 fHistCutoff = content;
873 if (fIsHistStrict) {
874 sum += content;
875 i--;
876 break;
877 } else {
878 i++;
879 break;
880 }
881 }
882 sum += content;
883 }
884
885 if (fIsHistStrict) {
886 // keep going to find the sum
887 for ( ; i >= 0; i--) {
888 fDataHist->get(bins[i]);
889 content = fDataHist->weight();
890 if (content == fHistCutoff) {
891 sum += content;
892 } else {
893 break; // content must be < fHistCutoff
894 }
895 }
896 } else {
897 // backtrack to find the cutoff and sum
898 for ( ; i < numBins; i++) {
899 fDataHist->get(bins[i]);
900 content = fDataHist->weight();
901 if (content > fHistCutoff) {
902 fHistCutoff = content;
903 break;
904 } else // content == fHistCutoff
905 sum -= content;
906 if (i == numBins - 1) {
907 // still haven't set fHistCutoff correctly yet, and we have no bins
908 // left, so set fHistCutoff to something higher than the tallest bin
909 fHistCutoff = content + 1.0;
910 }
911 }
912 }
913
914 fHistConfLevel = sum / nEntries;
915}
916
917////////////////////////////////////////////////////////////////////////////////
918
920{
921 if (fIntervalType == kShortest) {
922 if (fUseKeys) {
923 return fKeysConfLevel;
924 } else {
925 return fHistConfLevel;
926 }
927 } else if (fIntervalType == kTailFraction) {
928 return fTFConfLevel;
929 } else {
930 coutE(InputArguments) << "MCMCInterval::GetActualConfidenceLevel: "
931 << "not implemented for this type of interval. Returning 0." << endl;
932 return 0;
933 }
934}
935
936////////////////////////////////////////////////////////////////////////////////
937
939{
940 switch (fIntervalType) {
941 case kShortest:
942 return LowerLimitShortest(param);
943 case kTailFraction:
944 return LowerLimitTailFraction(param);
945 default:
946 coutE(InputArguments) << "MCMCInterval::LowerLimit(): " <<
947 "Error: Interval type not set" << endl;
948 return RooNumber::infinity();
949 }
950}
951
952////////////////////////////////////////////////////////////////////////////////
953
955{
956 switch (fIntervalType) {
957 case kShortest:
958 return UpperLimitShortest(param);
959 case kTailFraction:
960 return UpperLimitTailFraction(param);
961 default:
962 coutE(InputArguments) << "MCMCInterval::UpperLimit(): " <<
963 "Error: Interval type not set" << endl;
964 return RooNumber::infinity();
965 }
966}
967
968////////////////////////////////////////////////////////////////////////////////
969
971{
972 if (fTFLower == -1.0 * RooNumber::infinity())
974
975 return fTFLower;
976}
977
978////////////////////////////////////////////////////////////////////////////////
979
981{
984
985 return fTFUpper;
986}
987
988////////////////////////////////////////////////////////////////////////////////
989
991{
992 if (fUseKeys) {
993 return LowerLimitByKeys(param);
994 } else {
995 return LowerLimitByHist(param);
996 }
997}
998
999////////////////////////////////////////////////////////////////////////////////
1000
1002{
1003 if (fUseKeys) {
1004 return UpperLimitByKeys(param);
1005 } else {
1006 return UpperLimitByHist(param);
1007 }
1008}
1009
1010////////////////////////////////////////////////////////////////////////////////
1011/// Determine the lower limit for param on this interval
1012/// using the binned data set
1013
1015{
1016 if (fUseSparseHist) {
1017 return LowerLimitBySparseHist(param);
1018 } else {
1019 return LowerLimitByDataHist(param);
1020 }
1021}
1022
1023////////////////////////////////////////////////////////////////////////////////
1024/// Determine the upper limit for param on this interval
1025/// using the binned data set
1026
1028{
1029 if (fUseSparseHist) {
1030 return UpperLimitBySparseHist(param);
1031 } else {
1032 return UpperLimitByDataHist(param);
1033 }
1034}
1035
1036////////////////////////////////////////////////////////////////////////////////
1037/// Determine the lower limit for param on this interval
1038/// using the binned data set
1039
1041{
1042 if (fDimension != 1) {
1043 coutE(InputArguments) << "In MCMCInterval::LowerLimitBySparseHist: "
1044 << "Sorry, will not compute lower limit unless dimension == 1" << endl;
1045 return param.getMin();
1046 }
1047 if (fHistCutoff < 0)
1048 DetermineBySparseHist(); // this initializes fSparseHist
1049
1050 if (fHistCutoff < 0) {
1051 // if fHistCutoff < 0 still, then determination of interval failed
1052 coutE(Eval) << "In MCMCInterval::LowerLimitBySparseHist: "
1053 << "couldn't determine cutoff. Check that num burn in steps < num "
1054 << "steps in the Markov chain. Returning param.getMin()." << endl;
1055 return param.getMin();
1056 }
1057
1058 std::vector<Int_t> coord(fDimension);
1059 for (Int_t d = 0; d < fDimension; d++) {
1060 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1061 Long_t numBins = (Long_t)fSparseHist->GetNbins();
1062 double lowerLimit = param.getMax();
1063 double val;
1064 for (Long_t i = 0; i < numBins; i++) {
1065 if (fSparseHist->GetBinContent(i, &coord[0]) >= fHistCutoff) {
1066 val = fSparseHist->GetAxis(d)->GetBinCenter(coord[d]);
1067 if (val < lowerLimit)
1068 lowerLimit = val;
1069 }
1070 }
1071 return lowerLimit;
1072 }
1073 }
1074 return param.getMin();
1075}
1076
1077////////////////////////////////////////////////////////////////////////////////
1078/// Determine the lower limit for param on this interval
1079/// using the binned data set
1080
1082{
1083 if (fHistCutoff < 0)
1084 DetermineByDataHist(); // this initializes fDataHist
1085
1086 if (fHistCutoff < 0) {
1087 // if fHistCutoff < 0 still, then determination of interval failed
1088 coutE(Eval) << "In MCMCInterval::LowerLimitByDataHist: "
1089 << "couldn't determine cutoff. Check that num burn in steps < num "
1090 << "steps in the Markov chain. Returning param.getMin()." << endl;
1091 return param.getMin();
1092 }
1093
1094 for (Int_t d = 0; d < fDimension; d++) {
1095 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1096 Int_t numBins = fDataHist->numEntries();
1097 double lowerLimit = param.getMax();
1098 double val;
1099 for (Int_t i = 0; i < numBins; i++) {
1100 fDataHist->get(i);
1101 if (fDataHist->weight() >= fHistCutoff) {
1102 val = fDataHist->get()->getRealValue(param.GetName());
1103 if (val < lowerLimit)
1104 lowerLimit = val;
1105 }
1106 }
1107 return lowerLimit;
1108 }
1109 }
1110 return param.getMin();
1111}
1112
1113////////////////////////////////////////////////////////////////////////////////
1114/// Determine the upper limit for param on this interval
1115/// using the binned data set
1116
1118{
1119 if (fDimension != 1) {
1120 coutE(InputArguments) << "In MCMCInterval::UpperLimitBySparseHist: "
1121 << "Sorry, will not compute upper limit unless dimension == 1" << endl;
1122 return param.getMax();
1123 }
1124 if (fHistCutoff < 0)
1125 DetermineBySparseHist(); // this initializes fSparseHist
1126
1127 if (fHistCutoff < 0) {
1128 // if fHistCutoff < 0 still, then determination of interval failed
1129 coutE(Eval) << "In MCMCInterval::UpperLimitBySparseHist: "
1130 << "couldn't determine cutoff. Check that num burn in steps < num "
1131 << "steps in the Markov chain. Returning param.getMax()." << endl;
1132 return param.getMax();
1133 }
1134
1135 std::vector<Int_t> coord(fDimension);
1136 for (Int_t d = 0; d < fDimension; d++) {
1137 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1138 Long_t numBins = (Long_t)fSparseHist->GetNbins();
1139 double upperLimit = param.getMin();
1140 double val;
1141 for (Long_t i = 0; i < numBins; i++) {
1142 if (fSparseHist->GetBinContent(i, &coord[0]) >= fHistCutoff) {
1143 val = fSparseHist->GetAxis(d)->GetBinCenter(coord[d]);
1144 if (val > upperLimit)
1145 upperLimit = val;
1146 }
1147 }
1148 return upperLimit;
1149 }
1150 }
1151 return param.getMax();
1152}
1153
1154////////////////////////////////////////////////////////////////////////////////
1155/// Determine the upper limit for param on this interval
1156/// using the binned data set
1157
1159{
1160 if (fHistCutoff < 0)
1161 DetermineByDataHist(); // this initializes fDataHist
1162
1163 if (fHistCutoff < 0) {
1164 // if fHistCutoff < 0 still, then determination of interval failed
1165 coutE(Eval) << "In MCMCInterval::UpperLimitByDataHist: "
1166 << "couldn't determine cutoff. Check that num burn in steps < num "
1167 << "steps in the Markov chain. Returning param.getMax()." << endl;
1168 return param.getMax();
1169 }
1170
1171 for (Int_t d = 0; d < fDimension; d++) {
1172 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1173 Int_t numBins = fDataHist->numEntries();
1174 double upperLimit = param.getMin();
1175 double val;
1176 for (Int_t i = 0; i < numBins; i++) {
1177 fDataHist->get(i);
1178 if (fDataHist->weight() >= fHistCutoff) {
1179 val = fDataHist->get()->getRealValue(param.GetName());
1180 if (val > upperLimit)
1181 upperLimit = val;
1182 }
1183 }
1184 return upperLimit;
1185 }
1186 }
1187 return param.getMax();
1188}
1189
1190////////////////////////////////////////////////////////////////////////////////
1191/// Determine the lower limit for param on this interval
1192/// using the keys pdf
1193
1195{
1196 if (fKeysCutoff < 0)
1198
1199 if (fKeysDataHist == nullptr)
1201
1202 if (fKeysCutoff < 0 || fKeysDataHist == nullptr) {
1203 // failure in determination of cutoff and/or creation of histogram
1204 coutE(Eval) << "in MCMCInterval::LowerLimitByKeys(): "
1205 << "couldn't find lower limit, check that the number of burn in "
1206 << "steps < number of total steps in the Markov chain. Returning "
1207 << "param.getMin()" << endl;
1208 return param.getMin();
1209 }
1210
1211 for (Int_t d = 0; d < fDimension; d++) {
1212 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1213 Int_t numBins = fKeysDataHist->numEntries();
1214 double lowerLimit = param.getMax();
1215 double val;
1216 for (Int_t i = 0; i < numBins; i++) {
1217 fKeysDataHist->get(i);
1218 if (fKeysDataHist->weight() >= fKeysCutoff) {
1219 val = fKeysDataHist->get()->getRealValue(param.GetName());
1220 if (val < lowerLimit)
1221 lowerLimit = val;
1222 }
1223 }
1224 return lowerLimit;
1225 }
1226 }
1227 return param.getMin();
1228}
1229
1230////////////////////////////////////////////////////////////////////////////////
1231/// Determine the upper limit for param on this interval
1232/// using the keys pdf
1233
1235{
1236 if (fKeysCutoff < 0)
1238
1239 if (fKeysDataHist == nullptr)
1241
1242 if (fKeysCutoff < 0 || fKeysDataHist == nullptr) {
1243 // failure in determination of cutoff and/or creation of histogram
1244 coutE(Eval) << "in MCMCInterval::UpperLimitByKeys(): "
1245 << "couldn't find upper limit, check that the number of burn in "
1246 << "steps < number of total steps in the Markov chain. Returning "
1247 << "param.getMax()" << endl;
1248 return param.getMax();
1249 }
1250
1251 for (Int_t d = 0; d < fDimension; d++) {
1252 if (strcmp(fAxes[d]->GetName(), param.GetName()) == 0) {
1253 Int_t numBins = fKeysDataHist->numEntries();
1254 double upperLimit = param.getMin();
1255 double val;
1256 for (Int_t i = 0; i < numBins; i++) {
1257 fKeysDataHist->get(i);
1258 if (fKeysDataHist->weight() >= fKeysCutoff) {
1259 val = fKeysDataHist->get()->getRealValue(param.GetName());
1260 if (val > upperLimit)
1261 upperLimit = val;
1262 }
1263 }
1264 return upperLimit;
1265 }
1266 }
1267 return param.getMax();
1268}
1269
1270////////////////////////////////////////////////////////////////////////////////
1271/// Determine the approximate maximum value of the Keys PDF
1272
1274{
1275 if (fKeysCutoff < 0)
1277
1278 if (fKeysDataHist == nullptr)
1280
1281 if (fKeysDataHist == nullptr) {
1282 // failure in determination of cutoff and/or creation of histogram
1283 coutE(Eval) << "in MCMCInterval::KeysMax(): "
1284 << "couldn't find Keys max value, check that the number of burn in "
1285 << "steps < number of total steps in the Markov chain. Returning 0"
1286 << endl;
1287 return 0;
1288 }
1289
1290 Int_t numBins = fKeysDataHist->numEntries();
1291 double max = 0;
1292 double w;
1293 for (Int_t i = 0; i < numBins; i++) {
1294 fKeysDataHist->get(i);
1295 w = fKeysDataHist->weight();
1296 if (w > max)
1297 max = w;
1298 }
1299
1300 return max;
1301}
1302
1303////////////////////////////////////////////////////////////////////////////////
1304
1306{
1307 if (fHistCutoff < 0)
1309
1310 return fHistCutoff;
1311}
1312
1313////////////////////////////////////////////////////////////////////////////////
1314
1316{
1317 if (fKeysCutoff < 0)
1319
1320 // kbelasco: if fFull hasn't been set (because Keys creation failed because
1321 // fNumBurnInSteps >= fChain->Size()) then this will return infinity, which
1322 // seems ok to me since it will indicate error
1323
1324 return fKeysCutoff / fFull;
1325}
1326
1327////////////////////////////////////////////////////////////////////////////////
1328
1329double MCMCInterval::CalcConfLevel(double cutoff, double full)
1330{
1331 fCutoffVar->setVal(cutoff);
1332 std::unique_ptr<RooAbsReal> integral{fProduct->createIntegral(fParameters, NormSet(fParameters))};
1333 double confLevel = integral->getVal(fParameters) / full;
1334 coutI(Eval) << "cutoff = " << cutoff << ", conf = " << confLevel << endl;
1335 return confLevel;
1336}
1337
1338////////////////////////////////////////////////////////////////////////////////
1339
1341{
1342 if (fConfidenceLevel == 0) {
1343 coutE(InputArguments) << "Error in MCMCInterval::GetPosteriorHist: "
1344 << "confidence level not set " << endl;
1345 }
1346 if (fHist == nullptr)
1347 CreateHist();
1348
1349 if (fHist == nullptr) {
1350 // if fHist is still nullptr, then CreateHist failed
1351 return nullptr;
1352 }
1353
1354 return static_cast<TH1*>(fHist->Clone("MCMCposterior_hist"));
1355}
1356
1357////////////////////////////////////////////////////////////////////////////////
1358
1360{
1361 if (fConfidenceLevel == 0) {
1362 coutE(InputArguments) << "Error in MCMCInterval::GetPosteriorKeysPdf: "
1363 << "confidence level not set " << endl;
1364 }
1365 if (fKeysPdf == nullptr)
1366 CreateKeysPdf();
1367
1368 if (fKeysPdf == nullptr) {
1369 // if fKeysPdf is still nullptr, then it means CreateKeysPdf failed
1370 return nullptr;
1371 }
1372
1373 return static_cast<RooNDKeysPdf*>(fKeysPdf->Clone("MCMCPosterior_keys"));
1374}
1375
1376////////////////////////////////////////////////////////////////////////////////
1377
1379{
1380 if (fConfidenceLevel == 0) {
1381 coutE(InputArguments) << "MCMCInterval::GetPosteriorKeysProduct: "
1382 << "confidence level not set " << endl;
1383 }
1384 if (fProduct == nullptr) {
1385 CreateKeysPdf();
1387 }
1388
1389 if (fProduct == nullptr) {
1390 // if fProduct is still nullptr, then it means CreateKeysPdf failed
1391 return nullptr;
1392 }
1393
1394 return static_cast<RooProduct*>(fProduct->Clone("MCMCPosterior_keysproduct"));
1395}
1396
1397////////////////////////////////////////////////////////////////////////////////
1398
1400{
1401 // returns list of parameters
1402 return new RooArgSet(fParameters);
1403}
1404
1405////////////////////////////////////////////////////////////////////////////////
1406
1408{
1409 return (std::abs(confLevel - fConfidenceLevel) < fEpsilon);
1410}
1411
1412////////////////////////////////////////////////////////////////////////////////
1413
1415{
1416 return (std::abs(a - b) < std::abs(fDelta * (a + b)/2));
1417}
1418
1419////////////////////////////////////////////////////////////////////////////////
1420
1422{
1423 if (fAxes == nullptr)
1424 return;
1425 if (fProduct == nullptr)
1427 if (fProduct == nullptr) {
1428 // if fProduct still nullptr, then creation failed
1429 return;
1430 }
1431
1432 //RooAbsBinning** savedBinning = new RooAbsBinning*[fDimension];
1433 std::vector<Int_t> savedBins(fDimension);
1434 Int_t i;
1435 double numBins;
1436 RooRealVar* var;
1437
1438 // kbelasco: Note - the accuracy is only increased here if the binning for
1439 // each RooRealVar is uniform
1440
1441 // kbelasco: look into why saving the binnings and replacing them doesn't
1442 // work (replaces with 1 bin always).
1443 // Note: this code modifies the binning for the parameters (if they are
1444 // uniform) and sets them back to what they were. If the binnings are not
1445 // uniform, this code does nothing.
1446
1447 // first scan through fAxes to make sure all binnings are uniform, or else
1448 // we can't change the number of bins because there seems to be an error
1449 // when setting the binning itself rather than just the number of bins
1450 bool tempChangeBinning = true;
1451 for (i = 0; i < fDimension; i++) {
1452 if (!fAxes[i]->getBinning(nullptr, false, false).isUniform()) {
1453 tempChangeBinning = false;
1454 break;
1455 }
1456 }
1457
1458 // kbelasco: for 1 dimension this should be fine, but for more dimensions
1459 // the total number of bins in the histogram increases exponentially with
1460 // the dimension, so don't do this above 1-D for now.
1461 if (fDimension >= 2)
1462 tempChangeBinning = false;
1463
1464 if (tempChangeBinning) {
1465 // set high number of bins for high accuracy on lower/upper limit by keys
1466 for (i = 0; i < fDimension; i++) {
1467 var = fAxes[i];
1468 //savedBinning[i] = &var->getBinning("__binning_clone", false, true);
1469 savedBins[i] = var->getBinning(nullptr, false, false).numBins();
1470 numBins = (var->getMax() - var->getMin()) / fEpsilon;
1471 var->setBins((Int_t)numBins);
1472 }
1473 }
1474
1475 fKeysDataHist = new RooDataHist("_productDataHist",
1476 "Keys PDF & Heaviside Product Data Hist", fParameters);
1478
1479 if (tempChangeBinning) {
1480 // set the binning back to normal
1481 for (i = 0; i < fDimension; i++) {
1482 //fAxes[i]->setBinning(*savedBinning[i], nullptr);
1483 //fAxes[i]->setBins(savedBinning[i]->numBins(), nullptr);
1484 fAxes[i]->setBins(savedBins[i], nullptr);
1485 }
1486 }
1487}
1488
1489////////////////////////////////////////////////////////////////////////////////
1490
1491bool MCMCInterval::CheckParameters(const RooArgSet& parameterPoint) const
1492{
1493 // check that the parameters are correct
1494
1495 if (parameterPoint.size() != fParameters.size() ) {
1496 coutE(Eval) << "MCMCInterval: size is wrong, parameters don't match" << std::endl;
1497 return false;
1498 }
1499 if ( ! parameterPoint.equals( fParameters ) ) {
1500 coutE(Eval) << "MCMCInterval: size is ok, but parameters don't match" << std::endl;
1501 return false;
1502 }
1503 return true;
1504}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define coutW(a)
#define coutE(a)
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:377
char name[80]
Definition TGX11.cxx:110
THnSparseT< TArrayF > THnSparseF
Definition THnSparse.h:221
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:89
Int_t numBins() const
Return number of bins.
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.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
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.
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.
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:55
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
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.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:82
Generic N-dimensional implementation of a kernel estimation p.d.f.
Provides numeric constants used in RooFit.
Definition RooNumber.h:22
static constexpr double infinity()
Return internal infinity representation.
Definition RooNumber.h:25
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
ConfInterval is an interface class for a generic interval in the RooStats framework.
Represents the Heaviside function.
Definition Heaviside.h:21
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
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
virtual void CreateVector(RooRealVar *param)
RooDataHist * fDataHist
the binned Markov Chain data
virtual double LowerLimitByDataHist(RooRealVar &param)
determine lower limit using histogram
double fDelta
topCutoff (a) considered == bottomCutoff (b) iff (std::abs(a - b) < std::abs(fDelta * (a + b)/2)); Th...
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
TH1 * fHist
the binned Markov Chain data
virtual double UpperLimitBySparseHist(RooRealVar &param)
determine upper limit using histogram
enum IntervalType fIntervalType
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
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void DetermineByKeys()
double fTFLower
lower limit of the tail-fraction interval
bool fUseSparseHist
whether to use sparse hist (vs. RooDataHist)
double fVecWeight
sum of weights of all entries in fVector
double fLeftSideTF
left side tail-fraction for interval
bool AcceptableConfLevel(double confLevel)
virtual void DetermineBySparseHist()
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
RooProduct * fProduct
the (keysPdf * heaviside) product
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
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
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,...
virtual void DetermineTailFractionInterval()
double fTFConfLevel
the actual conf level of tail-fraction interval
double fHistConfLevel
the actual conf level determined by hist
virtual void CreateSparseHist()
Heaviside * fHeaviside
the Heaviside function
double fFull
Value of intergral of fProduct.
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
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
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
double fHistCutoff
cutoff bin size to be in interval
MarkovChain * fChain
the markov chain
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 floor for the actual confidence level (strict),...
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
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
RooDataHist * fKeysDataHist
data hist representing product
RooArgSet fParameters
parameters of interest for this interval
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
double fEpsilon
acceptable error for Keys interval determination
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 RooFit::OwningPtr< RooDataHist > GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual RooFit::OwningPtr< RooDataSet > GetAsDataSet(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataSet 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 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:621
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TAxis * GetZaxis()
Definition TH1.h:326
TAxis * GetXaxis()
Definition TH1.h:324
TAxis * GetYaxis()
Definition TH1.h:325
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2752
2-D histogram with a float per channel (see TH1 documentation)
Definition TH2.h:307
3-D histogram with a float per channel (see TH1 documentation)
Definition TH3.h:317
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:37
Double_t GetBinContent(const Int_t *idx) const
Forwards to THnBase::GetBinContent() overload.
Definition THnSparse.h:121
Long64_t GetBin(const Int_t *idx) const override
Definition THnSparse.h:96
void Sumw2() override
Enable calculation of errors.
Long64_t GetNbins() const override
Definition THnSparse.h:93
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 JSONIO.h:26
@ InputArguments
Namespace for the RooStats classes.
Definition Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
RooDataHist * fDataHist
CompareDataHistBins(RooDataHist *hist)
CompareSparseHistBins(THnSparse *hist)
CompareVectorIndices(MarkovChain *chain, RooRealVar *param)
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345