Logo ROOT   6.10/09
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 
62 #include "RooStats/MCMCInterval.h"
63 #include "RooStats/MarkovChain.h"
64 #include "RooStats/Heaviside.h"
65 #include "RooDataHist.h"
66 #include "RooNDKeysPdf.h"
67 #include "RooProduct.h"
68 #include "RooStats/RooStatsUtils.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 
91 using namespace RooFit;
92 using namespace RooStats;
93 using namespace std;
94 
95 static const Double_t DEFAULT_EPSILON = 0.01;
96 static const Double_t DEFAULT_DELTA = 10e-6;
97 
98 ////////////////////////////////////////////////////////////////////////////////
99 
100 MCMCInterval::MCMCInterval(const char* name)
101  : ConfInterval(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;
117  fCutoffVar = NULL;
118  fHist = NULL;
119  fNumBurnInSteps = 0;
120  fHistCutoff = -1;
121  fKeysCutoff = -1;
122  fDimension = 1;
123  fUseKeys = kFALSE;
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;
155  fCutoffVar = NULL;
156  fHist = NULL;
157  fHistCutoff = -1;
158  fKeysCutoff = -1;
159  fUseKeys = kFALSE;
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 
190 struct 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  }
200 };
201 
202 struct 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  }
210 };
211 
212 struct 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  }
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));
241  return fKeysPdf->getVal(&fParameters) >= fKeysCutoff;
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
253  Double_t* x = new Double_t[fDimension];
254  for (Int_t i = 0; i < fDimension; i++)
255  x[i] = fAxes[i]->getVal();
256  bin = fSparseHist->GetBin(x, kFALSE);
257  Double_t weight = fSparseHist->GetBinContent((Long64_t)bin);
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)
442  fHist->GetXaxis()->SetTitle(fAxes[0]->GetName());
443  if (fDimension >= 2)
444  fHist->GetYaxis()->SetTitle(fAxes[1]->GetName());
445  if (fDimension >= 3)
446  fHist->GetZaxis()->SetTitle(fAxes[2]->GetName());
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
481  fSparseHist->Sumw2();
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;
493  Double_t* x = new Double_t[fDimension];
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());
498  fSparseHist->Fill(x, fChain->Weight());
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 
566 void MCMCInterval::SetParameters(const RooArgSet& parameters)
567 {
569  fParameters.add(parameters);
571  if (fAxes != NULL)
572  delete[] fAxes;
573  fAxes = new RooRealVar*[fDimension];
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)
611  DetermineByKeys();
612  else
613  DetermineByHist();
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)
648  CreateVector(fAxes[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++) {
685  x = fChain->Get(fVector[i])->getRealValue(name);
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--) {
700  x = fChain->Get(fVector[i])->getRealValue(name);
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)
723  CreateKeysPdf();
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)
922  CreateDataHist();
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 {
1052  if (fTFUpper == RooNumber::infinity())
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)
1263  DetermineByKeys();
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)
1303  DetermineByKeys();
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)
1342  DetermineByKeys();
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)
1374  DetermineByHist();
1375 
1376  return fHistCutoff;
1377 }
1378 
1379 ////////////////////////////////////////////////////////////////////////////////
1380 
1382 {
1383  if (fKeysCutoff < 0)
1384  DetermineByKeys();
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();
1451  DetermineByKeys();
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)
1490  DetermineByKeys();
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 
1556 Bool_t MCMCInterval::CheckParameters(const RooArgSet& parameterPoint) const
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 }
virtual Double_t getMin(const char *name=0) const
virtual Bool_t IsInInterval(const RooArgSet &point) const
determine whether this point is in the confidence interval
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
MarkovChain * fChain
Definition: MCMCInterval.h:276
TIterator * createIterator(Bool_t dir=kIterForward) const
static long int sum(long int i)
Definition: Factory.cxx:2162
#define coutE(a)
Definition: RooMsgService.h:34
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
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 UpperLimitBySparseHist(RooRealVar &param)
determine upper limit using histogram
virtual Double_t LowerLimitTailFraction(RooRealVar &param)
determine lower limit of the lower confidence interval
virtual Double_t getMax(const char *name=0) const
long long Long64_t
Definition: RtypesCore.h:69
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:86
virtual const RooArgSet * get() const
Definition: RooDataHist.h:77
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins); ...
virtual RooArgSet * GetParameters() const
return a set containing the parameters of this interval the caller owns the returned RooArgSet* ...
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
virtual void CreateKeysPdf()
tomato 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:143
#define coutI(a)
Definition: RooMsgService.h:31
virtual void DetermineByKeys()
virtual void CreateVector(RooRealVar *param)
virtual Int_t numBins(const char *rangeName=0) const
virtual Double_t GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
std::vector< Int_t > fVector
Definition: MCMCInterval.h:295
RooRealVar * fCutoffVar
Definition: MCMCInterval.h:288
Long64_t GetBin(const Int_t *idx) const
Definition: THnSparse.h:94
THist< 1, float, THistStatContent, THistStatUncertainty > TH1F
Definition: THist.hxx:311
RooCmdArg NormSet(const RooArgSet &nset)
const RooAbsBinning & getBinning(const char *name=0, Bool_t verbose=kTRUE, Bool_t createOnTheFly=kFALSE) const
Return binning definition with name.
Definition: RooRealVar.cxx:267
Int_t getIndex(const RooArgSet &coord, Bool_t fast=kFALSE)
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:58
virtual Double_t UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:49
virtual void DetermineInterval()
tomato 1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:551
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically named contents.
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:33
#define NULL
Definition: RtypesCore.h:88
virtual void CreateKeysDataHist()
Long64_t GetNbins() const
Definition: THnSparse.h:91
TRObject operator()(const T1 &t1) const
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Iterator abstract base class.
Definition: TIterator.h:30
virtual void DetermineShortestInterval()
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual void CreateHist()
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.
virtual Double_t LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
RooNDKeysPdf * fKeysPdf
Definition: MCMCInterval.h:284
Double_t x[n]
Definition: legend1.C:17
Bool_t WithinDeltaFraction(Double_t a, Double_t b)
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
Efficient multidimensional histogram.
Definition: THnSparse.h:36
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual Double_t UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
Double_t GetSumw() const
Definition: THnBase.h:183
virtual Double_t weight() const
Definition: RooDataHist.h:96
void setBins(Int_t nBins, const char *name=0)
Definition: RooRealVar.h:77
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:217
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
RooDataHist * fKeysDataHist
Definition: MCMCInterval.h:287
Double_t GetBinContent(const Int_t *idx) const
Definition: THnSparse.h:116
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
virtual void DetermineByDataHist()
MCMCInterval(const char *name=0)
default constructor
Int_t getSize() const
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:501
enum IntervalType fIntervalType
Definition: MCMCInterval.h:323
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:125
RooCmdArg SelectVars(const RooArgSet &vars)
virtual Double_t GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval ...
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual Double_t LowerLimitBySparseHist(RooRealVar &param)
determine lower limit using histogram
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
tomato 2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:249
void Sumw2()
Enable calculation of errors.
Definition: THnSparse.cxx:949
virtual void DetermineByHist()
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
TAxis * GetYaxis()
Definition: TH1.h:301
RooDataHist * fDataHist
Definition: MCMCInterval.h:279
static const Double_t DEFAULT_EPSILON
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
virtual void SetConfidenceLevel(Double_t cl)
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Represents the Heaviside function.
Definition: Heaviside.h:18
Generic N-dimensional implementation of a kernel estimation p.d.f.
Definition: RooNDKeysPdf.h:45
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
const Bool_t kFALSE
Definition: RtypesCore.h:92
Double_t sum(Bool_t correctForBinSize, Bool_t inverseCorr=kFALSE) const
Return the sum of the weights of all hist bins.
virtual Double_t GetActualConfidenceLevel()
virtual Double_t GetKeysPdfCutoff() { return fKeysCutoff; }
static const Double_t DEFAULT_DELTA
long Long_t
Definition: RtypesCore.h:50
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooCmdArg EventRange(Int_t nStart, Int_t nStop)
#define ClassImp(name)
Definition: Rtypes.h:336
virtual Int_t numEntries() const
Return the number of bins.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
virtual void CreateDataHist()
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
The TH1 histogram class.
Definition: TH1.h:56
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Int_t numBins() const
Definition: RooAbsBinning.h:37
Bool_t AcceptableConfLevel(Double_t confLevel)
Double_t GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
virtual void CreateSparseHist()
THist< 3, float, THistStatContent, THistStatUncertainty > TH3F
Definition: THist.hxx:323
TAxis * GetZaxis()
Definition: TH1.h:302
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
Mother of all ROOT objects.
Definition: TObject.h:37
virtual void DetermineTailFractionInterval()
virtual Double_t UpperLimitByDataHist(RooRealVar &param)
determine upper limit using histogram
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:527
virtual void SetParameters(const RooArgSet &parameters)
Set the parameters of interest for this interval and change other internal data members accordingly...
virtual TObject * Next()=0
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2544
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Double_t UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
virtual Double_t LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void DetermineBySparseHist()
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:364
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
THist< 2, float, THistStatContent, THistStatUncertainty > TH2F
Definition: THist.hxx:317
const Bool_t kTRUE
Definition: RtypesCore.h:91
virtual Double_t UpperLimitShortest(RooRealVar &param)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
const Int_t n
Definition: legend1.C:16
virtual Double_t CalcConfLevel(Double_t cutoff, Double_t full)
static void Fill(TTree *tree, int init, int count)
Bool_t CheckParameters(const RooArgSet &point) const
check if parameters are correct. (dummy implementation to start)
virtual Double_t LowerLimitByDataHist(RooRealVar &param)
determine lower limit using histogram
TAxis * GetXaxis()
Definition: TH1.h:300