Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 std::endl;
60using namespace RooStats;
61
62////////////////////////////////////////////////////////////////////////////////
63
65
66////////////////////////////////////////////////////////////////////////////////
67
69{
70 SetMCMCInterval(interval);
71}
72
73////////////////////////////////////////////////////////////////////////////////
74
76{
77 delete fParameters;
78 // kbelasco: why does deleting fPosteriorHist remove the graphics
79 // but deleting TGraphs doesn't?
80 //delete fPosteriorHist;
81 // can we delete fNLLHist and fWeightHist?
82 //delete fNLLHist;
83 //delete fWeightHist;
84
85 // kbelasco: should we delete fPosteriorKeysPdf and fPosteriorKeysProduct?
86 delete fPosteriorKeysPdf;
88
89 delete fWalk;
90 delete fBurnIn;
91 delete fFirst;
92 delete fParamGraph;
93 delete fNLLGraph;
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
99{
100 fInterval = &interval;
103}
104
105////////////////////////////////////////////////////////////////////////////////
106
108{
109 DrawInterval(options);
110}
111
112////////////////////////////////////////////////////////////////////////////////
113
115{
116 if (fInterval->GetUseKeys()) {
117 DrawPosteriorKeysPdf(options);
118 } else {
119 DrawPosteriorHist(options);
120 }
121}
122
123////////////////////////////////////////////////////////////////////////////////
124
126 const char* title, bool scale)
127{
128 if (fPosteriorHist == nullptr)
130
131 if (fPosteriorHist == nullptr) {
132 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorHist: "
133 << "Couldn't get posterior histogram." << endl;
134 return nullptr;
135 }
136
137 // kbelasco: annoying hack because histogram drawing fails when it sees
138 // an un-recognized option like POSTERIOR_HIST, etc.
139 //const Option_t* myOpt = nullptr;
140
141 //TString tmpOpt(options);
142 //if (tmpOpt.Contains("same"))
143 // myOpt = "same";
144
145 // scale so highest bin has height 1
146 if (scale) {
148 }
149
150 TString ourTitle(GetTitle());
151 if (ourTitle.CompareTo("") == 0) {
152 if (title)
153 fPosteriorHist->SetTitle(title);
154 } else
156
157 //fPosteriorHist->Draw(myOpt);
158
159 return (void*)fPosteriorHist;
160}
161
162////////////////////////////////////////////////////////////////////////////////
163
165{
166 if (fPosteriorKeysPdf == nullptr)
168
169 if (fPosteriorKeysPdf == nullptr) {
170 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
171 << "Couldn't get posterior Keys PDF." << endl;
172 return nullptr;
173 }
174
175 TString title(GetTitle());
176 bool isEmpty = (title.CompareTo("") == 0);
177
178 if (fDimension == 1) {
179 RooRealVar* v = static_cast<RooRealVar*>(fParameters->first());
180 RooPlot* frame = v->frame();
181 if (frame == nullptr) {
182 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
183 << "Invalid parameter" << endl;
184 return nullptr;
185 }
186 if (isEmpty) {
187 frame->SetTitle(("Posterior Keys PDF for " + std::string(v->GetName())).c_str());
188 } else {
189 frame->SetTitle(GetTitle());
190 }
191 //fPosteriorKeysPdf->plotOn(frame);
192 //fPosteriorKeysPdf->plotOn(frame,
193 // RooFit::Normalization(1, RooAbsReal::Raw));
194 //frame->Draw(options);
195 return (void*)frame;
196 } else if (fDimension == 2) {
197 RooArgList* axes = fInterval->GetAxes();
198 RooRealVar* xVar = static_cast<RooRealVar*>(axes->at(0));
199 RooRealVar* yVar = static_cast<RooRealVar*>(axes->at(1));
200 TH2F* keysHist = static_cast<TH2F*>(fPosteriorKeysPdf->createHistogram(
201 "keysPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false)));
202 if (isEmpty) {
203 keysHist->SetTitle(
204 Form("MCMC histogram of posterior Keys PDF for %s, %s",
205 axes->at(0)->GetName(), axes->at(1)->GetName()));
206 } else {
207 keysHist->SetTitle(GetTitle());
208 }
209
210 keysHist->Draw(options);
211 delete axes;
212 return nullptr;
213 }
214 return nullptr;
215}
216
217////////////////////////////////////////////////////////////////////////////////
218
220{
221 switch (fInterval->GetIntervalType()) {
223 DrawShortestInterval(options);
224 break;
227 break;
228 default:
229 coutE(InputArguments) << "MCMCIntervalPlot::DrawInterval(): " <<
230 "Interval type not supported" << endl;
231 break;
232 }
233}
234
235////////////////////////////////////////////////////////////////////////////////
236
238{
239 if (fInterval->GetUseKeys()) {
240 DrawKeysPdfInterval(options);
241 } else {
242 DrawHistInterval(options);
243 }
244}
245
246////////////////////////////////////////////////////////////////////////////////
247
249{
250 TString title(GetTitle());
251 bool isEmpty = (title.CompareTo("") == 0);
252
253 if (fDimension == 1) {
254 // Draw the posterior keys PDF as well so the user can see where the
255 // limit bars line up
256 // fDimension == 1, so we know we will receive a RooPlot
257 RooPlot* frame = reinterpret_cast<RooPlot*>(DrawPosteriorKeysPdf(options));
258
259 //double height = 1;
260 //double height = 2.0 * fInterval->GetKeysPdfCutoff();
261 double height = fInterval->GetKeysMax();
262
263 RooRealVar* p = static_cast<RooRealVar*>(fParameters->first());
264 double ul = fInterval->UpperLimitByKeys(*p);
265 double ll = fInterval->LowerLimitByKeys(*p);
266
267 if (frame != nullptr && fPosteriorKeysPdf != nullptr) {
268 // draw shading in interval
269 if (isEmpty) {
270 frame->SetTitle(nullptr);
271 } else {
272 frame->SetTitle(GetTitle());
273 }
274 frame->GetYaxis()->SetTitle(("Posterior for parameter " + std::string(p->GetName())).c_str());
277 RooFit::Range(ll, ul, false),
282
283 // hack - this is drawn twice now:
284 // once by DrawPosteriorKeysPdf (which also configures things and sets
285 // the title), and once again here so the shading shows up behind.
288 }
289 if (frame) {
290 frame->Draw(options);
291 }
292
293 TLine* llLine = new TLine(ll, 0, ll, height);
294 TLine* ulLine = new TLine(ul, 0, ul, height);
295 llLine->SetLineColor(fLineColor);
296 ulLine->SetLineColor(fLineColor);
297 llLine->SetLineWidth(fLineWidth);
298 ulLine->SetLineWidth(fLineWidth);
299 llLine->Draw(options);
300 ulLine->Draw(options);
301 } else if (fDimension == 2) {
302 if (fPosteriorKeysPdf == nullptr)
304
305 if (fPosteriorKeysPdf == nullptr) {
306 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
307 << "Couldn't get posterior Keys PDF." << endl;
308 return;
309 }
310
311 RooArgList* axes = fInterval->GetAxes();
312 RooRealVar* xVar = static_cast<RooRealVar*>(axes->at(0));
313 RooRealVar* yVar = static_cast<RooRealVar*>(axes->at(1));
314 TH2F* contHist = static_cast<TH2F*>(fPosteriorKeysPdf->createHistogram(
315 "keysContour2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false)));
316 //if (isEmpty)
317 // contHist->SetTitle(Form("MCMC Keys conf. interval for %s, %s",
318 // axes->at(0)->GetName(), axes->at(1)->GetName()));
319 //else
320 // contHist->SetTitle(GetTitle());
321 if (!isEmpty) {
322 contHist->SetTitle(GetTitle());
323 } else {
324 contHist->SetTitle(nullptr);
325 }
326
327 contHist->SetStats(false);
328
329 TString tmpOpt(options);
330 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
331
332 double cutoff = fInterval->GetKeysPdfCutoff();
333 contHist->SetContour(1, &cutoff);
334 contHist->SetLineColor(fLineColor);
335 contHist->SetLineWidth(fLineWidth);
336 contHist->Draw(tmpOpt.Data());
337 delete axes;
338 } else {
339 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
340 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
341 }
342}
343
344////////////////////////////////////////////////////////////////////////////////
345
347{
348 TString title(GetTitle());
349 bool isEmpty = (title.CompareTo("") == 0);
350
351 if (fDimension == 1) {
352 // draw lower and upper limits
353 RooRealVar* p = static_cast<RooRealVar*>(fParameters->first());
354 double ul = fInterval->UpperLimitByHist(*p);
355 double ll = fInterval->LowerLimitByHist(*p);
356
357 // Draw the posterior histogram as well so the user can see where the
358 // limit bars line up
359 // fDimension == 1, so we know will get a TH1F*
360 TH1F* hist = reinterpret_cast<TH1F*>(DrawPosteriorHist(options, nullptr, false));
361 if (hist == nullptr) return;
362 if (isEmpty) {
363 hist->SetTitle(nullptr);
364 } else {
365 hist->SetTitle(GetTitle());
366 }
367 hist->GetYaxis()->SetTitle(("Posterior for parameter " + std::string(p->GetName())).c_str());
368 hist->SetStats(false);
369 TH1F* copy = static_cast<TH1F*>(hist->Clone((std::string(hist->GetTitle()) + "_copy").c_str()));
370 double histCutoff = fInterval->GetHistCutoff();
371
372 Int_t i;
373 Int_t nBins = copy->GetNbinsX();
374 double height;
375 for (i = 1; i <= nBins; i++) {
376 // remove bins with height < cutoff
377 height = copy->GetBinContent(i);
378 if (height < histCutoff) {
379 copy->SetBinContent(i, 0);
380 copy->SetBinError(i, 0);
381 }
382 }
383
384 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
385 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
386
387 copy->SetFillStyle(1001);
389 hist->Draw(options);
390 // to show the interval area
391 copy->Draw("HIST SAME");
392
394
395 TLine* llLine = new TLine(ll, 0, ll, 1);
396 TLine* ulLine = new TLine(ul, 0, ul, 1);
397 llLine->SetLineColor(fLineColor);
398 ulLine->SetLineColor(fLineColor);
399 llLine->SetLineWidth(fLineWidth);
400 ulLine->SetLineWidth(fLineWidth);
401 llLine->Draw(options);
402 ulLine->Draw(options);
403
404 } else if (fDimension == 2) {
405 if (fPosteriorHist == nullptr)
407
408 if (fPosteriorHist == nullptr) {
409 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
410 << "Couldn't get posterior histogram." << endl;
411 return;
412 }
413
414 RooArgList* axes = fInterval->GetAxes();
415 //if (isEmpty)
416 // fPosteriorHist->SetTitle(
417 // Form("MCMC histogram conf. interval for %s, %s",
418 // axes->at(0)->GetName(), axes->at(1)->GetName()));
419 //else
420 // fPosteriorHist->SetTitle(GetTitle());
421 if (!isEmpty) {
423 } else {
424 fPosteriorHist->SetTitle(nullptr);
425 }
426 delete axes;
427
428 fPosteriorHist->SetStats(false);
429
430 TString tmpOpt(options);
431 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
432
433 double cutoff = fInterval->GetHistCutoff();
434 fPosteriorHist->SetContour(1, &cutoff);
437 fPosteriorHist->Draw(tmpOpt.Data());
438 } else {
439 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
440 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
441 }
442}
443
444////////////////////////////////////////////////////////////////////////////////
445
447{
448 TString title(GetTitle());
449 bool isEmpty = (title.CompareTo("") == 0);
450
451 if (fDimension == 1) {
452 // Draw the posterior histogram as well so the user can see where the
453 // limit bars line up
454 RooRealVar* p = static_cast<RooRealVar*>(fParameters->first());
455 double ul = fInterval->UpperLimitTailFraction(*p);
456 double ll = fInterval->LowerLimitTailFraction(*p);
457
458 TH1F* hist = reinterpret_cast<TH1F*>(DrawPosteriorHist(options, nullptr, false));
459 if (hist == nullptr) return;
460 if (isEmpty) {
461 hist->SetTitle(nullptr);
462 } else {
463 hist->SetTitle(GetTitle());
464 }
465 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
466 p->GetName()));
467 hist->SetStats(false);
468 TH1F* copy = static_cast<TH1F*>(hist->Clone(Form("%s_copy", hist->GetTitle())));
469
470 Int_t i;
471 Int_t nBins = copy->GetNbinsX();
472 double center;
473 for (i = 1; i <= nBins; i++) {
474 // remove bins outside interval
475 center = copy->GetBinCenter(i);
476 if (center < ll || center > ul) {
477 copy->SetBinContent(i, 0);
478 copy->SetBinError(i, 0);
479 }
480 }
481
482 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
483 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
484
485 copy->SetFillStyle(1001);
487 hist->Draw(options);
488 copy->Draw("hist same");
489
490 // draw lower and upper limits
491 TLine* llLine = new TLine(ll, 0, ll, 1);
492 TLine* ulLine = new TLine(ul, 0, ul, 1);
493 llLine->SetLineColor(fLineColor);
494 ulLine->SetLineColor(fLineColor);
495 llLine->SetLineWidth(fLineWidth);
496 ulLine->SetLineWidth(fLineWidth);
497 llLine->Draw(options);
498 ulLine->Draw(options);
499 } else {
500 coutE(InputArguments) << "MCMCIntervalPlot::DrawTailFractionInterval: "
501 << " Sorry: " << fDimension << "-D plots not currently supported"
502 << endl;
503 }
504}
505
506////////////////////////////////////////////////////////////////////////////////
507
509{
510 if (fPosteriorKeysProduct == nullptr)
512
513 if (fPosteriorKeysProduct == nullptr) {
514 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysProduct: "
515 << "Couldn't get posterior Keys product." << endl;
516 return nullptr;
517 }
518
519 RooArgList* axes = fInterval->GetAxes();
520
521 TString title(GetTitle());
522 bool isEmpty = (title.CompareTo("") == 0);
523
524 if (fDimension == 1) {
525 RooPlot* frame = (static_cast<RooRealVar*>(fParameters->first()))->frame();
526 if (!frame) return nullptr;
527 if (isEmpty) {
528 frame->SetTitle(Form("Posterior Keys PDF * Heaviside product for %s",
529 axes->at(0)->GetName()));
530 } else {
531 frame->SetTitle(GetTitle());
532 }
533 //fPosteriorKeysProduct->plotOn(frame);
536 frame->Draw(options);
537 return (void*)frame;
538 } else if (fDimension == 2) {
539 RooRealVar* xVar = static_cast<RooRealVar*>(axes->at(0));
540 RooRealVar* yVar = static_cast<RooRealVar*>(axes->at(1));
541 TH2F* productHist = static_cast<TH2F*>(fPosteriorKeysProduct->createHistogram(
542 "prodPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(false)));
543 if (isEmpty) {
544 productHist->SetTitle(
545 Form("MCMC Posterior Keys Product Hist. for %s, %s",
546 axes->at(0)->GetName(), axes->at(1)->GetName()));
547 } else {
548 productHist->SetTitle(GetTitle());
549 }
550 productHist->Draw(options);
551 return nullptr;
552 }
553 delete axes;
554 return nullptr;
555}
556
557////////////////////////////////////////////////////////////////////////////////
558
560{
561 const MarkovChain* markovChain = fInterval->GetChain();
562
563 Int_t size = markovChain->Size();
564 Int_t burnInSteps;
565 if (fShowBurnIn) {
566 burnInSteps = fInterval->GetNumBurnInSteps();
567 } else {
568 burnInSteps = 0;
569 }
570
571 double* x = new double[size - burnInSteps];
572 double* y = new double[size - burnInSteps];
573 double* burnInX = nullptr;
574 double* burnInY = nullptr;
575 if (burnInSteps > 0) {
576 burnInX = new double[burnInSteps];
577 burnInY = new double[burnInSteps];
578 }
579 double firstX;
580 double firstY;
581
582 for (Int_t i = burnInSteps; i < size; i++) {
583 x[i - burnInSteps] = markovChain->Get(i)->getRealValue(xVar.GetName());
584 y[i - burnInSteps] = markovChain->Get(i)->getRealValue(yVar.GetName());
585 }
586
587 for (Int_t i = 0; i < burnInSteps; i++) {
588 burnInX[i] = markovChain->Get(i)->getRealValue(xVar.GetName());
589 burnInY[i] = markovChain->Get(i)->getRealValue(yVar.GetName());
590 }
591
592 firstX = markovChain->Get(0)->getRealValue(xVar.GetName());
593 firstY = markovChain->Get(0)->getRealValue(yVar.GetName());
594
595 TString title(GetTitle());
596 bool isEmpty = (title.CompareTo("") == 0);
597
598 TGraph* walk = new TGraph(size - burnInSteps, x, y);
599 if (isEmpty) {
600 walk->SetTitle(Form("2-D Scatter Plot of Markov chain for %s, %s",
601 xVar.GetName(), yVar.GetName()));
602 } else {
603 walk->SetTitle(GetTitle());
604 }
605 // kbelasco: figure out how to set TGraph variable ranges
606 walk->GetXaxis()->Set(xVar.numBins(), xVar.getMin(), xVar.getMax());
607 walk->GetXaxis()->SetTitle(xVar.GetName());
608 walk->GetYaxis()->Set(yVar.numBins(), yVar.getMin(), yVar.getMax());
609 walk->GetYaxis()->SetTitle(yVar.GetName());
610 walk->SetLineColor(kGray);
611 walk->SetMarkerStyle(6);
612 walk->SetMarkerColor(kViolet);
613 walk->Draw("A,L,P,same");
614
615 TGraph* burnIn = nullptr;
616 if (burnInX != nullptr && burnInY != nullptr) {
617 burnIn = new TGraph(burnInSteps - 1, burnInX, burnInY);
618 burnIn->SetLineColor(kPink);
619 burnIn->SetMarkerStyle(6);
620 burnIn->SetMarkerColor(kPink);
621 burnIn->Draw("L,P,same");
622 }
623
624 TGraph* first = new TGraph(1, &firstX, &firstY);
625 first->SetLineColor(kGreen);
626 first->SetMarkerStyle(3);
627 first->SetMarkerSize(2);
628 first->SetMarkerColor(kGreen);
629 first->Draw("L,P,same");
630
631 //walkCanvas->Update();
632 delete [] x;
633 delete [] y;
634 if (burnInX != nullptr) delete [] burnInX;
635 if (burnInY != nullptr) delete [] burnInY;
636 //delete walk;
637 //delete burnIn;
638 //delete first;
639}
640
641////////////////////////////////////////////////////////////////////////////////
642
644{
645 const MarkovChain* markovChain = fInterval->GetChain();
646 Int_t size = markovChain->Size();
647 Int_t numEntries = 2 * size;
648 double* value = new double[numEntries];
649 double* time = new double[numEntries];
650 double val;
651 Int_t weight;
652 Int_t t = 0;
653 for (Int_t i = 0; i < size; i++) {
654 val = markovChain->Get(i)->getRealValue(param.GetName());
655 weight = (Int_t)markovChain->Weight();
656 value[2*i] = val;
657 value[2*i + 1] = val;
658 time[2*i] = t;
659 t += weight;
660 time[2*i + 1] = t;
661 }
662
663 TString title(GetTitle());
664 bool isEmpty = (title.CompareTo("") == 0);
665
666 TGraph* paramGraph = new TGraph(numEntries, time, value);
667 if (isEmpty) {
668 paramGraph->SetTitle(Form("%s vs. time in Markov chain",param.GetName()));
669 } else {
670 paramGraph->SetTitle(GetTitle());
671 }
672 paramGraph->GetXaxis()->SetTitle("Time (discrete steps)");
673 paramGraph->GetYaxis()->SetTitle(param.GetName());
674 //paramGraph->SetLineColor(fLineColor);
675 paramGraph->Draw("A,L,same");
676 delete [] value;
677 delete [] time;
678 //gPad->Update();
679}
680
681////////////////////////////////////////////////////////////////////////////////
682
684{
685 const MarkovChain* markovChain = fInterval->GetChain();
686 Int_t size = markovChain->Size();
687 Int_t numEntries = 2 * size;
688 double* nllValue = new double[numEntries];
689 double* time = new double[numEntries];
690 double nll;
691 Int_t weight;
692 Int_t t = 0;
693 for (Int_t i = 0; i < size; i++) {
694 nll = markovChain->NLL(i);
695 weight = (Int_t)markovChain->Weight();
696 nllValue[2*i] = nll;
697 nllValue[2*i + 1] = nll;
698 time[2*i] = t;
699 t += weight;
700 time[2*i + 1] = t;
701 }
702
703 TString title(GetTitle());
704 bool isEmpty = (title.CompareTo("") == 0);
705
706 TGraph* nllGraph = new TGraph(numEntries, time, nllValue);
707 if (isEmpty) {
708 nllGraph->SetTitle("NLL value vs. time in Markov chain");
709 } else {
710 nllGraph->SetTitle(GetTitle());
711 }
712 nllGraph->GetXaxis()->SetTitle("Time (discrete steps)");
713 nllGraph->GetYaxis()->SetTitle("NLL (-log(likelihood))");
714 //nllGraph->SetLineColor(fLineColor);
715 nllGraph->Draw("A,L,same");
716 //gPad->Update();
717 delete [] nllValue;
718 delete [] time;
719}
720
721////////////////////////////////////////////////////////////////////////////////
722
724{
725 if (fNLLHist == nullptr) {
726 const MarkovChain* markovChain = fInterval->GetChain();
727 // find the max NLL value
728 double maxNLL = 0;
729 Int_t size = markovChain->Size();
730 for (Int_t i = 0; i < size; i++) {
731 if (markovChain->NLL(i) > maxNLL)
732 maxNLL = markovChain->NLL(i);
733 }
734 RooRealVar* nllVar = fInterval->GetNLLVar();
735 fNLLHist = new TH1F("mcmc_nll_hist", "MCMC NLL Histogram",
736 nllVar->getBins(), 0, maxNLL);
737 TString title(GetTitle());
738 bool isEmpty = (title.CompareTo("") == 0);
739 if (!isEmpty)
741 fNLLHist->GetXaxis()->SetTitle("-log(likelihood)");
742 for (Int_t i = 0; i < size; i++)
743 fNLLHist->Fill(markovChain->NLL(i), markovChain->Weight());
744 }
745 fNLLHist->Draw(options);
746}
747
748////////////////////////////////////////////////////////////////////////////////
749
751{
752 if (fWeightHist == nullptr) {
753 const MarkovChain* markovChain = fInterval->GetChain();
754 // find the max weight value
755 double maxWeight = 0;
756 Int_t size = markovChain->Size();
757 for (Int_t i = 0; i < size; i++) {
758 if (markovChain->Weight(i) > maxWeight)
759 maxWeight = markovChain->Weight(i);
760 }
761 fWeightHist = new TH1F("mcmc_weight_hist", "MCMC Weight Histogram",
762 (Int_t)(maxWeight + 1), 0, maxWeight * 1.02);
763 for (Int_t i = 0; i < size; i++)
764 fWeightHist->Fill(markovChain->Weight(i));
765 }
766 fWeightHist->Draw(options);
767}
768
769/*
770/////////////////////////////////////////////////////////////////////
771 // 3-d plot of the parameter points
772 dataCanvas->cd(2);
773 // also plot the points in the markov chain
774 RooDataSet* markovChainData = ((MCMCInterval*)mcmcint)->GetChainAsDataSet();
775
776 TTree& chain = ((RooTreeDataStore*) markovChainData->store())->tree();
777 chain.SetMarkerStyle(6);
778 chain.SetMarkerColor(kRed);
779 chain.Draw("s:ratioSigEff:ratioBkgEff","","box"); // 3-d box proportional to posterior
780
781 // the points used in the profile construction
782 TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
783 parameterScan.SetMarkerStyle(24);
784 parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
785
786 chain.SetMarkerStyle(6);
787 chain.SetMarkerColor(kRed);
788 //chain.Draw("s:ratioSigEff:ratioBkgEff", "_MarkovChain_local_nll","box");
789 //chain.Draw("_MarkovChain_local_nll");
790////////////////////////////////////////////////////////////////////
791*/
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutE(a)
int Int_t
Definition RtypesCore.h:45
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
@ kGray
Definition Rtypes.h:65
@ kPink
Definition Rtypes.h:67
@ 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:2489
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={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:124
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(RooStringView 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={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) 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
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:237
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
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual bool GetUseKeys()
get whether we used kernel estimation to determine the interval
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
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
virtual Int_t GetNumBurnInSteps()
get the number of steps in the chain to discard as burn-in,
virtual RooRealVar * GetNLLVar() const
Get a clone of the NLL variable from the markov chain.
virtual const MarkovChain * GetChain()
Get the markov chain on which this interval is based You do not own the returned MarkovChain*.
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.
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 SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:45
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition TAxis.cxx:794
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:809
TAxis * GetXaxis() const
Get x axis of the graph.
Definition TGraph.cxx:1544
TAxis * GetYaxis() const
Get y axis of the graph.
Definition TGraph.cxx:1553
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2374
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:9105
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6682
TAxis * GetXaxis()
Definition TH1.h:324
virtual Int_t GetNbinsX() const
Definition TH1.h:297
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:9170
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3340
TAxis * GetYaxis()
Definition TH1.h:325
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
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:9186
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition TH1.cxx:8541
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:8447
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6568
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2748
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition TH1.cxx:8954
2-D histogram with a float per channel (see TH1 documentation)
Definition TH2.h:295
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:439
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
Basic string class.
Definition TString.h:139
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition TString.cxx:457
const char * Data() const
Definition TString.h:376
TString & Append(const char *cs)
Definition TString.h:572
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg={})
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
Namespace for the RooStats classes.
Definition Asimov.h:19