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 "TList.h"
38#include "TGraph.h"
39#include "RooRealVar.h"
40#include "RooPlot.h"
41#include "TH2.h"
42#include "TH1F.h"
43#include "RooArgList.h"
44#include "TAxis.h"
45#include "RooGlobalFunc.h"
46
47#include <iostream>
48
49// Extra draw commands
50//static const char* POSTERIOR_HIST = "posterior_hist";
51//static const char* POSTERIOR_KEYS_PDF = "posterior_keys_pdf";
52//static const char* POSTERIOR_KEYS_PRODUCT = "posterior_keys_product";
53//static const char* HIST_INTERVAL = "hist_interval";
54//static const char* KEYS_PDF_INTERVAL = "keys_pdf_interval";
55//static const char* TAIL_FRACTION_INTERVAL = "tail_fraction_interval";
56//static const char* OPTION_SEP = ":";
57
59
60using namespace std;
61using namespace RooStats;
62
63////////////////////////////////////////////////////////////////////////////////
64
66{
67 fInterval = NULL;
68 fParameters = NULL;
69 fPosteriorHist = NULL;
70 fPosteriorKeysPdf = NULL;
72 fDimension = 0;
75 fLineWidth = 1;
76 //fContourColor = kBlack;
78 fWalk = NULL;
79 fBurnIn = NULL;
80 fFirst = NULL;
81 fParamGraph = NULL;
82 fNLLGraph = NULL;
83 fNLLHist = NULL;
84 fWeightHist = NULL;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
92{
93 SetMCMCInterval(interval);
94 fPosteriorHist = NULL;
95 fPosteriorKeysPdf = NULL;
99 fLineWidth = 1;
100 //fContourColor = kBlack;
102 fWalk = NULL;
103 fBurnIn = NULL;
104 fFirst = NULL;
105 fParamGraph = NULL;
106 fNLLGraph = NULL;
107 fNLLHist = NULL;
108 fWeightHist = NULL;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114
116{
117 delete fParameters;
118 // kbelasco: why does deleting fPosteriorHist remove the graphics
119 // but deleting TGraphs doesn't?
120 //delete fPosteriorHist;
121 // can we delete fNLLHist and fWeightHist?
122 //delete fNLLHist;
123 //delete fWeightHist;
124
125 // kbelasco: should we delete fPosteriorKeysPdf and fPosteriorKeysProduct?
126 delete fPosteriorKeysPdf;
128
129 delete fWalk;
130 delete fBurnIn;
131 delete fFirst;
132 delete fParamGraph;
133 delete fNLLGraph;
134}
135
136////////////////////////////////////////////////////////////////////////////////
137
139{
140 fInterval = &interval;
143}
144
145////////////////////////////////////////////////////////////////////////////////
146
148{
149 DrawInterval(options);
150}
151
152////////////////////////////////////////////////////////////////////////////////
153
155{
156 if (fInterval->GetUseKeys())
157 DrawPosteriorKeysPdf(options);
158 else
159 DrawPosteriorHist(options);
160}
161
162////////////////////////////////////////////////////////////////////////////////
163
165 const char* title, Bool_t scale)
166{
167 if (fPosteriorHist == NULL)
169
170 if (fPosteriorHist == NULL) {
171 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorHist: "
172 << "Couldn't get posterior histogram." << endl;
173 return NULL;
174 }
175
176 // kbelasco: annoying hack because histogram drawing fails when it sees
177 // an un-recognized option like POSTERIOR_HIST, etc.
178 //const Option_t* myOpt = NULL;
179
180 //TString tmpOpt(options);
181 //if (tmpOpt.Contains("same"))
182 // myOpt = "same";
183
184 // scale so highest bin has height 1
185 if (scale)
188
189 TString ourTitle(GetTitle());
190 if (ourTitle.CompareTo("") == 0) {
191 if (title)
192 fPosteriorHist->SetTitle(title);
193 } else
195
196 //fPosteriorHist->Draw(myOpt);
197
198 return (void*)fPosteriorHist;
199}
200
201////////////////////////////////////////////////////////////////////////////////
202
204{
205 if (fPosteriorKeysPdf == NULL)
207
208 if (fPosteriorKeysPdf == NULL) {
209 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
210 << "Couldn't get posterior Keys PDF." << endl;
211 return NULL;
212 }
213
214 TString title(GetTitle());
215 Bool_t isEmpty = (title.CompareTo("") == 0);
216
217 if (fDimension == 1) {
219 RooPlot* frame = v->frame();
220 if (frame == NULL) {
221 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysPdf: "
222 << "Invalid parameter" << endl;
223 return NULL;
224 }
225 if (isEmpty)
226 frame->SetTitle(Form("Posterior Keys PDF for %s", v->GetName()));
227 else
228 frame->SetTitle(GetTitle());
229 //fPosteriorKeysPdf->plotOn(frame);
230 //fPosteriorKeysPdf->plotOn(frame,
231 // RooFit::Normalization(1, RooAbsReal::Raw));
232 //frame->Draw(options);
233 return (void*)frame;
234 } else if (fDimension == 2) {
235 RooArgList* axes = fInterval->GetAxes();
236 RooRealVar* xVar = (RooRealVar*)axes->at(0);
237 RooRealVar* yVar = (RooRealVar*)axes->at(1);
239 "keysPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
240 if (isEmpty)
241 keysHist->SetTitle(
242 Form("MCMC histogram of posterior Keys PDF for %s, %s",
243 axes->at(0)->GetName(), axes->at(1)->GetName()));
244 else
245 keysHist->SetTitle(GetTitle());
246
247 keysHist->Draw(options);
248 delete axes;
249 return NULL;
250 }
251 return NULL;
252}
253
254////////////////////////////////////////////////////////////////////////////////
255
257{
258 switch (fInterval->GetIntervalType()) {
260 DrawShortestInterval(options);
261 break;
264 break;
265 default:
266 coutE(InputArguments) << "MCMCIntervalPlot::DrawInterval(): " <<
267 "Interval type not supported" << endl;
268 break;
269 }
270}
271
272////////////////////////////////////////////////////////////////////////////////
273
275{
276 if (fInterval->GetUseKeys())
277 DrawKeysPdfInterval(options);
278 else
279 DrawHistInterval(options);
280}
281
282////////////////////////////////////////////////////////////////////////////////
283
285{
286 TString title(GetTitle());
287 Bool_t isEmpty = (title.CompareTo("") == 0);
288
289 if (fDimension == 1) {
290 // Draw the posterior keys PDF as well so the user can see where the
291 // limit bars line up
292 // fDimension == 1, so we know we will receive a RooPlot
293 RooPlot* frame = (RooPlot*)DrawPosteriorKeysPdf(options);
294
295 //Double_t height = 1;
296 //Double_t height = 2.0 * fInterval->GetKeysPdfCutoff();
297 Double_t height = fInterval->GetKeysMax();
298
302
303 if (frame != NULL && fPosteriorKeysPdf != NULL) {
304 // draw shading in interval
305 if (isEmpty)
306 frame->SetTitle(NULL);
307 else
308 frame->SetTitle(GetTitle());
309 frame->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
310 p->GetName()));
313 RooFit::Range(ll, ul, kFALSE),
318
319 // hack - this is drawn twice now:
320 // once by DrawPosteriorKeysPdf (which also configures things and sets
321 // the title), and once again here so the shading shows up behind.
324 }
325 if (frame) {
326 frame->Draw(options);
327 }
328
329 TLine* llLine = new TLine(ll, 0, ll, height);
330 TLine* ulLine = new TLine(ul, 0, ul, height);
331 llLine->SetLineColor(fLineColor);
332 ulLine->SetLineColor(fLineColor);
333 llLine->SetLineWidth(fLineWidth);
334 ulLine->SetLineWidth(fLineWidth);
335 llLine->Draw(options);
336 ulLine->Draw(options);
337 } else if (fDimension == 2) {
338 if (fPosteriorKeysPdf == NULL)
340
341 if (fPosteriorKeysPdf == NULL) {
342 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
343 << "Couldn't get posterior Keys PDF." << endl;
344 return;
345 }
346
347 RooArgList* axes = fInterval->GetAxes();
348 RooRealVar* xVar = (RooRealVar*)axes->at(0);
349 RooRealVar* yVar = (RooRealVar*)axes->at(1);
351 "keysContour2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
352 //if (isEmpty)
353 // contHist->SetTitle(Form("MCMC Keys conf. interval for %s, %s",
354 // axes->at(0)->GetName(), axes->at(1)->GetName()));
355 //else
356 // contHist->SetTitle(GetTitle());
357 if (!isEmpty)
358 contHist->SetTitle(GetTitle());
359 else
360 contHist->SetTitle(NULL);
361
362 contHist->SetStats(kFALSE);
363
364 TString tmpOpt(options);
365 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
366
368 contHist->SetContour(1, &cutoff);
369 contHist->SetLineColor(fLineColor);
370 contHist->SetLineWidth(fLineWidth);
371 contHist->Draw(tmpOpt.Data());
372 delete axes;
373 } else {
374 coutE(InputArguments) << "MCMCIntervalPlot::DrawKeysPdfInterval: "
375 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
376 }
377}
378
379////////////////////////////////////////////////////////////////////////////////
380
382{
383 TString title(GetTitle());
384 Bool_t isEmpty = (title.CompareTo("") == 0);
385
386 if (fDimension == 1) {
387 // draw lower and upper limits
391
392 // Draw the posterior histogram as well so the user can see where the
393 // limit bars line up
394 // fDimension == 1, so we know will get a TH1F*
395 TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
396 if (hist == NULL) return;
397 if (isEmpty)
398 hist->SetTitle(NULL);
399 else
400 hist->SetTitle(GetTitle());
401 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
402 p->GetName()));
403 hist->SetStats(kFALSE);
404 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
405 Double_t histCutoff = fInterval->GetHistCutoff();
406
407 Int_t i;
408 Int_t nBins = copy->GetNbinsX();
409 Double_t height;
410 for (i = 1; i <= nBins; i++) {
411 // remove bins with height < cutoff
412 height = copy->GetBinContent(i);
413 if (height < histCutoff) {
414 copy->SetBinContent(i, 0);
415 copy->SetBinError(i, 0);
416 }
417 }
418
419 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
420 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
421
422 copy->SetFillStyle(1001);
424 hist->Draw(options);
425 // to show the interval area
426 copy->Draw("HIST SAME");
427
429
430 TLine* llLine = new TLine(ll, 0, ll, 1);
431 TLine* ulLine = new TLine(ul, 0, ul, 1);
432 llLine->SetLineColor(fLineColor);
433 ulLine->SetLineColor(fLineColor);
434 llLine->SetLineWidth(fLineWidth);
435 ulLine->SetLineWidth(fLineWidth);
436 llLine->Draw(options);
437 ulLine->Draw(options);
438
439 } else if (fDimension == 2) {
440 if (fPosteriorHist == NULL)
442
443 if (fPosteriorHist == NULL) {
444 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
445 << "Couldn't get posterior histogram." << endl;
446 return;
447 }
448
449 RooArgList* axes = fInterval->GetAxes();
450 //if (isEmpty)
451 // fPosteriorHist->SetTitle(
452 // Form("MCMC histogram conf. interval for %s, %s",
453 // axes->at(0)->GetName(), axes->at(1)->GetName()));
454 //else
455 // fPosteriorHist->SetTitle(GetTitle());
456 if (!isEmpty)
458 else
460 delete axes;
461
463
464 TString tmpOpt(options);
465 if (!tmpOpt.Contains("CONT2")) tmpOpt.Append("CONT2");
466
467 Double_t cutoff = fInterval->GetHistCutoff();
468 fPosteriorHist->SetContour(1, &cutoff);
471 fPosteriorHist->Draw(tmpOpt.Data());
472 } else {
473 coutE(InputArguments) << "MCMCIntervalPlot::DrawHistInterval: "
474 << " Sorry: " << fDimension << "-D plots not currently supported" << endl;
475 }
476}
477
478////////////////////////////////////////////////////////////////////////////////
479
481{
482 TString title(GetTitle());
483 Bool_t isEmpty = (title.CompareTo("") == 0);
484
485 if (fDimension == 1) {
486 // Draw the posterior histogram as well so the user can see where the
487 // limit bars line up
491
492 TH1F* hist = (TH1F*)DrawPosteriorHist(options, NULL, false);
493 if (hist == NULL) return;
494 if (isEmpty)
495 hist->SetTitle(NULL);
496 else
497 hist->SetTitle(GetTitle());
498 hist->GetYaxis()->SetTitle(Form("Posterior for parameter %s",
499 p->GetName()));
500 hist->SetStats(kFALSE);
501 TH1F* copy = (TH1F*)hist->Clone(Form("%s_copy", hist->GetTitle()));
502
503 Int_t i;
504 Int_t nBins = copy->GetNbinsX();
505 Double_t center;
506 for (i = 1; i <= nBins; i++) {
507 // remove bins outside interval
508 center = copy->GetBinCenter(i);
509 if (center < ll || center > ul) {
510 copy->SetBinContent(i, 0);
511 copy->SetBinError(i, 0);
512 }
513 }
514
515 hist->Scale(1/hist->GetBinContent(hist->GetMaximumBin()));
516 copy->Scale(1/copy->GetBinContent(hist->GetMaximumBin()));
517
518 copy->SetFillStyle(1001);
520 hist->Draw(options);
521 copy->Draw("hist same");
522
523 // draw lower and upper limits
524 TLine* llLine = new TLine(ll, 0, ll, 1);
525 TLine* ulLine = new TLine(ul, 0, ul, 1);
526 llLine->SetLineColor(fLineColor);
527 ulLine->SetLineColor(fLineColor);
528 llLine->SetLineWidth(fLineWidth);
529 ulLine->SetLineWidth(fLineWidth);
530 llLine->Draw(options);
531 ulLine->Draw(options);
532 } else {
533 coutE(InputArguments) << "MCMCIntervalPlot::DrawTailFractionInterval: "
534 << " Sorry: " << fDimension << "-D plots not currently supported"
535 << endl;
536 }
537}
538
539////////////////////////////////////////////////////////////////////////////////
540
542{
543 if (fPosteriorKeysProduct == NULL)
545
546 if (fPosteriorKeysProduct == NULL) {
547 coutE(InputArguments) << "MCMCIntervalPlot::DrawPosteriorKeysProduct: "
548 << "Couldn't get posterior Keys product." << endl;
549 return NULL;
550 }
551
552 RooArgList* axes = fInterval->GetAxes();
553
554 TString title(GetTitle());
555 Bool_t isEmpty = (title.CompareTo("") == 0);
556
557 if (fDimension == 1) {
558 RooPlot* frame = ((RooRealVar*)fParameters->first())->frame();
559 if (!frame) return NULL;
560 if (isEmpty)
561 frame->SetTitle(Form("Posterior Keys PDF * Heaviside product for %s",
562 axes->at(0)->GetName()));
563 else
564 frame->SetTitle(GetTitle());
565 //fPosteriorKeysProduct->plotOn(frame);
568 frame->Draw(options);
569 return (void*)frame;
570 } else if (fDimension == 2) {
571 RooRealVar* xVar = (RooRealVar*)axes->at(0);
572 RooRealVar* yVar = (RooRealVar*)axes->at(1);
574 "prodPlot2D", *xVar, RooFit::YVar(*yVar), RooFit::Scaling(kFALSE));
575 if (isEmpty)
576 productHist->SetTitle(
577 Form("MCMC Posterior Keys Product Hist. for %s, %s",
578 axes->at(0)->GetName(), axes->at(1)->GetName()));
579 else
580 productHist->SetTitle(GetTitle());
581 productHist->Draw(options);
582 return NULL;
583 }
584 delete axes;
585 return NULL;
586}
587
588////////////////////////////////////////////////////////////////////////////////
589
591{
592 const MarkovChain* markovChain = fInterval->GetChain();
593
594 Int_t size = markovChain->Size();
595 Int_t burnInSteps;
596 if (fShowBurnIn)
597 burnInSteps = fInterval->GetNumBurnInSteps();
598 else
599 burnInSteps = 0;
600
601 Double_t* x = new Double_t[size - burnInSteps];
602 Double_t* y = new Double_t[size - burnInSteps];
603 Double_t* burnInX = NULL;
604 Double_t* burnInY = NULL;
605 if (burnInSteps > 0) {
606 burnInX = new Double_t[burnInSteps];
607 burnInY = new Double_t[burnInSteps];
608 }
609 Double_t firstX;
610 Double_t firstY;
611
612 for (Int_t i = burnInSteps; i < size; i++) {
613 x[i - burnInSteps] = markovChain->Get(i)->getRealValue(xVar.GetName());
614 y[i - burnInSteps] = markovChain->Get(i)->getRealValue(yVar.GetName());
615 }
616
617 for (Int_t i = 0; i < burnInSteps; i++) {
618 burnInX[i] = markovChain->Get(i)->getRealValue(xVar.GetName());
619 burnInY[i] = markovChain->Get(i)->getRealValue(yVar.GetName());
620 }
621
622 firstX = markovChain->Get(0)->getRealValue(xVar.GetName());
623 firstY = markovChain->Get(0)->getRealValue(yVar.GetName());
624
625 TString title(GetTitle());
626 Bool_t isEmpty = (title.CompareTo("") == 0);
627
628 TGraph* walk = new TGraph(size - burnInSteps, x, y);
629 if (isEmpty)
630 walk->SetTitle(Form("2-D Scatter Plot of Markov chain for %s, %s",
631 xVar.GetName(), yVar.GetName()));
632 else
633 walk->SetTitle(GetTitle());
634 // kbelasco: figure out how to set TGraph variable ranges
635 walk->GetXaxis()->Set(xVar.numBins(), xVar.getMin(), xVar.getMax());
636 walk->GetXaxis()->SetTitle(xVar.GetName());
637 walk->GetYaxis()->Set(yVar.numBins(), yVar.getMin(), yVar.getMax());
638 walk->GetYaxis()->SetTitle(yVar.GetName());
639 walk->SetLineColor(kGray);
640 walk->SetMarkerStyle(6);
641 walk->SetMarkerColor(kViolet);
642 walk->Draw("A,L,P,same");
643
644 TGraph* burnIn = NULL;
645 if (burnInX != NULL && burnInY != NULL) {
646 burnIn = new TGraph(burnInSteps - 1, burnInX, burnInY);
647 burnIn->SetLineColor(kPink);
648 burnIn->SetMarkerStyle(6);
649 burnIn->SetMarkerColor(kPink);
650 burnIn->Draw("L,P,same");
651 }
652
653 TGraph* first = new TGraph(1, &firstX, &firstY);
654 first->SetLineColor(kGreen);
655 first->SetMarkerStyle(3);
656 first->SetMarkerSize(2);
657 first->SetMarkerColor(kGreen);
658 first->Draw("L,P,same");
659
660 //walkCanvas->Update();
661 delete [] x;
662 delete [] y;
663 if (burnInX != NULL) delete [] burnInX;
664 if (burnInY != NULL) delete [] burnInY;
665 //delete walk;
666 //delete burnIn;
667 //delete first;
668}
669
670////////////////////////////////////////////////////////////////////////////////
671
673{
674 const MarkovChain* markovChain = fInterval->GetChain();
675 Int_t size = markovChain->Size();
676 Int_t numEntries = 2 * size;
677 Double_t* value = new Double_t[numEntries];
678 Double_t* time = new Double_t[numEntries];
679 Double_t val;
680 Int_t weight;
681 Int_t t = 0;
682 for (Int_t i = 0; i < size; i++) {
683 val = markovChain->Get(i)->getRealValue(param.GetName());
684 weight = (Int_t)markovChain->Weight();
685 value[2*i] = val;
686 value[2*i + 1] = val;
687 time[2*i] = t;
688 t += weight;
689 time[2*i + 1] = t;
690 }
691
692 TString title(GetTitle());
693 Bool_t isEmpty = (title.CompareTo("") == 0);
694
695 TGraph* paramGraph = new TGraph(numEntries, time, value);
696 if (isEmpty)
697 paramGraph->SetTitle(Form("%s vs. time in Markov chain",param.GetName()));
698 else
699 paramGraph->SetTitle(GetTitle());
700 paramGraph->GetXaxis()->SetTitle("Time (discrete steps)");
701 paramGraph->GetYaxis()->SetTitle(param.GetName());
702 //paramGraph->SetLineColor(fLineColor);
703 paramGraph->Draw("A,L,same");
704 delete [] value;
705 delete [] time;
706 //gPad->Update();
707}
708
709////////////////////////////////////////////////////////////////////////////////
710
712{
713 const MarkovChain* markovChain = fInterval->GetChain();
714 Int_t size = markovChain->Size();
715 Int_t numEntries = 2 * size;
716 Double_t* nllValue = new Double_t[numEntries];
717 Double_t* time = new Double_t[numEntries];
718 Double_t nll;
719 Int_t weight;
720 Int_t t = 0;
721 for (Int_t i = 0; i < size; i++) {
722 nll = markovChain->NLL(i);
723 weight = (Int_t)markovChain->Weight();
724 nllValue[2*i] = nll;
725 nllValue[2*i + 1] = nll;
726 time[2*i] = t;
727 t += weight;
728 time[2*i + 1] = t;
729 }
730
731 TString title(GetTitle());
732 Bool_t isEmpty = (title.CompareTo("") == 0);
733
734 TGraph* nllGraph = new TGraph(numEntries, time, nllValue);
735 if (isEmpty)
736 nllGraph->SetTitle("NLL value vs. time in Markov chain");
737 else
738 nllGraph->SetTitle(GetTitle());
739 nllGraph->GetXaxis()->SetTitle("Time (discrete steps)");
740 nllGraph->GetYaxis()->SetTitle("NLL (-log(likelihood))");
741 //nllGraph->SetLineColor(fLineColor);
742 nllGraph->Draw("A,L,same");
743 //gPad->Update();
744 delete [] nllValue;
745 delete [] time;
746}
747
748////////////////////////////////////////////////////////////////////////////////
749
751{
752 if (fNLLHist == NULL) {
753 const MarkovChain* markovChain = fInterval->GetChain();
754 // find the max NLL value
755 Double_t maxNLL = 0;
756 Int_t size = markovChain->Size();
757 for (Int_t i = 0; i < size; i++)
758 if (markovChain->NLL(i) > maxNLL)
759 maxNLL = markovChain->NLL(i);
760 RooRealVar* nllVar = fInterval->GetNLLVar();
761 fNLLHist = new TH1F("mcmc_nll_hist", "MCMC NLL Histogram",
762 nllVar->getBins(), 0, maxNLL);
763 TString title(GetTitle());
764 Bool_t isEmpty = (title.CompareTo("") == 0);
765 if (!isEmpty)
767 fNLLHist->GetXaxis()->SetTitle("-log(likelihood)");
768 for (Int_t i = 0; i < size; i++)
769 fNLLHist->Fill(markovChain->NLL(i), markovChain->Weight());
770 }
771 fNLLHist->Draw(options);
772}
773
774////////////////////////////////////////////////////////////////////////////////
775
777{
778 if (fWeightHist == NULL) {
779 const MarkovChain* markovChain = fInterval->GetChain();
780 // find the max weight value
781 Double_t maxWeight = 0;
782 Int_t size = markovChain->Size();
783 for (Int_t i = 0; i < size; i++)
784 if (markovChain->Weight(i) > maxWeight)
785 maxWeight = markovChain->Weight(i);
786 fWeightHist = new TH1F("mcmc_weight_hist", "MCMC Weight Histogram",
787 (Int_t)(maxWeight + 1), 0, maxWeight * 1.02);
788 for (Int_t i = 0; i < size; i++)
789 fWeightHist->Fill(markovChain->Weight(i));
790 }
791 fWeightHist->Draw(options);
792}
793
794/*
795/////////////////////////////////////////////////////////////////////
796 // 3-d plot of the parameter points
797 dataCanvas->cd(2);
798 // also plot the points in the markov chain
799 RooDataSet* markovChainData = ((MCMCInterval*)mcmcint)->GetChainAsDataSet();
800
801 TTree& chain = ((RooTreeDataStore*) markovChainData->store())->tree();
802 chain.SetMarkerStyle(6);
803 chain.SetMarkerColor(kRed);
804 chain.Draw("s:ratioSigEff:ratioBkgEff","","box"); // 3-d box proportional to posterior
805
806 // the points used in the profile construction
807 TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
808 parameterScan.SetMarkerStyle(24);
809 parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
810
811 chain.SetMarkerStyle(6);
812 chain.SetMarkerColor(kRed);
813 //chain.Draw("s:ratioSigEff:ratioBkgEff", "_MarkovChain_local_nll","box");
814 //chain.Draw("_MarkovChain_local_nll");
815////////////////////////////////////////////////////////////////////
816*/
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutE(a)
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:101
const Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
@ 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
char * Form(const char *fmt,...)
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.
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:121
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.
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:44
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition RooPlot.cxx:1257
TAxis * GetYaxis() const
Definition RooPlot.cxx:1278
static RooPlot * frame(const RooAbsRealLValue &var, Double_t xmin, Double_t xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:249
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:661
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
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.
void DrawNLLHist(const Option_t *options=NULL)
void Draw(const Option_t *options=NULL)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
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
virtual enum IntervalType GetIntervalType()
Return the type of this interval.
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
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_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.
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:731
A TGraph is an 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:2353
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition TGraph.cxx:769
TAxis * GetXaxis() const
Get x axis of the graph.
Definition TGraph.cxx:1634
TAxis * GetYaxis() const
Get y axis of the graph.
Definition TGraph.cxx:1644
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:575
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition TH1.cxx:6667
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:8971
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition TH1.cxx:2741
virtual Int_t GetNbinsX() const
Definition TH1.h:296
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:9036
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3351
TAxis * GetYaxis()
Definition TH1.h:321
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition TH1.cxx:8313
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:9052
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074
virtual Int_t GetMaximumBin() const
Return location of bin with maximum value in the range.
Definition TH1.cxx:8407
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4994
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6553
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition TH1.cxx:8820
2-D histogram with a float per channel (see TH1 documentation)}
Definition TH2.h:251
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
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:267
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:442
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 Scaling(Bool_t flag)
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
RooCmdArg FillColor(Color_t color)
RooCmdArg DrawOption(const char *opt)
RooCmdArg MoveToBack()
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
RooCmdArg Normalization(Double_t scaleFactor)
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
Definition first.py:1