Logo ROOT  
Reference Guide
MCMCIntervalPlot.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 * Project: RooStats *
6 * Package: RooFit/RooStats *
7 *************************************************************************
8 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
9 * All rights reserved. *
10 * *
11 * For the licensing terms see $ROOTSYS/LICENSE. *
12 * For the list of contributors see $ROOTSYS/README/CREDITS. *
13 *************************************************************************/
14
15/** \class RooStats::MCMCIntervalPlot
16 \ingroup Roostats
17
18This class provides simple and straightforward utilities to plot a MCMCInterval
19object. Basic use only requires a few lines once you have an MCMCInterval*:
20
21~~~ {.cpp}
22 MCMCIntervalPlot plot(*interval);
23 plot.Draw();
24~~~
25
26The standard Draw() function will currently draw the confidence interval
27range with bars if 1-D and a contour if 2-D. The MCMC posterior will also be
28plotted for the 1-D case.
29
30*/
31
35
36#include "TLine.h"
37#include "TGraph.h"
38#include "RooRealVar.h"
39#include "RooPlot.h"
40#include "TH2.h"
41#include "TH1F.h"
42#include "RooArgList.h"
43#include "TAxis.h"
44#include "RooGlobalFunc.h"
45
46#include <iostream>
47
48// Extra draw commands
49//static const char* POSTERIOR_HIST = "posterior_hist";
50//static const char* POSTERIOR_KEYS_PDF = "posterior_keys_pdf";
51//static const char* POSTERIOR_KEYS_PRODUCT = "posterior_keys_product";
52//static const char* HIST_INTERVAL = "hist_interval";
53//static const char* KEYS_PDF_INTERVAL = "keys_pdf_interval";
54//static const char* TAIL_FRACTION_INTERVAL = "tail_fraction_interval";
55//static const char* OPTION_SEP = ":";
56
58
59using namespace std;
60using namespace RooStats;
61
62////////////////////////////////////////////////////////////////////////////////
63
64MCMCIntervalPlot::MCMCIntervalPlot()
65{
66 fInterval = nullptr;
67 fParameters = nullptr;
68 fPosteriorHist = nullptr;
69 fPosteriorKeysPdf = nullptr;
70 fPosteriorKeysProduct = nullptr;
71 fDimension = 0;
74 fLineWidth = 1;
75 //fContourColor = kBlack;
76 fShowBurnIn = true;
77 fWalk = nullptr;
78 fBurnIn = nullptr;
79 fFirst = nullptr;
80 fParamGraph = nullptr;
81 fNLLGraph = nullptr;
82 fNLLHist = nullptr;
83 fWeightHist = nullptr;
84 fPosteriorHistHistCopy = nullptr;
85 fPosteriorHistTFCopy = nullptr;
86}
87
88////////////////////////////////////////////////////////////////////////////////
89
91{
92 SetMCMCInterval(interval);
93 fPosteriorHist = nullptr;
94 fPosteriorKeysPdf = nullptr;
95 fPosteriorKeysProduct = nullptr;
98 fLineWidth = 1;
99 //fContourColor = kBlack;
100 fShowBurnIn = true;
101 fWalk = nullptr;
102 fBurnIn = nullptr;
103 fFirst = nullptr;
104 fParamGraph = nullptr;
105 fNLLGraph = nullptr;
106 fNLLHist = nullptr;
107 fWeightHist = nullptr;
108 fPosteriorHistHistCopy = nullptr;
109 fPosteriorHistTFCopy = nullptr;
110}
111
112////////////////////////////////////////////////////////////////////////////////
113
115{
116 delete fParameters;
117 // kbelasco: why does deleting fPosteriorHist remove the graphics
118 // but deleting TGraphs doesn't?
119 //delete fPosteriorHist;
120 // can we delete fNLLHist and fWeightHist?
121 //delete fNLLHist;
122 //delete fWeightHist;
123
124 // kbelasco: should we delete fPosteriorKeysPdf and fPosteriorKeysProduct?
125 delete fPosteriorKeysPdf;
127
128 delete fWalk;
129 delete fBurnIn;
130 delete fFirst;
131 delete fParamGraph;
132 delete fNLLGraph;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136
138{
139 fInterval = &interval;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145
147{
148 DrawInterval(options);
149}
150
151////////////////////////////////////////////////////////////////////////////////
152
154{
155 if (fInterval->GetUseKeys())
156 DrawPosteriorKeysPdf(options);
157 else
158 DrawPosteriorHist(options);
159}
160
161////////////////////////////////////////////////////////////////////////////////
162
164 const char* title, bool scale)
165{
166 if (fPosteriorHist == nullptr)
168
169 if (fPosteriorHist == nullptr) {
170 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorHist: "
171 << "Couldn't get posterior histogram." << endl;
172 return nullptr;
173 }
174
175 // kbelasco: annoying hack because histogram drawing fails when it sees
176 // an un-recognized option like POSTERIOR_HIST, etc.
177 //const Option_t* myOpt = nullptr;
178
179 //TString tmpOpt(options);
180 //if (tmpOpt.Contains("same"))
181 // myOpt = "same";
182
183 // scale so highest bin has height 1
184 if (scale)
187
188 TString ourTitle(GetTitle());
189 if (ourTitle.CompareTo("") == 0) {
190 if (title)
191 fPosteriorHist->SetTitle(title);
192 } else
194
195 //fPosteriorHist->Draw(myOpt);
196
197 return (void*)fPosteriorHist;
198}
199
200////////////////////////////////////////////////////////////////////////////////
201
203{
204 if (fPosteriorKeysPdf == nullptr)
206
207 if (fPosteriorKeysPdf == nullptr) {
208 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
209 << "Couldn't get posterior Keys PDF." << endl;
210 return nullptr;
211 }
212
213 TString title(GetTitle());
214 bool isEmpty = (title.CompareTo("") == 0);
215
216 if (fDimension == 1) {
218 RooPlot* frame = v->frame();
219 if (frame == nullptr) {
220 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
221 << "Invalid parameter" << endl;
222 return nullptr;
223 }
224 if (isEmpty)
225 frame->SetTitle(Form("Posterior Keys PDF for %s", v->GetName()));
226 else
227 frame->SetTitle(GetTitle());
228 //fPosteriorKeysPdf->plotOn(frame);
229 //fPosteriorKeysPdf->plotOn(frame,
230 // RooFit::Normalization(1, RooAbsReal::Raw));
231 //frame->Draw(options);
232 return (void*)frame;
233 } else if (fDimension == 2) {
234 RooArgList* axes = fInterval->GetAxes();
235 RooRealVar* xVar = (RooRealVar*)axes->at(0);
236 RooRealVar* yVar = (RooRealVar*)axes->at(1);
238 "keysPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false));
239 if (isEmpty)
240 keysHist->SetTitle(
241 Form("MCMC histogram of posterior Keys PDF for %s, %s",
242 axes->at(0)->GetName(), axes->at(1)->GetName()));
243 else
244 keysHist->SetTitle(GetTitle());
245
246 keysHist->Draw(options);
247 delete axes;
248 return nullptr;
249 }
250 return nullptr;
251}
252
253////////////////////////////////////////////////////////////////////////////////
254
256{
257 switch (fInterval->GetIntervalType()) {
259 DrawShortestInterval(options);
260 break;
263 break;
264 default:
265 coutE(InputArguments) << "MCMCIntervalPlot::DrawInterval(): " <<
266 "Interval type not supported" << endl;
267 break;
268 }
269}
270
271////////////////////////////////////////////////////////////////////////////////
272
274{
275 if (fInterval->GetUseKeys())
276 DrawKeysPdfInterval(options);
277 else
278 DrawHistInterval(options);
279}
280
281////////////////////////////////////////////////////////////////////////////////
282
284{
285 TString title(GetTitle());
286 bool isEmpty = (title.CompareTo("") == 0);
287
288 if (fDimension == 1) {
289 // Draw the posterior keys PDF as well so the user can see where the
290 // limit bars line up
291 // fDimension == 1, so we know we will receive a RooPlot
292 RooPlot* frame = (RooPlot*)DrawPosteriorKeysPdf(options);
293
294 //double height = 1;
295 //double height = 2.0 * fInterval->GetKeysPdfCutoff();
296 double height = fInterval->GetKeysMax();
297
299 double ul = fInterval->UpperLimitByKeys(*p);
300 double ll = fInterval->LowerLimitByKeys(*p);
301
302 if (frame != nullptr && fPosteriorKeysPdf != nullptr) {
303 // draw shading in interval
304 if (isEmpty)
305 frame->SetTitle(nullptr);
306 else
307 frame->SetTitle(GetTitle());
308 frame->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
309 p->GetName()));
312 RooFit::Range(ll, ul, false),
317
318 // hack - this is drawn twice now:
319 // once by DrawPosteriorKeysPdf (which also configures things and sets
320 // the title), and once again here so the shading shows up behind.
323 }
324 if (frame) {
325 frame->Draw(options);
326 }
327
328 TLine* llLine = new TLine(ll, 0, ll, height);
329 TLine* ulLine = new TLine(ul, 0, ul, height);
330 llLine->SetLineColor(fLineColor);
331 ulLine->SetLineColor(fLineColor);
332 llLine->SetLineWidth(fLineWidth);
333 ulLine->SetLineWidth(fLineWidth);
334 llLine->Draw(options);
335 ulLine->Draw(options);
336 } else if (fDimension == 2) {
337 if (fPosteriorKeysPdf == nullptr)
339
340 if (fPosteriorKeysPdf == nullptr) {
341 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
342 << "Couldn't get posterior Keys PDF." << endl;
343 return;
344 }
345
346 RooArgList* axes = fInterval->GetAxes();
347 RooRealVar* xVar = (RooRealVar*)axes->at(0);
348 RooRealVar* yVar = (RooRealVar*)axes->at(1);
350 "keysContour2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false));
351 //if (isEmpty)
352 // contHist->SetTitle(Form("MCMC Keys conf. interval for %s, %s",
353 // axes->at(0)->GetName(), axes->at(1)->GetName()));
354 //else
355 // contHist->SetTitle(GetTitle());
356 if (!isEmpty)
357 contHist->SetTitle(GetTitle());
358 else
359 contHist->SetTitle(nullptr);
360
361 contHist->SetStats(false);
362
363 TString tmpOpt(options);
364 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
365
366 double cutoff = fInterval->GetKeysPdfCutoff();
367 contHist->SetContour(1, &cutoff);
368 contHist->SetLineColor(fLineColor);
369 contHist->SetLineWidth(fLineWidth);
370 contHist->Draw(tmpOpt.Data());
371 delete axes;
372 } else {
373 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
374 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
375 }
376}
377
378////////////////////////////////////////////////////////////////////////////////
379
381{
382 TString title(GetTitle());
383 bool isEmpty = (title.CompareTo("") == 0);
384
385 if (fDimension == 1) {
386 // draw lower and upper limits
388 double ul = fInterval->UpperLimitByHist(*p);
389 double ll = fInterval->LowerLimitByHist(*p);
390
391 // Draw the posterior histogram as well so the user can see where the
392 // limit bars line up
393 // fDimension == 1, so we know will get a TH1F*
394 TH1F* hist = (TH1F*)DrawPosteriorHist(options, nullptr, false);
395 if (hist == nullptr) return;
396 if (isEmpty)
397 hist->SetTitle(nullptr);
398 else
399 hist->SetTitle(GetTitle());
400 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
401 p->GetName()));
402 hist->SetStats(false);
403 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
404 double histCutoff = fInterval->GetHistCutoff();
405
406 Int_t i;
407 Int_t nBins = copy->GetNbinsX();
408 double height;
409 for (i = 1; i <= nBins; i++) {
410 // remove bins with height < cutoff
411 height = copy->GetBinContent(i);
412 if (height < histCutoff) {
413 copy->SetBinContent(i, 0);
414 copy->SetBinError(i, 0);
415 }
416 }
417
418 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
419 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
420
421 copy->SetFillStyle(1001);
423 hist->Draw(options);
424 // to show the interval area
425 copy->Draw("HIST SAME");
426
428
429 TLine* llLine = new TLine(ll, 0, ll, 1);
430 TLine* ulLine = new TLine(ul, 0, ul, 1);
431 llLine->SetLineColor(fLineColor);
432 ulLine->SetLineColor(fLineColor);
433 llLine->SetLineWidth(fLineWidth);
434 ulLine->SetLineWidth(fLineWidth);
435 llLine->Draw(options);
436 ulLine->Draw(options);
437
438 } else if (fDimension == 2) {
439 if (fPosteriorHist == nullptr)
441
442 if (fPosteriorHist == nullptr) {
443 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
444 << "Couldn't get posterior histogram." << endl;
445 return;
446 }
447
448 RooArgList* axes = fInterval->GetAxes();
449 //if (isEmpty)
450 // fPosteriorHist->SetTitle(
451 // Form("MCMC histogram conf. interval for %s, %s",
452 // axes->at(0)->GetName(), axes->at(1)->GetName()));
453 //else
454 // fPosteriorHist->SetTitle(GetTitle());
455 if (!isEmpty)
457 else
458 fPosteriorHist->SetTitle(nullptr);
459 delete axes;
460
461 fPosteriorHist->SetStats(false);
462
463 TString tmpOpt(options);
464 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
465
466 double cutoff = fInterval->GetHistCutoff();
467 fPosteriorHist->SetContour(1, &cutoff);
470 fPosteriorHist->Draw(tmpOpt.Data());
471 } else {
472 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
473 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
474 }
475}
476
477////////////////////////////////////////////////////////////////////////////////
478
480{
481 TString title(GetTitle());
482 bool isEmpty = (title.CompareTo("") == 0);
483
484 if (fDimension == 1) {
485 // Draw the posterior histogram as well so the user can see where the
486 // limit bars line up
488 double ul = fInterval->UpperLimitTailFraction(*p);
489 double ll = fInterval->LowerLimitTailFraction(*p);
490
491 TH1F* hist = (TH1F*)DrawPosteriorHist(options, nullptr, false);
492 if (hist == nullptr) return;
493 if (isEmpty)
494 hist->SetTitle(nullptr);
495 else
496 hist->SetTitle(GetTitle());
497 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
498 p->GetName()));
499 hist->SetStats(false);
500 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
501
502 Int_t i;
503 Int_t nBins = copy->GetNbinsX();
504 double center;
505 for (i = 1; i <= nBins; i++) {
506 // remove bins outside interval
507 center = copy->GetBinCenter(i);
508 if (center < ll || center > ul) {
509 copy->SetBinContent(i, 0);
510 copy->SetBinError(i, 0);
511 }
512 }
513
514 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
515 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
516
517 copy->SetFillStyle(1001);
519 hist->Draw(options);
520 copy->Draw("hist same");
521
522 // draw lower and upper limits
523 TLine* llLine = new TLine(ll, 0, ll, 1);
524 TLine* ulLine = new TLine(ul, 0, ul, 1);
525 llLine->SetLineColor(fLineColor);
526 ulLine->SetLineColor(fLineColor);
527 llLine->SetLineWidth(fLineWidth);
528 ulLine->SetLineWidth(fLineWidth);
529 llLine->Draw(options);
530 ulLine->Draw(options);
531 } else {
532 coutE(InputArguments) << "MCMCIntervalPlot::DrawTailFractionInterval: "
533 << " Sorry: " << fDimension << "-D plots not currently supported"
534 << endl;
535 }
536}
537
538////////////////////////////////////////////////////////////////////////////////
539
541{
542 if (fPosteriorKeysProduct == nullptr)
544
545 if (fPosteriorKeysProduct == nullptr) {
546 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysProduct: "
547 << "Couldn't get posterior Keys product." << endl;
548 return nullptr;
549 }
550
551 RooArgList* axes = fInterval->GetAxes();
552
553 TString title(GetTitle());
554 bool isEmpty = (title.CompareTo("") == 0);
555
556 if (fDimension == 1) {
557 RooPlot* frame = ((RooRealVar*)fParameters->first())->frame();
558 if (!frame) return nullptr;
559 if (isEmpty)
560 frame->SetTitle(Form("Posterior Keys PDF * Heaviside product for %s",
561 axes->at(0)->GetName()));
562 else
563 frame->SetTitle(GetTitle());
564 //fPosteriorKeysProduct->plotOn(frame);
567 frame->Draw(options);
568 return (void*)frame;
569 } else if (fDimension == 2) {
570 RooRealVar* xVar = (RooRealVar*)axes->at(0);
571 RooRealVar* yVar = (RooRealVar*)axes->at(1);
573 "prodPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false));
574 if (isEmpty)
575 productHist->SetTitle(
576 Form("MCMC Posterior Keys Product Hist. for %s, %s",
577 axes->at(0)->GetName(), axes->at(1)->GetName()));
578 else
579 productHist->SetTitle(GetTitle());
580 productHist->Draw(options);
581 return nullptr;
582 }
583 delete axes;
584 return nullptr;
585}
586
587////////////////////////////////////////////////////////////////////////////////
588
590{
591 const MarkovChain* markovChain = fInterval->GetChain();
592
593 Int_t size = markovChain->Size();
594 Int_t burnInSteps;
595 if (fShowBurnIn)
596 burnInSteps = fInterval->GetNumBurnInSteps();
597 else
598 burnInSteps = 0;
599
600 double* x = new double[size - burnInSteps];
601 double* y = new double[size - burnInSteps];
602 double* burnInX = nullptr;
603 double* burnInY = nullptr;
604 if (burnInSteps > 0) {
605 burnInX = new double[burnInSteps];
606 burnInY = new double[burnInSteps];
607 }
608 double firstX;
609 double firstY;
610
611 for (Int_t i = burnInSteps; i < size; i++) {
612 x[i - burnInSteps] = markovChain->Get(i)->getRealValue(xVar.GetName());
613 y[i - burnInSteps] = markovChain->Get(i)->getRealValue(yVar.GetName());
614 }
615
616 for (Int_t i = 0; i < burnInSteps; i++) {
617 burnInX[i] = markovChain->Get(i)->getRealValue(xVar.GetName());
618 burnInY[i] = markovChain->Get(i)->getRealValue(yVar.GetName());
619 }
620
621 firstX = markovChain->Get(0)->getRealValue(xVar.GetName());
622 firstY = markovChain->Get(0)->getRealValue(yVar.GetName());
623
624 TString title(GetTitle());
625 bool isEmpty = (title.CompareTo("") == 0);
626
627 TGraph* walk = new TGraph(size - burnInSteps, x, y);
628 if (isEmpty)
629 walk->SetTitle(Form("2-D Scatter Plot of Markov chain for %s, %s",
630 xVar.GetName(), yVar.GetName()));
631 else
632 walk->SetTitle(GetTitle());
633 // kbelasco: figure out how to set TGraph variable ranges
634 walk->GetXaxis()->Set(xVar.numBins(), xVar.getMin(), xVar.getMax());
635 walk->GetXaxis()->SetTitle(xVar.GetName());
636 walk->GetYaxis()->Set(yVar.numBins(), yVar.getMin(), yVar.getMax());
637 walk->GetYaxis()->SetTitle(yVar.GetName());
638 walk->SetLineColor(kGray);
639 walk->SetMarkerStyle(6);
640 walk->SetMarkerColor(kViolet);
641 walk->Draw("A,L,P,same");
642
643 TGraph* burnIn = nullptr;
644 if (burnInX != nullptr && burnInY != nullptr) {
645 burnIn = new TGraph(burnInSteps - 1, burnInX, burnInY);
646 burnIn->SetLineColor(kPink);
647 burnIn->SetMarkerStyle(6);
648 burnIn->SetMarkerColor(kPink);
649 burnIn->Draw("L,P,same");
650 }
651
652 TGraph* first = new TGraph(1, &firstX, &firstY);
653 first->SetLineColor(kGreen);
654 first->SetMarkerStyle(3);
655 first->SetMarkerSize(2);
656 first->SetMarkerColor(kGreen);
657 first->Draw("L,P,same");
658
659 //walkCanvas->Update();
660 delete [] x;
661 delete [] y;
662 if (burnInX != nullptr) delete [] burnInX;
663 if (burnInY != nullptr) delete [] burnInY;
664 //delete walk;
665 //delete burnIn;
666 //delete first;
667}
668
669////////////////////////////////////////////////////////////////////////////////
670
672{
673 const MarkovChain* markovChain = fInterval->GetChain();
674 Int_t size = markovChain->Size();
675 Int_t numEntries = 2 * size;
676 double* value = new double[numEntries];
677 double* time = new double[numEntries];
678 double val;
679 Int_t weight;
680 Int_t t = 0;
681 for (Int_t i = 0; i < size; i++) {
682 val = markovChain->Get(i)->getRealValue(param.GetName());
683 weight = (Int_t)markovChain->Weight();
684 value[2*i] = val;
685 value[2*i + 1] = val;
686 time[2*i] = t;
687 t += weight;
688 time[2*i + 1] = t;
689 }
690
691 TString title(GetTitle());
692 bool isEmpty = (title.CompareTo("") == 0);
693
694 TGraph* paramGraph = new TGraph(numEntries, time, value);
695 if (isEmpty)
696 paramGraph->SetTitle(Form("%s vs. time in Markov chain",param.GetName()));
697 else
698 paramGraph->SetTitle(GetTitle());
699 paramGraph->GetXaxis()->SetTitle("Time (discrete steps)");
700 paramGraph->GetYaxis()->SetTitle(param.GetName());
701 //paramGraph->SetLineColor(fLineColor);
702 paramGraph->Draw("A,L,same");
703 delete [] value;
704 delete [] time;
705 //gPad->Update();
706}
707
708////////////////////////////////////////////////////////////////////////////////
709
711{
712 const MarkovChain* markovChain = fInterval->GetChain();
713 Int_t size = markovChain->Size();
714 Int_t numEntries = 2 * size;
715 double* nllValue = new double[numEntries];
716 double* time = new double[numEntries];
717 double nll;
718 Int_t weight;
719 Int_t t = 0;
720 for (Int_t i = 0; i < size; i++) {
721 nll = markovChain->NLL(i);
722 weight = (Int_t)markovChain->Weight();
723 nllValue[2*i] = nll;
724 nllValue[2*i + 1] = nll;
725 time[2*i] = t;
726 t += weight;
727 time[2*i + 1] = t;
728 }
729
730 TString title(GetTitle());
731 bool isEmpty = (title.CompareTo("") == 0);
732
733 TGraph* nllGraph = new TGraph(numEntries, time, nllValue);
734 if (isEmpty)
735 nllGraph->SetTitle("NLL value vs. time in Markov chain");
736 else
737 nllGraph->SetTitle(GetTitle());
738 nllGraph->GetXaxis()->SetTitle("Time (discrete steps)");
739 nllGraph->GetYaxis()->SetTitle("NLL (-log(likelihood))");
740 //nllGraph->SetLineColor(fLineColor);
741 nllGraph->Draw("A,L,same");
742 //gPad->Update();
743 delete [] nllValue;
744 delete [] time;
745}
746
747////////////////////////////////////////////////////////////////////////////////
748
750{
751 if (fNLLHist == nullptr) {
752 const MarkovChain* markovChain = fInterval->GetChain();
753 // find the max NLL value
754 double maxNLL = 0;
755 Int_t size = markovChain->Size();
756 for (Int_t i = 0; i < size; i++)
757 if (markovChain->NLL(i) > maxNLL)
758 maxNLL = markovChain->NLL(i);
759 RooRealVar* nllVar = fInterval->GetNLLVar();
760 fNLLHist = new TH1F("mcmc_nll_hist", "MCMC NLL Histogram",
761 nllVar->getBins(), 0, maxNLL);
762 TString title(GetTitle());
763 bool isEmpty = (title.CompareTo("") == 0);
764 if (!isEmpty)
766 fNLLHist->GetXaxis()->SetTitle("-log(likelihood)");
767 for (Int_t i = 0; i < size; i++)
768 fNLLHist->Fill(markovChain->NLL(i), markovChain->Weight());
769 }
770 fNLLHist->Draw(options);
771}
772
773////////////////////////////////////////////////////////////////////////////////
774
776{
777 if (fWeightHist == nullptr) {
778 const MarkovChain* markovChain = fInterval->GetChain();
779 // find the max weight value
780 double maxWeight = 0;
781 Int_t size = markovChain->Size();
782 for (Int_t i = 0; i < size; i++)
783 if (markovChain->Weight(i) > maxWeight)
784 maxWeight = markovChain->Weight(i);
785 fWeightHist = new TH1F("mcmc_weight_hist", "MCMC Weight Histogram",
786 (Int_t)(maxWeight + 1), 0, maxWeight * 1.02);
787 for (Int_t i = 0; i < size; i++)
788 fWeightHist->Fill(markovChain->Weight(i));
789 }
790 fWeightHist->Draw(options);
791}
792
793/*
794/////////////////////////////////////////////////////////////////////
795 // 3-d plot of the parameter points
796 dataCanvas->cd(2);
797 // also plot the points in the markov chain
798 RooDataSet* markovChainData = ((MCMCInterval*)mcmcint)->GetChainAsDataSet();
799
800 TTree& chain = ((RooTreeDataStore*) markovChainData->store())->tree();
801 chain.SetMarkerStyle(6);
802 chain.SetMarkerColor(kRed);
803 chain.Draw("s:ratioSigEff:ratioBkgEff","","box"); // 3-d box proportional to posterior
804
805 // the points used in the profile construction
806 TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
807 parameterScan.SetMarkerStyle(24);
808 parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
809
810 chain.SetMarkerStyle(6);
811 chain.SetMarkerColor(kRed);
812 //chain.Draw("s:ratioSigEff:ratioBkgEff", "_MarkovChain_local_nll","box");
813 //chain.Draw("_MarkovChain_local_nll");
814////////////////////////////////////////////////////////////////////
815*/
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutE(a)
Definition: RooMsgService.h:37
int Int_t
Definition: RtypesCore.h:45
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:375
@ kGray
Definition: Rtypes.h:65
@ kPink
Definition: Rtypes.h:67
@ kBlack
Definition: Rtypes.h:65
@ kGreen
Definition: Rtypes.h:66
@ kViolet
Definition: Rtypes.h:67
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t height
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2456
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
RooAbsArg * first() const
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), 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 RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition: RooAbsPdf.h:127
Int_t numBins(const char *rangeName=nullptr) const override
virtual Int_t getBins(const char *name=nullptr) const
Get number of bins of currently defined range.
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
TH1 * createHistogram(const char *varNameList, Int_t xbins=0, Int_t ybins=0, Int_t zbins=0) const
Create and fill a ROOT histogram TH1, TH2 or TH3 with the values of this function for the variables w...
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg(), const RooCmdArg &arg10=RooCmdArg()) const
Plot (project) PDF on specified frame.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:43
void SetTitle(const char *name) override
Set the title of the RooPlot to 'title'.
Definition: RooPlot.cxx:1255
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1276
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:649
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void DrawNLLHist(const Option_t *options=nullptr)
void Draw(const Option_t *options=nullptr) override
void * DrawPosteriorHist(const Option_t *options=nullptr, const char *title=nullptr, bool scale=true)
void DrawChainScatter(RooRealVar &xVar, RooRealVar &yVar)
void DrawTailFractionInterval(const Option_t *options=nullptr)
void SetMCMCInterval(MCMCInterval &interval)
~MCMCIntervalPlot() override
Destructor of SamplingDistribution.
void DrawPosterior(const Option_t *options=nullptr)
void DrawHistInterval(const Option_t *options=nullptr)
void * DrawPosteriorKeysProduct(const Option_t *options=nullptr)
void DrawWeightHist(const Option_t *options=nullptr)
void DrawInterval(const Option_t *options=nullptr)
void DrawParameterVsTime(RooRealVar &param)
void * DrawPosteriorKeysPdf(const Option_t *options=nullptr)
void DrawShortestInterval(const Option_t *options=nullptr)
void DrawKeysPdfInterval(const Option_t *options=nullptr)
RooNDKeysPdf * fPosteriorKeysPdf
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:33
virtual bool GetUseKeys()
get whether we used kernel estimation to determine the interval
Definition: MCMCInterval.h:173
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
Definition: MCMCInterval.h:247
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
virtual double LowerLimitTailFraction(RooRealVar &param)
determine lower limit of the lower confidence interval
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList
Definition: MCMCInterval.h:102
virtual Int_t GetNumBurnInSteps()
get the number of steps in the chain to discard as burn-in,
Definition: MCMCInterval.h:179
virtual RooRealVar * GetNLLVar() const
Get a clone of the NLL variable from the markov chain.
Definition: MCMCInterval.h:221
virtual const MarkovChain * GetChain()
Get the markov chain on which this interval is based You do not own the returned MarkovChain*.
Definition: MCMCInterval.h:198
virtual double LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
virtual double LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
virtual double UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual double UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual double UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
virtual Int_t GetDimension() const
Get the number of parameters of interest in this interval.
Definition: MCMCInterval.h:194
virtual double GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual double GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
virtual double NLL(Int_t i) const
get the NLL value of entry at position i
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
virtual double Weight() const
get the weight of the current (last indexed) entry
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:49
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:759
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition: TGraph.cxx:808
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1553
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1563
void SetTitle(const char *title="") override
Change (i.e.
Definition: TGraph.cxx:2338
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:577
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:9008
void SetTitle(const char *title) override
Change (i.e.
Definition: TH1.cxx:6700
TAxis * GetXaxis()
Definition: TH1.h:322
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:9073
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3338
TAxis * GetYaxis()
Definition: TH1.h:323
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3060
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9089
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition: TH1.cxx:8444
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
virtual void SetContour(Int_t nlevels, const Double_t *levels=nullptr)
Set the number and values of contour levels.
Definition: TH1.cxx:8350
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6586
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition: TH1.cxx:2727
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8857
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:257
Use the TLine constructor to create a simple line.
Definition: TLine.h:22
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:440
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:274
Basic string class.
Definition: TString.h:136
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:451
const char * Data() const
Definition: TString.h:369
TString & Append(const char *cs)
Definition: TString.h:564
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
RooCmdArg Scaling(bool flag)
RooCmdArg FillColor(Color_t color)
RooCmdArg DrawOption(const char *opt)
RooCmdArg Range(const char *rangeName, bool adjustNorm=true)
RooCmdArg Normalization(double scaleFactor)
RooCmdArg MoveToBack()
RooCmdArg VLines()
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
Definition: first.py:1