Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
89 ClassImp(RooStats::MCMCInterval);
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;
116  fKeysDataHist = NULL;
117  fCutoffVar = NULL;
118  fHist = NULL;
119  fNumBurnInSteps = 0;
120  fHistCutoff = -1;
121  fKeysCutoff = -1;
122  fDimension = 1;
123  fUseKeys = kFALSE;
124  fUseSparseHist = kFALSE;
125  fIsHistStrict = kTRUE;
126  fEpsilon = DEFAULT_EPSILON;
127  fDelta = DEFAULT_DELTA;
128  fIntervalType = kShortest;
129  fTFLower = -1.0 * RooNumber::infinity();
130  fTFUpper = RooNumber::infinity();
131  fVecWeight = 0;
132  fLeftSideTF = -1;
133 }
134 
135 ////////////////////////////////////////////////////////////////////////////////
136 
137 MCMCInterval::MCMCInterval(const char* name,
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;
159  fUseKeys = kFALSE;
160  fUseSparseHist = kFALSE;
161  fIsHistStrict = kTRUE;
162  fEpsilon = DEFAULT_EPSILON;
163  SetParameters(parameters);
164  fDelta = DEFAULT_DELTA;
165  fIntervalType = kShortest;
166  fTFLower = -1.0 * RooNumber::infinity();
167  fTFUpper = RooNumber::infinity();
168  fVecWeight = 0;
169  fLeftSideTF = -1;
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 
174 MCMCInterval::~MCMCInterval()
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  }
199  RooDataHist* fDataHist;
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  }
209  THnSparse* fSparseHist;
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  }
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 
232 Bool_t MCMCInterval::IsInInterval(const RooArgSet& point) const
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
249  RooStats::SetParameters(&point,
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 
291 void MCMCInterval::SetConfidenceLevel(Double_t cl)
292 {
293  fConfidenceLevel = cl;
294  DetermineInterval();
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 
321 void MCMCInterval::SetAxes(RooArgList& axes)
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 
337 void MCMCInterval::CreateKeysPdf()
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()) {
349  coutE(InputArguments) <<
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 
364  RooDataSet* chain = fChain->GetAsDataSet(SelectVars(fParameters),
365  EventRange(fNumBurnInSteps, fChain->Size()));
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",
376  RooArgSet(*fKeysPdf, *fHeaviside));
377 }
378 
379 ////////////////////////////////////////////////////////////////////////////////
380 
381 void MCMCInterval::CreateHist()
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()) {
393  coutE(InputArguments) <<
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 
451 void MCMCInterval::CreateSparseHist()
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()) {
484  coutE(InputArguments) <<
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 
505 void MCMCInterval::CreateDataHist()
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()) {
515  coutE(InputArguments) <<
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 
523  fDataHist = fChain->GetAsDataHist(SelectVars(fParameters),
524  EventRange(fNumBurnInSteps, fChain->Size()));
525 }
526 
527 ////////////////////////////////////////////////////////////////////////////////
528 
529 void MCMCInterval::CreateVector(RooRealVar* param)
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()) {
543  coutE(InputArguments) <<
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 {
568  fParameters.removeAll();
569  fParameters.add(parameters);
570  fDimension = fParameters.getSize();
571  if (fAxes != NULL)
572  delete[] fAxes;
573  fAxes = new RooRealVar*[fDimension];
574  TIterator* it = fParameters.createIterator();
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 
590 void MCMCInterval::DetermineInterval()
591 {
592  switch (fIntervalType) {
593  case kShortest:
594  DetermineShortestInterval();
595  break;
596  case kTailFraction:
597  DetermineTailFractionInterval();
598  break;
599  default:
600  coutE(InputArguments) << "MCMCInterval::DetermineInterval(): " <<
601  "Error: Interval type not set" << endl;
602  break;
603  }
604 }
605 
606 ////////////////////////////////////////////////////////////////////////////////
607 
608 void MCMCInterval::DetermineShortestInterval()
609 {
610  if (fUseKeys)
611  DetermineByKeys();
612  else
613  DetermineByHist();
614 }
615 
616 ////////////////////////////////////////////////////////////////////////////////
617 
618 void MCMCInterval::DetermineTailFractionInterval()
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();
658  fTFUpper = RooNumber::infinity();
659  fTFConfLevel = 0.0;
660  fVecWeight = 0;
661  return;
662  }
663 
664  RooRealVar* param = fAxes[0];
665 
666  Double_t c = fConfidenceLevel;
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 
720 void MCMCInterval::DetermineByKeys()
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);
739  RooAbsReal* integral = fProduct->createIntegral(fParameters, NormSet(fParameters));
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;
755  TIterator* it = fParameters.createIterator();
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 
837 void MCMCInterval::DetermineByHist()
838 {
839  if (fUseSparseHist)
840  DetermineBySparseHist();
841  else
842  DetermineByDataHist();
843 }
844 
845 ////////////////////////////////////////////////////////////////////////////////
846 
847 void MCMCInterval::DetermineBySparseHist()
848 {
849  Long_t numBins;
850  if (fSparseHist == NULL)
851  CreateSparseHist();
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 
918 void MCMCInterval::DetermineByDataHist()
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 
990 Double_t MCMCInterval::GetActualConfidenceLevel()
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 
1008 Double_t MCMCInterval::LowerLimit(RooRealVar& param)
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 
1024 Double_t MCMCInterval::UpperLimit(RooRealVar& param)
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 
1040 Double_t MCMCInterval::LowerLimitTailFraction(RooRealVar& /*param*/)
1041 {
1042  if (fTFLower == -1.0 * RooNumber::infinity())
1043  DetermineTailFractionInterval();
1044 
1045  return fTFLower;
1046 }
1047 
1048 ////////////////////////////////////////////////////////////////////////////////
1049 
1050 Double_t MCMCInterval::UpperLimitTailFraction(RooRealVar& /*param*/)
1051 {
1052  if (fTFUpper == RooNumber::infinity())
1053  DetermineTailFractionInterval();
1054 
1055  return fTFUpper;
1056 }
1057 
1058 ////////////////////////////////////////////////////////////////////////////////
1059 
1060 Double_t MCMCInterval::LowerLimitShortest(RooRealVar& param)
1061 {
1062  if (fUseKeys)
1063  return LowerLimitByKeys(param);
1064  else
1065  return LowerLimitByHist(param);
1066 }
1067 
1068 ////////////////////////////////////////////////////////////////////////////////
1069 
1070 Double_t MCMCInterval::UpperLimitShortest(RooRealVar& param)
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 
1082 Double_t MCMCInterval::LowerLimitByHist(RooRealVar& param)
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 
1094 Double_t MCMCInterval::UpperLimitByHist(RooRealVar& param)
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 
1106 Double_t MCMCInterval::LowerLimitBySparseHist(RooRealVar& param)
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 
1147 Double_t MCMCInterval::LowerLimitByDataHist(RooRealVar& param)
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 
1183 Double_t MCMCInterval::UpperLimitBySparseHist(RooRealVar& param)
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 
1224 Double_t MCMCInterval::UpperLimitByDataHist(RooRealVar& param)
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 
1260 Double_t MCMCInterval::LowerLimitByKeys(RooRealVar& param)
1261 {
1262  if (fKeysCutoff < 0)
1263  DetermineByKeys();
1264 
1265  if (fKeysDataHist == NULL)
1266  CreateKeysDataHist();
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 
1300 Double_t MCMCInterval::UpperLimitByKeys(RooRealVar& param)
1301 {
1302  if (fKeysCutoff < 0)
1303  DetermineByKeys();
1304 
1305  if (fKeysDataHist == NULL)
1306  CreateKeysDataHist();
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 
1339 Double_t MCMCInterval::GetKeysMax()
1340 {
1341  if (fKeysCutoff < 0)
1342  DetermineByKeys();
1343 
1344  if (fKeysDataHist == NULL)
1345  CreateKeysDataHist();
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 
1371 Double_t MCMCInterval::GetHistCutoff()
1372 {
1373  if (fHistCutoff < 0)
1374  DetermineByHist();
1375 
1376  return fHistCutoff;
1377 }
1378 
1379 ////////////////////////////////////////////////////////////////////////////////
1380 
1381 Double_t MCMCInterval::GetKeysPdfCutoff()
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 
1395 Double_t MCMCInterval::CalcConfLevel(Double_t cutoff, Double_t full)
1396 {
1397  RooAbsReal* integral;
1398  Double_t confLevel;
1399  fCutoffVar->setVal(cutoff);
1400  integral = fProduct->createIntegral(fParameters, NormSet(fParameters));
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 
1410 TH1* MCMCInterval::GetPosteriorHist()
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 
1427 RooNDKeysPdf* MCMCInterval::GetPosteriorKeysPdf()
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 
1444 RooProduct* MCMCInterval::GetPosteriorKeysProduct()
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 
1463 RooArgSet* MCMCInterval::GetParameters() const
1464 {
1465  // returns list of parameters
1466  return new RooArgSet(fParameters);
1467 }
1468 
1469 ////////////////////////////////////////////////////////////////////////////////
1470 
1471 Bool_t MCMCInterval::AcceptableConfLevel(Double_t confLevel)
1472 {
1473  return (TMath::Abs(confLevel - fConfidenceLevel) < fEpsilon);
1474 }
1475 
1476 ////////////////////////////////////////////////////////////////////////////////
1477 
1478 Bool_t MCMCInterval::WithinDeltaFraction(Double_t a, Double_t b)
1479 {
1480  return (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2));
1481 }
1482 
1483 ////////////////////////////////////////////////////////////////////////////////
1484 
1485 void MCMCInterval::CreateKeysDataHist()
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);
1540  fKeysDataHist = fProduct->fillDataHist(fKeysDataHist, &fParameters, 1.);
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 }