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 "TROOT.h"
37#include "TMath.h"
38#include "TLine.h"
39#include "TObjArray.h"
40#include "TList.h"
41#include "TGraph.h"
42#include "TPad.h"
43#include "RooRealVar.h"
44#include "RooPlot.h"
45#include "TH2.h"
46#include "TH1F.h"
47#include "RooArgList.h"
48#include "TAxis.h"
49#include "RooGlobalFunc.h"
50
51#include <iostream>
52
53// Extra draw commands
54//static const char* POSTERIOR_HIST = "posterior_hist";
55//static const char* POSTERIOR_KEYS_PDF = "posterior_keys_pdf";
56//static const char* POSTERIOR_KEYS_PRODUCT = "posterior_keys_product";
57//static const char* HIST_INTERVAL = "hist_interval";
58//static const char* KEYS_PDF_INTERVAL = "keys_pdf_interval";
59//static const char* TAIL_FRACTION_INTERVAL = "tail_fraction_interval";
60//static const char* OPTION_SEP = ":";
61
63
64using namespace std;
65using namespace RooStats;
66
67////////////////////////////////////////////////////////////////////////////////
68
69MCMCIntervalPlot::MCMCIntervalPlot()
70{
71 fInterval = NULL;
72 fParameters = NULL;
73 fPosteriorHist = NULL;
74 fPosteriorKeysPdf = NULL;
76 fDimension = 0;
79 fLineWidth = 1;
80 //fContourColor = kBlack;
82 fWalk = NULL;
83 fBurnIn = NULL;
84 fFirst = NULL;
85 fParamGraph = NULL;
86 fNLLGraph = NULL;
87 fNLLHist = NULL;
88 fWeightHist = NULL;
91}
92
93////////////////////////////////////////////////////////////////////////////////
94
96{
97 SetMCMCInterval(interval);
98 fPosteriorHist = NULL;
99 fPosteriorKeysPdf = NULL;
103 fLineWidth = 1;
104 //fContourColor = kBlack;
106 fWalk = NULL;
107 fBurnIn = NULL;
108 fFirst = NULL;
109 fParamGraph = NULL;
110 fNLLGraph = NULL;
111 fNLLHist = NULL;
112 fWeightHist = NULL;
115}
116
117////////////////////////////////////////////////////////////////////////////////
118
120{
121 delete fParameters;
122 // kbelasco: why does deleting fPosteriorHist remove the graphics
123 // but deleting TGraphs doesn't?
124 //delete fPosteriorHist;
125 // can we delete fNLLHist and fWeightHist?
126 //delete fNLLHist;
127 //delete fWeightHist;
128
129 // kbelasco: should we delete fPosteriorKeysPdf and fPosteriorKeysProduct?
130 delete fPosteriorKeysPdf;
132
133 delete fWalk;
134 delete fBurnIn;
135 delete fFirst;
136 delete fParamGraph;
137 delete fNLLGraph;
138}
139
140////////////////////////////////////////////////////////////////////////////////
141
143{
144 fInterval = &interval;
147}
148
149////////////////////////////////////////////////////////////////////////////////
150
152{
153 DrawInterval(options);
154}
155
156////////////////////////////////////////////////////////////////////////////////
157
159{
160 if (fInterval->GetUseKeys())
161 DrawPosteriorKeysPdf(options);
162 else
163 DrawPosteriorHist(options);
164}
165
166////////////////////////////////////////////////////////////////////////////////
167
169 const char* title, Bool_t scale)
170{
171 if (fPosteriorHist == NULL)
173
174 if (fPosteriorHist == NULL) {
175 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorHist: "
176 << "Couldn't get posterior histogram." << endl;
177 return NULL;
178 }
179
180 // kbelasco: annoying hack because histogram drawing fails when it sees
181 // an un-recognized option like POSTERIOR_HIST, etc.
182 //const Option_t* myOpt = NULL;
183
184 //TString tmpOpt(options);
185 //if (tmpOpt.Contains("same"))
186 // myOpt = "same";
187
188 // scale so highest bin has height 1
189 if (scale)
192
193 TString ourTitle(GetTitle());
194 if (ourTitle.CompareTo("") == 0) {
195 if (title)
196 fPosteriorHist->SetTitle(title);
197 } else
199
200 //fPosteriorHist->Draw(myOpt);
201
202 return (void*)fPosteriorHist;
203}
204
205////////////////////////////////////////////////////////////////////////////////
206
208{
209 if (fPosteriorKeysPdf == NULL)
211
212 if (fPosteriorKeysPdf == NULL) {
213 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
214 << "Couldn't get posterior Keys PDF." << endl;
215 return NULL;
216 }
217
218 TString title(GetTitle());
219 Bool_t isEmpty = (title.CompareTo("") == 0);
220
221 if (fDimension == 1) {
223 RooPlot* frame = v->frame();
224 if (frame == NULL) {
225 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
226 << "Invalid parameter" << endl;
227 return NULL;
228 }
229 if (isEmpty)
230 frame->SetTitle(Form("Posterior Keys PDF for %s", v->GetName()));
231 else
232 frame->SetTitle(GetTitle());
233 //fPosteriorKeysPdf->plotOn(frame);
234 //fPosteriorKeysPdf->plotOn(frame,
235 // RooFit::Normalization(1, RooAbsReal::Raw));
236 //frame->Draw(options);
237 return (void*)frame;
238 } else if (fDimension == 2) {
239 RooArgList* axes = fInterval->GetAxes();
240 RooRealVar* xVar = (RooRealVar*)axes->at(0);
241 RooRealVar* yVar = (RooRealVar*)axes->at(1);
243 "keysPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
244 if (isEmpty)
245 keysHist->SetTitle(
246 Form("MCMC histogram of posterior Keys PDF for %s, %s",
247 axes->at(0)->GetName(), axes->at(1)->GetName()));
248 else
249 keysHist->SetTitle(GetTitle());
250
251 keysHist->Draw(options);
252 delete axes;
253 return NULL;
254 }
255 return NULL;
256}
257
258////////////////////////////////////////////////////////////////////////////////
259
261{
262 switch (fInterval->GetIntervalType()) {
264 DrawShortestInterval(options);
265 break;
268 break;
269 default:
270 coutE(InputArguments) << "MCMCIntervalPlot::DrawInterval(): " <<
271 "Interval type not supported" << endl;
272 break;
273 }
274}
275
276////////////////////////////////////////////////////////////////////////////////
277
279{
280 if (fInterval->GetUseKeys())
281 DrawKeysPdfInterval(options);
282 else
283 DrawHistInterval(options);
284}
285
286////////////////////////////////////////////////////////////////////////////////
287
289{
290 TString title(GetTitle());
291 Bool_t isEmpty = (title.CompareTo("") == 0);
292
293 if (fDimension == 1) {
294 // Draw the posterior keys PDF as well so the user can see where the
295 // limit bars line up
296 // fDimension == 1, so we know we will receive a RooPlot
297 RooPlot* frame = (RooPlot*)DrawPosteriorKeysPdf(options);
298
299 //Double_t height = 1;
300 //Double_t height = 2.0 * fInterval->GetKeysPdfCutoff();
301 Double_t height = fInterval->GetKeysMax();
302
306
307 if (frame != NULL && fPosteriorKeysPdf != NULL) {
308 // draw shading in interval
309 if (isEmpty)
310 frame->SetTitle(NULL);
311 else
312 frame->SetTitle(GetTitle());
313 frame->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
314 p->GetName()));
317 RooFit::Range(ll, ul, kFALSE),
322
323 // hack - this is drawn twice now:
324 // once by DrawPosteriorKeysPdf (which also configures things and sets
325 // the title), and once again here so the shading shows up behind.
328 }
329 if (frame) {
330 frame->Draw(options);
331 }
332
333 TLine* llLine = new TLine(ll, 0, ll, height);
334 TLine* ulLine = new TLine(ul, 0, ul, height);
335 llLine->SetLineColor(fLineColor);
336 ulLine->SetLineColor(fLineColor);
337 llLine->SetLineWidth(fLineWidth);
338 ulLine->SetLineWidth(fLineWidth);
339 llLine->Draw(options);
340 ulLine->Draw(options);
341 } else if (fDimension == 2) {
342 if (fPosteriorKeysPdf == NULL)
344
345 if (fPosteriorKeysPdf == NULL) {
346 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
347 << "Couldn't get posterior Keys PDF." << endl;
348 return;
349 }
350
351 RooArgList* axes = fInterval->GetAxes();
352 RooRealVar* xVar = (RooRealVar*)axes->at(0);
353 RooRealVar* yVar = (RooRealVar*)axes->at(1);
355 "keysContour2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
356 //if (isEmpty)
357 // contHist->SetTitle(Form("MCMC Keys conf. interval for %s, %s",
358 // axes->at(0)->GetName(), axes->at(1)->GetName()));
359 //else
360 // contHist->SetTitle(GetTitle());
361 if (!isEmpty)
362 contHist->SetTitle(GetTitle());
363 else
364 contHist->SetTitle(NULL);
365
366 contHist->SetStats(kFALSE);
367
368 TString tmpOpt(options);
369 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
370
372 contHist->SetContour(1, &cutoff);
373 contHist->SetLineColor(fLineColor);
374 contHist->SetLineWidth(fLineWidth);
375 contHist->Draw(tmpOpt.Data());
376 delete axes;
377 } else {
378 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
379 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
380 }
381}
382
383////////////////////////////////////////////////////////////////////////////////
384
386{
387 TString title(GetTitle());
388 Bool_t isEmpty = (title.CompareTo("") == 0);
389
390 if (fDimension == 1) {
391 // draw lower and upper limits
395
396 // Draw the posterior histogram as well so the user can see where the
397 // limit bars line up
398 // fDimension == 1, so we know will get a TH1F*
399 TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
400 if (hist == NULL) return;
401 if (isEmpty)
402 hist->SetTitle(NULL);
403 else
404 hist->SetTitle(GetTitle());
405 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
406 p->GetName()));
407 hist->SetStats(kFALSE);
408 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
409 Double_t histCutoff = fInterval->GetHistCutoff();
410
411 Int_t i;
412 Int_t nBins = copy->GetNbinsX();
413 Double_t height;
414 for (i = 1; i <= nBins; i++) {
415 // remove bins with height < cutoff
416 height = copy->GetBinContent(i);
417 if (height < histCutoff) {
418 copy->SetBinContent(i, 0);
419 copy->SetBinError(i, 0);
420 }
421 }
422
423 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
424 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
425
426 copy->SetFillStyle(1001);
428 hist->Draw(options);
429 // to show the interval area
430 copy->Draw("HIST SAME");
431
433
434 TLine* llLine = new TLine(ll, 0, ll, 1);
435 TLine* ulLine = new TLine(ul, 0, ul, 1);
436 llLine->SetLineColor(fLineColor);
437 ulLine->SetLineColor(fLineColor);
438 llLine->SetLineWidth(fLineWidth);
439 ulLine->SetLineWidth(fLineWidth);
440 llLine->Draw(options);
441 ulLine->Draw(options);
442
443 } else if (fDimension == 2) {
444 if (fPosteriorHist == NULL)
446
447 if (fPosteriorHist == NULL) {
448 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
449 << "Couldn't get posterior histogram." << endl;
450 return;
451 }
452
453 RooArgList* axes = fInterval->GetAxes();
454 //if (isEmpty)
455 // fPosteriorHist->SetTitle(
456 // Form("MCMC histogram conf. interval for %s, %s",
457 // axes->at(0)->GetName(), axes->at(1)->GetName()));
458 //else
459 // fPosteriorHist->SetTitle(GetTitle());
460 if (!isEmpty)
462 else
464 delete axes;
465
467
468 TString tmpOpt(options);
469 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
470
471 Double_t cutoff = fInterval->GetHistCutoff();
472 fPosteriorHist->SetContour(1, &cutoff);
475 fPosteriorHist->Draw(tmpOpt.Data());
476 } else {
477 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
478 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
479 }
480}
481
482////////////////////////////////////////////////////////////////////////////////
483
485{
486 TString title(GetTitle());
487 Bool_t isEmpty = (title.CompareTo("") == 0);
488
489 if (fDimension == 1) {
490 // Draw the posterior histogram as well so the user can see where the
491 // limit bars line up
495
496 TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
497 if (hist == NULL) return;
498 if (isEmpty)
499 hist->SetTitle(NULL);
500 else
501 hist->SetTitle(GetTitle());
502 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
503 p->GetName()));
504 hist->SetStats(kFALSE);
505 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
506
507 Int_t i;
508 Int_t nBins = copy->GetNbinsX();
509 Double_t center;
510 for (i = 1; i <= nBins; i++) {
511 // remove bins outside interval
512 center = copy->GetBinCenter(i);
513 if (center < ll || center > ul) {
514 copy->SetBinContent(i, 0);
515 copy->SetBinError(i, 0);
516 }
517 }
518
519 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
520 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
521
522 copy->SetFillStyle(1001);
524 hist->Draw(options);
525 copy->Draw("hist same");
526
527 // draw lower and upper limits
528 TLine* llLine = new TLine(ll, 0, ll, 1);
529 TLine* ulLine = new TLine(ul, 0, ul, 1);
530 llLine->SetLineColor(fLineColor);
531 ulLine->SetLineColor(fLineColor);
532 llLine->SetLineWidth(fLineWidth);
533 ulLine->SetLineWidth(fLineWidth);
534 llLine->Draw(options);
535 ulLine->Draw(options);
536 } else {
537 coutE(InputArguments) << "MCMCIntervalPlot::DrawTailFractionInterval: "
538 << " Sorry: " << fDimension << "-D plots not currently supported"
539 << endl;
540 }
541}
542
543////////////////////////////////////////////////////////////////////////////////
544
546{
547 if (fPosteriorKeysProduct == NULL)
549
550 if (fPosteriorKeysProduct == NULL) {
551 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysProduct: "
552 << "Couldn't get posterior Keys product." << endl;
553 return NULL;
554 }
555
556 RooArgList* axes = fInterval->GetAxes();
557
558 TString title(GetTitle());
559 Bool_t isEmpty = (title.CompareTo("") == 0);
560
561 if (fDimension == 1) {
562 RooPlot* frame = ((RooRealVar*)fParameters->first())->frame();
563 if (!frame) return NULL;
564 if (isEmpty)
565 frame->SetTitle(Form("Posterior Keys PDF * Heaviside product for %s",
566 axes->at(0)->GetName()));
567 else
568 frame->SetTitle(GetTitle());
569 //fPosteriorKeysProduct->plotOn(frame);
572 frame->Draw(options);
573 return (void*)frame;
574 } else if (fDimension == 2) {
575 RooRealVar* xVar = (RooRealVar*)axes->at(0);
576 RooRealVar* yVar = (RooRealVar*)axes->at(1);
578 "prodPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
579 if (isEmpty)
580 productHist->SetTitle(
581 Form("MCMC Posterior Keys Product Hist. for %s, %s",
582 axes->at(0)->GetName(), axes->at(1)->GetName()));
583 else
584 productHist->SetTitle(GetTitle());
585 productHist->Draw(options);
586 return NULL;
587 }
588 delete axes;
589 return NULL;
590}
591
592////////////////////////////////////////////////////////////////////////////////
593
595{
596 const MarkovChain* markovChain = fInterval->GetChain();
597
598 Int_t size = markovChain->Size();
599 Int_t burnInSteps;
600 if (fShowBurnIn)
601 burnInSteps = fInterval->GetNumBurnInSteps();
602 else
603 burnInSteps = 0;
604
605 Double_t* x = new Double_t[size - burnInSteps];
606 Double_t* y = new Double_t[size - burnInSteps];
607 Double_t* burnInX = NULL;
608 Double_t* burnInY = NULL;
609 if (burnInSteps > 0) {
610 burnInX = new Double_t[burnInSteps];
611 burnInY = new Double_t[burnInSteps];
612 }
613 Double_t firstX;
614 Double_t firstY;
615
616 for (Int_t i = burnInSteps; i < size; i++) {
617 x[i - burnInSteps] = markovChain->Get(i)->getRealValue(xVar.GetName());
618 y[i - burnInSteps] = markovChain->Get(i)->getRealValue(yVar.GetName());
619 }
620
621 for (Int_t i = 0; i < burnInSteps; i++) {
622 burnInX[i] = markovChain->Get(i)->getRealValue(xVar.GetName());
623 burnInY[i] = markovChain->Get(i)->getRealValue(yVar.GetName());
624 }
625
626 firstX = markovChain->Get(0)->getRealValue(xVar.GetName());
627 firstY = markovChain->Get(0)->getRealValue(yVar.GetName());
628
629 TString title(GetTitle());
630 Bool_t isEmpty = (title.CompareTo("") == 0);
631
632 TGraph* walk = new TGraph(size - burnInSteps, x, y);
633 if (isEmpty)
634 walk->SetTitle(Form("2-D Scatter Plot of Markov chain for %s, %s",
635 xVar.GetName(), yVar.GetName()));
636 else
637 walk->SetTitle(GetTitle());
638 // kbelasco: figure out how to set TGraph variable ranges
639 walk->GetXaxis()->Set(xVar.numBins(), xVar.getMin(), xVar.getMax());
640 walk->GetXaxis()->SetTitle(xVar.GetName());
641 walk->GetYaxis()->Set(yVar.numBins(), yVar.getMin(), yVar.getMax());
642 walk->GetYaxis()->SetTitle(yVar.GetName());
643 walk->SetLineColor(kGray);
644 walk->SetMarkerStyle(6);
645 walk->SetMarkerColor(kViolet);
646 walk->Draw("A,L,P,same");
647
648 TGraph* burnIn = NULL;
649 if (burnInX != NULL && burnInY != NULL) {
650 burnIn = new TGraph(burnInSteps - 1, burnInX, burnInY);
651 burnIn->SetLineColor(kPink);
652 burnIn->SetMarkerStyle(6);
653 burnIn->SetMarkerColor(kPink);
654 burnIn->Draw("L,P,same");
655 }
656
657 TGraph* first = new TGraph(1, &firstX, &firstY);
658 first->SetLineColor(kGreen);
659 first->SetMarkerStyle(3);
660 first->SetMarkerSize(2);
661 first->SetMarkerColor(kGreen);
662 first->Draw("L,P,same");
663
664 //walkCanvas->Update();
665 delete [] x;
666 delete [] y;
667 if (burnInX != NULL) delete [] burnInX;
668 if (burnInY != NULL) delete [] burnInY;
669 //delete walk;
670 //delete burnIn;
671 //delete first;
672}
673
674////////////////////////////////////////////////////////////////////////////////
675
677{
678 const MarkovChain* markovChain = fInterval->GetChain();
679 Int_t size = markovChain->Size();
680 Int_t numEntries = 2 * size;
681 Double_t* value = new Double_t[numEntries];
682 Double_t* time = new Double_t[numEntries];
683 Double_t val;
684 Int_t weight;
685 Int_t t = 0;
686 for (Int_t i = 0; i < size; i++) {
687 val = markovChain->Get(i)->getRealValue(param.GetName());
688 weight = (Int_t)markovChain->Weight();
689 value[2*i] = val;
690 value[2*i + 1] = val;
691 time[2*i] = t;
692 t += weight;
693 time[2*i + 1] = t;
694 }
695
696 TString title(GetTitle());
697 Bool_t isEmpty = (title.CompareTo("") == 0);
698
699 TGraph* paramGraph = new TGraph(numEntries, time, value);
700 if (isEmpty)
701 paramGraph->SetTitle(Form("%s vs. time in Markov chain",param.GetName()));
702 else
703 paramGraph->SetTitle(GetTitle());
704 paramGraph->GetXaxis()->SetTitle("Time (discrete steps)");
705 paramGraph->GetYaxis()->SetTitle(param.GetName());
706 //paramGraph->SetLineColor(fLineColor);
707 paramGraph->Draw("A,L,same");
708 delete [] value;
709 delete [] time;
710 //gPad->Update();
711}
712
713////////////////////////////////////////////////////////////////////////////////
714
716{
717 const MarkovChain* markovChain = fInterval->GetChain();
718 Int_t size = markovChain->Size();
719 Int_t numEntries = 2 * size;
720 Double_t* nllValue = new Double_t[numEntries];
721 Double_t* time = new Double_t[numEntries];
722 Double_t nll;
723 Int_t weight;
724 Int_t t = 0;
725 for (Int_t i = 0; i < size; i++) {
726 nll = markovChain->NLL(i);
727 weight = (Int_t)markovChain->Weight();
728 nllValue[2*i] = nll;
729 nllValue[2*i + 1] = nll;
730 time[2*i] = t;
731 t += weight;
732 time[2*i + 1] = t;
733 }
734
735 TString title(GetTitle());
736 Bool_t isEmpty = (title.CompareTo("") == 0);
737
738 TGraph* nllGraph = new TGraph(numEntries, time, nllValue);
739 if (isEmpty)
740 nllGraph->SetTitle("NLL value vs. time in Markov chain");
741 else
742 nllGraph->SetTitle(GetTitle());
743 nllGraph->GetXaxis()->SetTitle("Time (discrete steps)");
744 nllGraph->GetYaxis()->SetTitle("NLL (-log(likelihood))");
745 //nllGraph->SetLineColor(fLineColor);
746 nllGraph->Draw("A,L,same");
747 //gPad->Update();
748 delete [] nllValue;
749 delete [] time;
750}
751
752////////////////////////////////////////////////////////////////////////////////
753
755{
756 if (fNLLHist == NULL) {
757 const MarkovChain* markovChain = fInterval->GetChain();
758 // find the max NLL value
759 Double_t maxNLL = 0;
760 Int_t size = markovChain->Size();
761 for (Int_t i = 0; i < size; i++)
762 if (markovChain->NLL(i) > maxNLL)
763 maxNLL = markovChain->NLL(i);
764 RooRealVar* nllVar = fInterval->GetNLLVar();
765 fNLLHist = new TH1F("mcmc_nll_hist", "MCMC NLL Histogram",
766 nllVar->getBins(), 0, maxNLL);
767 TString title(GetTitle());
768 Bool_t isEmpty = (title.CompareTo("") == 0);
769 if (!isEmpty)
771 fNLLHist->GetXaxis()->SetTitle("-log(likelihood)");
772 for (Int_t i = 0; i < size; i++)
773 fNLLHist->Fill(markovChain->NLL(i), markovChain->Weight());
774 }
775 fNLLHist->Draw(options);
776}
777
778////////////////////////////////////////////////////////////////////////////////
779
781{
782 if (fWeightHist == NULL) {
783 const MarkovChain* markovChain = fInterval->GetChain();
784 // find the max weight value
785 Double_t maxWeight = 0;
786 Int_t size = markovChain->Size();
787 for (Int_t i = 0; i < size; i++)
788 if (markovChain->Weight(i) > maxWeight)
789 maxWeight = markovChain->Weight(i);
790 fWeightHist = new TH1F("mcmc_weight_hist", "MCMC Weight Histogram",
791 (Int_t)(maxWeight + 1), 0, maxWeight * 1.02);
792 for (Int_t i = 0; i < size; i++)
793 fWeightHist->Fill(markovChain->Weight(i));
794 }
795 fWeightHist->Draw(options);
796}
797
798/*
799/////////////////////////////////////////////////////////////////////
800 // 3-d plot of the parameter points
801 dataCanvas->cd(2);
802 // also plot the points in the markov chain
803 RooDataSet* markovChainData = ((MCMCInterval*)mcmcint)->GetChainAsDataSet();
804
805 TTree& chain = ((RooTreeDataStore*) markovChainData->store())->tree();
806 chain.SetMarkerStyle(6);
807 chain.SetMarkerColor(kRed);
808 chain.Draw("s:ratioSigEff:ratioBkgEff","","box"); // 3-d box proportional to posterior
809
810 // the points used in the profile construction
811 TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
812 parameterScan.SetMarkerStyle(24);
813 parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
814
815 chain.SetMarkerStyle(6);
816 chain.SetMarkerColor(kRed);
817 //chain.Draw("s:ratioSigEff:ratioBkgEff", "_MarkovChain_local_nll","box");
818 //chain.Draw("_MarkovChain_local_nll");
819////////////////////////////////////////////////////////////////////
820*/
#define coutE(a)
Definition: RooMsgService.h:34
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
@ kGray
Definition: Rtypes.h:63
@ kPink
Definition: Rtypes.h:65
@ kBlack
Definition: Rtypes.h:63
@ kGreen
Definition: Rtypes.h:64
@ kViolet
Definition: Rtypes.h:65
char * Form(const char *fmt,...)
RooAbsArg * first() const
virtual 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
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition: RooAbsPdf.h:118
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Int_t numBins(const char *rangeName=0) const
virtual Int_t getBins(const char *name=0) const
Get number of bins of currently defined range.
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
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.
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...
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:472
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition: RooPlot.cxx:1258
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1277
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void DrawWeightHist(const Option_t *options=NULL)
void DrawKeysPdfInterval(const Option_t *options=NULL)
void DrawInterval(const Option_t *options=NULL)
void * DrawPosteriorHist(const Option_t *options=NULL, const char *title=NULL, Bool_t scale=kTRUE)
void DrawChainScatter(RooRealVar &xVar, RooRealVar &yVar)
void * DrawPosteriorKeysProduct(const Option_t *options=NULL)
void SetMCMCInterval(MCMCInterval &interval)
void * DrawPosteriorKeysPdf(const Option_t *options=NULL)
void DrawTailFractionInterval(const Option_t *options=NULL)
void DrawPosterior(const Option_t *options=NULL)
void DrawHistInterval(const Option_t *options=NULL)
void DrawShortestInterval(const Option_t *options=NULL)
void DrawParameterVsTime(RooRealVar &param)
virtual ~MCMCIntervalPlot()
Destructor of SamplingDistribution.
RooNDKeysPdf * fPosteriorKeysPdf
void DrawNLLHist(const Option_t *options=NULL)
void Draw(const Option_t *options=NULL)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual Double_t UpperLimitByKeys(RooRealVar &param)
determine upper limit in the shortest interval by using keys pdf
virtual Bool_t GetUseKeys()
get whether we used kernel estimation to determine the interval
Definition: MCMCInterval.h:170
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
Definition: MCMCInterval.h:244
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
virtual Double_t LowerLimitByKeys(RooRealVar &param)
determine lower limit in the shortest interval by using keys pdf
virtual Double_t GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
virtual Double_t UpperLimitTailFraction(RooRealVar &param)
determine upper limit of the lower confidence interval
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList
Definition: MCMCInterval.h:99
virtual Int_t GetNumBurnInSteps()
get the number of steps in the chain to discard as burn-in,
Definition: MCMCInterval.h:176
virtual RooRealVar * GetNLLVar() const
Get a clone of the NLL variable from the markov chain.
Definition: MCMCInterval.h:218
virtual const MarkovChain * GetChain()
Get the markov chain on which this interval is based You do not own the returned MarkovChain*.
Definition: MCMCInterval.h:195
virtual Double_t GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual Double_t LowerLimitTailFraction(RooRealVar &param)
determine lower limit of the lower confidence interval
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
virtual Double_t UpperLimitByHist(RooRealVar &param)
determine upper limit using histogram
virtual Double_t LowerLimitByHist(RooRealVar &param)
determine lower limit using histogram
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual Int_t GetDimension() const
Get the number of parameters of interest in this interval.
Definition: MCMCInterval.h:191
Double_t GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
virtual RooArgSet * GetParameters() const
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
virtual Double_t 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 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:717
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2312
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:753
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1619
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1629
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6333
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8585
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2665
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
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:8650
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
TAxis * GetYaxis()
Definition: TH1.h:317
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition: TH1.cxx:7935
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:8666
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition: TH1.cxx:8024
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4899
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6234
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8434
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
A simple line.
Definition: TLine.h:23
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
Basic string class.
Definition: TString.h:131
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:418
const char * Data() const
Definition: TString.h:364
TString & Append(const char *cs)
Definition: TString.h:559
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
RooCmdArg Scaling(Bool_t flag)
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
RooCmdArg DrawOption(const char *opt)
RooCmdArg FillColor(Color_t color)
@ InputArguments
Definition: RooGlobalFunc.h:68
RooCmdArg MoveToBack()
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
RooCmdArg Normalization(Double_t scaleFactor)
RooCmdArg VLines()
Namespace for the RooStats classes.
Definition: Asimov.h:20
Definition: first.py:1