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