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