BayesianAnalysisToolkit  0.9.3
Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
BCSummaryTool Class Reference

A class for summarizing the results of an analysis. More...

#include <BCSummaryTool.h>

Collaboration diagram for BCSummaryTool:
Collaboration graph
[legend]

Public Member Functions

Constructors and destructors
 BCSummaryTool ()
 
 BCSummaryTool (BCModel *model)
 
 ~BCSummaryTool ()
 
Member functions (get)
BCSummaryPriorModel * GetPriorModel ()
 
Member functions (set)
void SetModel (BCModel *model)
 
Member functions (misc)
int CalculatePriorModel ()
 
int CopySummaryData ()
 
int PrintParameterPlot (const char *filename="parameters.pdf")
 
int PrintCorrelationMatrix (const char *filename="matrix.pdf")
 
int PrintCorrelationPlot (const char *filename="correlation.pdf")
 
int DrawKnowledgeUpdatePlot1D (int index, std::string options_post="", std::string options_prior="")
 
int PrintKnowledgeUpdatePlot1D (int index, const char *filename, std::string options_post="", std::string options_prior="")
 
int PrintKnowledgeUpdatePlots (const char *filename="update.pdf", std::string options="")
 
int PrintParameterLatex (const char *filename)
 

Static Private Member Functions

static unsigned int getNextIndex ()
 

Private Attributes

BCModel * fModel
 
std::vector< std::string > fParName
 
std::vector< double > fParMin
 
std::vector< double > fParMax
 
std::vector< double > fCorrCoeff
 
std::vector< double > fMargMode
 
std::vector< double > fMean
 
std::vector< double > fGlobalMode
 
std::vector< double > fQuantiles
 
std::vector< double > fSmallInt
 
std::vector< double > fRMS
 
std::vector< double > fSumProb
 
BCSummaryPriorModel * fPriorModel
 
bool fFlagInfoMarg
 
bool fFlagInfoOpt
 

Static Private Attributes

static unsigned int fHCounter =0
 

Detailed Description

A class for summarizing the results of an analysis.

This class can be used to summarize the results of an analysis. The prior and posterior probabilities are compared.

Author
Daniel Kollar
Kevin Kröninger
Version
1.0.0
Date
15.02.2010

Definition at line 34 of file BCSummaryTool.h.

Constructor & Destructor Documentation

BCSummaryTool::BCSummaryTool ( )

The default constructor.

Definition at line 42 of file BCSummaryTool.cxx.

43  : fModel(0)
44  , fPriorModel(0)
45  , fFlagInfoMarg(false)
46  , fFlagInfoOpt(false)
47 {
48  // define sum of probabilities for quantiles
49  fSumProb.push_back(0.05);
50  fSumProb.push_back(0.10);
51  fSumProb.push_back(0.1587);
52  fSumProb.push_back(0.50);
53  fSumProb.push_back(0.8413);
54  fSumProb.push_back(0.90);
55  fSumProb.push_back(0.95);
56 
57  // set text style
58  gStyle->SetPaintTextFormat(".2g");
59 }
BCSummaryTool::BCSummaryTool ( BCModel *  model)

A constructor.

Definition at line 62 of file BCSummaryTool.cxx.

63  : fModel(model)
64  , fPriorModel(0)
65  , fFlagInfoMarg(false)
66  , fFlagInfoOpt(false)
67 {
68  // define sum of probabilities for quantiles
69  fSumProb.push_back(0.05);
70  fSumProb.push_back(0.10);
71  fSumProb.push_back(0.1587);
72  fSumProb.push_back(0.50);
73  fSumProb.push_back(0.8413);
74  fSumProb.push_back(0.90);
75  fSumProb.push_back(0.95);
76 
77  // set text style
78  gStyle->SetPaintTextFormat(".2g");
79 }
BCSummaryTool::~BCSummaryTool ( )

The default destructor.

Definition at line 82 of file BCSummaryTool.cxx.

83 {
84  delete fPriorModel;
85 }

Member Function Documentation

int BCSummaryTool::CalculatePriorModel ( )

Calculate the marginalized distributions using the prior knowledge alone.

Returns
An error flag.

Definition at line 931 of file BCSummaryTool.cxx.

932 {
933  // create new prior model
934  delete fPriorModel;
935 
936  fPriorModel = new BCSummaryPriorModel();
937 
938  // set model
939  fPriorModel->SetModel(fModel);
940 
941  // perform marginalization
942  fPriorModel->MarginalizeAll();
943 
944  // perform minimization
945  fPriorModel->FindMode( fPriorModel->GetBestFitParameters() );
946 
947  // no error
948  return 1;
949 }
int BCSummaryTool::CopySummaryData ( )

Copy the summary information from the model.

Returns
An error flag.

Definition at line 88 of file BCSummaryTool.cxx.

89 {
90  // check if model exists
91  if (!fModel)
92  return 0;
93 
94  // clear information
95  fParName.clear();
96  fParMin.clear();
97  fParMax.clear();
98  fMean.clear();
99  fMargMode.clear();
100  fGlobalMode.clear();
101  fQuantiles.clear();
102  fSmallInt.clear();
103  fRMS.clear();
104  fCorrCoeff.clear();
105 
106  // copy information from marginalized distributions
107  for (unsigned i = 0; i < fModel->GetNParameters(); ++i) {
108 
109  // copy parameter information
110  fParName.push_back( (fModel->GetParameter(i)->GetLatexName()) );
111  fParMin.push_back( fModel->GetParameter(i)->GetLowerLimit() );
112  fParMax.push_back( fModel->GetParameter(i)->GetUpperLimit() );
113 
114  // copy 1D marginalized information
115  BCH1D * bch1d_temp = fModel->GetMarginalized(i);
116  if (bch1d_temp) {
117  fFlagInfoMarg = true;
118  fMean.push_back( bch1d_temp->GetMean() );
119  fRMS.push_back( bch1d_temp->GetRMS() );
120  fMargMode.push_back( bch1d_temp->GetMode() );
121  for (unsigned j = 0; j < fSumProb.size(); ++j)
122  fQuantiles.push_back( bch1d_temp->GetQuantile( fSumProb.at(j) ) );
123  std::vector<double> intervals = bch1d_temp->GetSmallestIntervals();
124  int nintervals = int(intervals.size() / 5);
125  fSmallInt.push_back(nintervals);
126  fSmallInt.insert( fSmallInt.end(), intervals.begin(), intervals.end() );
127  }
128  else {
129  double tmpval = fModel->GetParameter(i)->GetUpperLimit() - fModel->GetParameter(i)->GetLowerLimit();
130  tmpval = fModel->GetParameter(i)->GetLowerLimit() - 2. * tmpval;
131  fMean.push_back( tmpval );
132  fRMS.push_back( tmpval );
133  fMargMode.push_back( tmpval );
134  for (unsigned j = 0; j < fSumProb.size(); ++j)
135  fQuantiles.push_back( tmpval );
136  fSmallInt.push_back( 0 );
137  }
138 
139  // copy 2D marginal information
140  for (unsigned j = 0; j < fModel->GetNParameters(); ++j) {
141  if (i == j)
142  fCorrCoeff.push_back(1.);
143  else {
144  BCH2D * bch2d_temp = fModel->GetMarginalized(i, j);
145  if ( bch2d_temp ) {
146  fFlagInfoMarg = true;
147  fCorrCoeff.push_back( bch2d_temp->GetHistogram()->GetCorrelationFactor() );
148  }
149  else
150  fCorrCoeff.push_back( 0. );
151  }
152  }
153 
154  // copy optimization information
155  if ((fModel->GetBestFitParameters()).size() > 0) {
156  fFlagInfoOpt = true;
157  fGlobalMode.push_back ( (fModel->GetBestFitParameters()).at(i) );
158  }
159  }
160 
161  // no error
162  return 1;
163 }
int BCSummaryTool::DrawKnowledgeUpdatePlot1D ( int  index,
std::string  options_post = "",
std::string  options_prior = "" 
)

Draw a comparison of the prior knowledge to the posterior knowledge for each parameter.

Returns
An error flag.

Definition at line 583 of file BCSummaryTool.cxx.

584 {
585  // option flags
586  bool flag_slice_post = false;
587  bool flag_slice_prior = false;
588 
589  // check content of options string
590  if (options_post.find("slice") < options_post.size()) {
591  flag_slice_post = true;
592  }
593  if (options_prior.find("slice") < options_prior.size()) {
594  flag_slice_prior = true;
595  }
596 
597  // prepare legend
598  TLegend* legend = new TLegend();
599  legend->SetBorderSize(0);
600  legend->SetFillColor(kWhite);
601  legend->SetTextAlign(12);
602  legend->SetTextFont(62);
603  legend->SetTextSize(0.03);
604 
605  // get histograms;
606  const BCParameter * par = fModel->GetParameter(index);
607  BCH1D* hist_prior = fPriorModel->GetMarginalized(par);
608  BCH1D* hist_posterior = 0;
609 
610  if (flag_slice_prior && fPriorModel->GetNParameters()==2) {
611  if (index == 0) {
612  TH1D* hist = fPriorModel->GetSlice(fPriorModel->GetParameter(0),fPriorModel->GetParameter(1))->GetHistogram()->ProjectionX(Form("projx_%i",BCLog::GetHIndex()));
613  hist->Scale(1.0/hist->Integral("width"));
614  for (int i = 1; i <= hist_prior->GetHistogram()->GetNbinsX(); ++i)
615  hist_prior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
616  }
617  else {
618  TH1D* hist = fPriorModel->GetSlice(fPriorModel->GetParameter(0),fPriorModel->GetParameter(1))->GetHistogram()->ProjectionY(Form("projy_%i",BCLog::GetHIndex()));
619  hist->Scale(1.0/hist->Integral("width"));
620  for (int i = 1; i <= hist_prior->GetHistogram()->GetNbinsX(); ++i)
621  hist_prior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
622  }
623  hist_prior->GetHistogram()->SetStats(kFALSE);
624  }
625  else if (flag_slice_prior && fPriorModel->GetNParameters()==1) {
626  hist_prior = fPriorModel->GetSlice(par);
627  hist_prior->GetHistogram()->SetStats(kFALSE);
628  }
629  // if marginal doesn't exist, skip ahead
630  if ( !hist_prior)
631  return 0;
632 
633  hist_prior->GetHistogram()->SetLineColor(kRed);
634 
635  hist_posterior = fModel->GetMarginalized(par);
636  if (flag_slice_post && fModel->GetNParameters()==2) {
637  if (index == 0) {
638  TH1D* hist = fModel->GetSlice(fModel->GetParameter(0),fModel->GetParameter(1))->GetHistogram()->ProjectionX(Form("projx_%i",BCLog::GetHIndex()));
639  hist->Scale(1.0/hist->Integral("width"));
640  for (int i = 1; i <= hist_posterior->GetHistogram()->GetNbinsX(); ++i)
641  hist_posterior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
642  }
643  else {
644  TH1D* hist = fModel->GetSlice(fModel->GetParameter(0),fModel->GetParameter(1))->GetHistogram()->ProjectionY(Form("projy_%i",BCLog::GetHIndex()));
645  hist->Scale(1.0/hist->Integral("width"));
646  for (int i = 1; i <= hist_posterior->GetHistogram()->GetNbinsX(); ++i)
647  hist_posterior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
648  }
649  hist_posterior->GetHistogram()->SetStats(kFALSE);
650  }
651  else if (flag_slice_post && fModel->GetNParameters()==1) {
652  hist_posterior = fModel->GetSlice(par);
653  hist_posterior->GetHistogram()->SetStats(kFALSE);
654  }
655 
656  // if marginal doesn't exist, skip ahead
657  if ( !hist_posterior)
658  return 0;
659 
660  legend->AddEntry(hist_prior->GetHistogram(), "prior", "L");
661  legend->AddEntry(hist_posterior->GetHistogram(), "posterior", "L");
662 
663  // scale histograms
664  hist_posterior->GetHistogram()->Scale(1./hist_posterior->GetHistogram()->Integral("width"));
665  hist_prior->GetHistogram()->Scale(1.0/hist_prior->GetHistogram()->Integral("width"));
666 
667  // get maximum
668  double max_prior = hist_prior->GetHistogram()->GetMaximum();
669  double max_posterior = hist_posterior->GetHistogram()->GetMaximum();
670  double maxy = 1.1 * TMath::Max(max_prior, max_posterior);
671 
672  double height = 0.03*legend->GetNRows();
673 
674  // plot
675  hist_prior->GetHistogram()->GetXaxis()->SetNdivisions(508);
676  hist_posterior->GetHistogram()->GetXaxis()->SetNdivisions(508);
677 
678  hist_prior->Draw(options_prior);
679  hist_posterior->Draw(std::string(options_post+"same").c_str());
680 
681  // scale axes
682  hist_prior->GetHistogram()->GetYaxis()->SetRangeUser(0.0, maxy);
683  hist_posterior->GetHistogram()->GetYaxis()->SetRangeUser(0.0, maxy);
684 
685  gPad->SetTopMargin(0.02);
686  double xlegend1 = gPad->GetLeftMargin() + 0.10 * (1.0 - gPad->GetRightMargin() - gPad->GetLeftMargin());
687  double xlegend2 = 1.0-gPad->GetRightMargin();
688  double ylegend1 = 1.-gPad->GetTopMargin()-height;
689  double ylegend2 = 1.-gPad->GetTopMargin();
690 
691  // place legend on top of histogram
692  legend->SetX1NDC(xlegend1);
693  legend->SetX2NDC(xlegend2);
694  legend->SetY1NDC(ylegend1);
695  legend->SetY2NDC(ylegend2);
696 
697  // draw legend
698  legend->Draw();
699 
700  // rescale top margin
701  gPad->SetTopMargin(1.-ylegend1+0.01);
702 
703  gPad->RedrawAxis();
704 
705  return 1;
706 }
static unsigned int BCSummaryTool::getNextIndex ( )
inlinestaticprivate

Helper method to get an unique number to be used in histogram naming

Definition at line 131 of file BCSummaryTool.h.

132  { return ++fHCounter; }
BCSummaryPriorModel* BCSummaryTool::GetPriorModel ( )
inline

Retrieve pointer to the Prior model to allow for its detailed setup

Definition at line 59 of file BCSummaryTool.h.

60  { return fPriorModel; }
int BCSummaryTool::PrintCorrelationMatrix ( const char *  filename = "matrix.pdf")

Print a correlation matrix for the parameters.

Returns
An error flag.

Definition at line 317 of file BCSummaryTool.cxx.

318 {
319  // copy summary data
320  if (!CopySummaryData())
321  return 0;
322 
323  // check if marginalized information is there
324  if (!fFlagInfoMarg)
325  return 0;
326 
327  // get number of parameters
328  int npar = fModel->GetNParameters();
329 
330  // create histogram
331  TH2D * hist_corr = new TH2D(
332  TString::Format("hist_corr_%d",getNextIndex()),
333  ";;",npar, -0.5, npar-0.5,npar, -0.5, npar-0.5);
334  hist_corr->SetStats(kFALSE);
335  hist_corr->GetXaxis()->SetTickLength(0.0);
336  hist_corr->GetYaxis()->SetTickLength(0.0);
337  hist_corr->GetZaxis()->SetRangeUser(-1.0, 1.0);
338 
339  for (int i = 0; i < npar; ++i) {
340  hist_corr->GetXaxis()->SetLabelSize(0.06);
341  hist_corr->GetYaxis()->SetLabelSize(0.06);
342  if (npar < 5) {
343  hist_corr->GetXaxis()->SetBinLabel( i+1, fParName.at(i).c_str() );
344  hist_corr->GetYaxis()->SetBinLabel( npar-i, fParName.at(i).c_str() );
345  }
346  else {
347  hist_corr->GetXaxis()->SetBinLabel( i+1, TString::Format("%d",i) );
348  hist_corr->GetYaxis()->SetBinLabel( npar-i, TString::Format("%d",i) );
349  }
350  }
351 
352  // fill plot
353  for (int i = 0; i < npar; ++i)
354  for (int j = 0; j < npar; ++j) {
355  int index = i * npar + j;
356  double corr = fCorrCoeff.at(index);
357  hist_corr->SetBinContent(i+1, npar-j, corr);
358  }
359 
360  // print to file
361  TCanvas * c_corr = new TCanvas(TString::Format("c_corr_matrix_%d",getNextIndex()));
362  c_corr->cd();
363  hist_corr->Draw("colz text");
364 
365  TF1 * f = new TF1("fUp","x",-0.5,npar-0.5);
366  TGaxis * A1 = new TGaxis(-0.5,npar-0.5,npar-0.5,npar-0.5,"fUp",100,"-");
367  A1->ImportAxisAttributes(hist_corr->GetXaxis());
368  A1->Draw();
369 
370  // redraw the histogram to overlay thetop axis tick marks since
371  // we don't know how to make them disappear
372  hist_corr->GetXaxis()->SetLabelSize(0.);
373  hist_corr->Draw("colz text same");
374 
375  for (int i = 0; i < npar; ++i)
376  for (int j = 0; j < npar; ++j) {
377  BCH2D * bch2d_temp = fModel->GetMarginalized(fModel->GetParameter(i),fModel->GetParameter(j));
378  if ( bch2d_temp || i==j )
379  continue;
380 
381  TBox * bempty = new TBox(
382  hist_corr->GetXaxis()->GetBinLowEdge(i+1),
383  hist_corr->GetYaxis()->GetBinLowEdge(npar-j),
384  hist_corr->GetXaxis()->GetBinLowEdge(i+2),
385  hist_corr->GetYaxis()->GetBinLowEdge(npar-j+1)
386  );
387  bempty->SetLineStyle(0);
388  bempty->SetLineWidth(0);
389  bempty->SetFillColor(kWhite);
390  bempty->Draw();
391  }
392 
393  // redraw top and right axes
394  TLine * lA1 = new TLine(-0.5,npar-0.5,npar-0.5,npar-0.5);
395  lA1->Draw("same");
396  TLine * lA2 = new TLine(npar-0.5,npar-0.5,npar-0.5,-0.5);
397  lA2->Draw("same");
398 
399  gPad->RedrawAxis();
400  c_corr->Print(filename);
401 
402  delete f;
403  delete A1;
404  delete lA1;
405  delete lA2;
406  delete hist_corr;
407  delete c_corr;
408 
409  // no error
410  return 1;
411 }
int BCSummaryTool::PrintCorrelationPlot ( const char *  filename = "correlation.pdf")

Print a correlation plot for the parameters.

Returns
An error flag.

Definition at line 414 of file BCSummaryTool.cxx.

415 {
416  // get number of parameters
417  int npar = fModel->GetNParameters();
418 
419  TCanvas * c = new TCanvas(TString::Format("c_corr_%d",getNextIndex()));
420  c->cd();
421 
422  double margin = .1;
423  double padsize = (1. - 2.*margin) / (double)npar;
424 
425  // array with pads holding the histograms
426  TPad ** pad = new TPad*[npar*npar];
427 
428  // position of pads
429  double xlow, xup, ylow, yup;
430  double marginleft, marginright, margintop, marginbottom;
431  marginleft=marginright=margintop=marginbottom=0.01;
432 
433  // rotate label to avoid overlap with long parameter names
434  size_t maxlength = 0;
435  for (int i = 0; i < npar ; ++i) {
436  maxlength = std::max(maxlength, fModel->GetParameter(i)->GetName().length());
437  }
438 
439  double rotation = 0;
440  short xalignment = 22;
441  short yalignment = 22;
442 
443  if (maxlength > 20) {
444  rotation = 10;
445  xalignment = 12;
446  yalignment = 32;
447  }
448 
449  // drawing all histograms
450  for (int i=npar-1;i>=0;--i) {
451  xlow = (double)i*padsize + margin;
452  xup = xlow+padsize;
453 
454  for (int j=i;j<=npar-1;j++) {
455  yup = 1. - (double)j*padsize - margin;
456  ylow = yup-padsize;
457 
458  // preparing the pad
459  int ipad=i*npar+j;
460  pad[ipad]=new TPad(TString::Format("pad_%d_%d",ipad,getNextIndex()), "", xlow, ylow, xup, yup);
461  pad[ipad]->SetTopMargin(margintop);
462  pad[ipad]->SetBottomMargin(marginbottom);
463  pad[ipad]->SetLeftMargin(marginleft);
464  pad[ipad]->SetRightMargin(marginright);
465  pad[ipad]->SetFillColor(kWhite);
466 
467  pad[ipad]->Draw();
468  pad[ipad]->cd();
469 
470  // get the histogram
471  TH1 * hh = 0;
472  BCH1D * bh1 = 0;
473  BCH2D * bh2 = 0;
474  if(i==j) {
475  bh1 = fModel->GetMarginalized(fModel->GetParameter(i));
476  if (bh1)
477  hh = (TH1D*)bh1->GetHistogram()->Clone();
478  }
479  else {
480  bh2 = fModel->GetMarginalized(fModel->GetParameter(i),fModel->GetParameter(j));
481  if (bh2)
482  hh = (TH2D*)bh2->GetHistogram()->Clone();
483  }
484 
485  // if the histogram is not available, draw N/A
486  if (!hh) {
487  pad[ipad]->SetFillColor(kGray);
488  TBox * box = new TBox(marginleft,marginbottom,1.-marginright,1.-margintop);
489  box->SetLineWidth(1);
490  box->SetLineColor(kGray+1);
491  box->SetFillColor(kWhite);
492  box->Draw();
493 
494  TText * text_na = new TText(.5,.5,"N/A");
495  text_na->SetTextFont(42);
496  text_na->SetTextAlign(22);
497  text_na->SetTextSize(.8/(double)npar);
498  text_na->SetTextColor(kGray+1);
499  text_na->Draw();
500  }
501  // otherwise draw the histogram
502  else {
503  if(i==j) {
504  bh1->Draw("BTsiB3CS1D0");
505  }
506  else {
507  bh2->Draw("BTfB3CS1nL");
508  }
509 
510  hh->GetXaxis()->SetLabelOffset(5500);
511  hh->GetYaxis()->SetLabelOffset(5500);
512  hh->GetXaxis()->SetTitleSize(10.00);
513  hh->GetYaxis()->SetTitleSize(10.00);
514  }
515 
516  c->cd();
517 
518  // draw axis label
519  double labelsize = .8/(double)npar/10.;
520  double xtext, ytext;
521 
522  // y axis
523  if(i==0) {
524  TLatex * label = new TLatex();
525  label->SetTextFont(62);
526  label->SetTextSize(labelsize);
527  label->SetTextAlign(yalignment);
528  label->SetNDC();
529 
530  label->SetTextAngle(90 - rotation);
531 
532  xtext = margin * (1. - 8. * labelsize);
533  ytext = yup - padsize / 2.;
534 
535  label->DrawLatex(xtext,ytext,fModel->GetParameter(j)->GetLatexName().c_str());
536  }
537 
538  // x axis
539  if(j==npar-1) {
540  TLatex * label = new TLatex();
541  label->SetTextFont(62);
542  label->SetTextSize(labelsize);
543  label->SetTextAlign(xalignment);
544  label->SetNDC();
545 
546  label->SetTextAngle(360 - rotation);
547 
548  xtext = xlow + padsize / 2.;
549  ytext = margin * (1. - 6. * labelsize);
550 
551  label->DrawLatex(xtext,ytext, fModel->GetParameter(i)->GetLatexName().c_str());
552  }
553  }
554  }
555 
556  gPad->RedrawAxis();
557  c->Print(filename);
558 
559  return 1;
560 }
int BCSummaryTool::PrintKnowledgeUpdatePlot1D ( int  index,
const char *  filename,
std::string  options_post = "",
std::string  options_prior = "" 
)

Print a comparison of the prior knowledge to the posterior knowledge for each parameter.

Returns
An error flag.

Definition at line 563 of file BCSummaryTool.cxx.

564 {
565  // perform analysis
566  CalculatePriorModel();
567 
568  // create canvas
569  TCanvas * c = new TCanvas();
570  c->cd();
571 
572  // draw
573  DrawKnowledgeUpdatePlot1D(index, options_post, options_prior);
574 
575  // print
576  c->Print(filename);
577 
578  // no error
579  return 1;
580 }
int BCSummaryTool::PrintKnowledgeUpdatePlots ( const char *  filename = "update.pdf",
std::string  options = "" 
)

Print a comparison of the prior knowledge to the posterior knowledge for each parameter.

Returns
An error flag.

Definition at line 709 of file BCSummaryTool.cxx.

710 {
711  // perform analysis
712  CalculatePriorModel();
713 
714  // option flags
715  bool flag_slice = false;
716 
717  // check content of options string
718  if (options.find("slice") < options.size()) {
719  flag_slice = true;
720  }
721 
722  std::string file(filename);
723 
724  // check if file extension is pdf
725  if ( (file.find_last_of(".") != std::string::npos) &&
726  (file.substr(file.find_last_of(".")+1) == "pdf") ) {
727  ; // it's a PDF file
728 
729  }
730  else if ( (file.find_last_of(".") != std::string::npos) &&
731  (file.substr(file.find_last_of(".")+1) == "ps") ) {
732  ; // it's a PS file
733  }
734  else {
735  ; // make it a PDF file
736  file += ".pdf";
737  }
738 
739  // create canvas and prepare postscript
740  TCanvas * c = new TCanvas(TString::Format("c_%d",getNextIndex()));
741  c->cd();
742 
743  // loop over all parameters and draw 1D plots
744  int npar = fModel->GetNParameters();
745  c->Print(std::string(file + "[").c_str());
746  for (int i = 0; i < npar; ++i) {
747  if ( !DrawKnowledgeUpdatePlot1D(i, options, options))
748  continue;
749  c->Print(file.c_str());
750  }
751 
752  // create legend
753  TLegend * legend2d = new TLegend();
754  legend2d->SetBorderSize(0);
755  legend2d->SetFillColor(0);
756  legend2d->SetTextAlign(12);
757  legend2d->SetTextFont(62);
758  legend2d->SetTextSize(0.03);
759 
760  // create markers and arrows
761  TMarker * marker_prior = new TMarker();
762  marker_prior->SetMarkerStyle(24);
763  marker_prior->SetMarkerColor(kRed);
764 
765  TMarker * marker_posterior = new TMarker();
766  marker_posterior->SetMarkerStyle(24);
767  marker_posterior->SetMarkerColor(kBlack);
768 
769  TArrow * arrow = new TArrow();
770  arrow->SetArrowSize(0.02);
771  arrow->SetLineColor(kBlue);
772  // arrow->SetLineStyle(2);
773 
774  // loop over all parameters
775  for (int i = 0; i < npar; ++i) {
776  for (int j = 0; j < i; ++j) {
777  c->cd();
778 
779  // get parameters
780  const BCParameter * par1 = fModel->GetParameter(i);
781  const BCParameter * par2 = fModel->GetParameter(j);
782 
783  // get 2-d histograms
784  BCH2D* bch2d_2dprior = 0;
785  BCH2D* bch2d_2dposterior = 0;
786  if (flag_slice && npar == 2) {
787  bch2d_2dprior = fPriorModel->GetSlice(par2, par1);
788  bch2d_2dposterior = fModel->GetSlice(par2, par1);
789  }
790  else {
791  bch2d_2dprior = fPriorModel->GetMarginalized(par1, par2);
792  bch2d_2dposterior = fModel->GetMarginalized(par1, par2);
793  }
794 
795  // can't draw anything
796  if ( !bch2d_2dprior || !bch2d_2dposterior)
797  continue;
798 
799  // get histograms
800  TH2D* hist_2dprior = bch2d_2dprior->GetHistogram();
801  hist_2dprior->SetLineColor(kRed);
802  TH2D* hist_2dposterior = bch2d_2dposterior->GetHistogram();
803 
804  // scale histograms
805  hist_2dprior->Scale(1.0/hist_2dprior->Integral("width"));
806  hist_2dposterior->Scale(1.0/hist_2dposterior->Integral("width"));
807 
808  // calculate contours
809  bch2d_2dprior->CalculateIntegratedHistogram();
810  bch2d_2dposterior->CalculateIntegratedHistogram();
811 
812  double level[1] = {bch2d_2dprior->GetLevel(0.32)};
813  hist_2dprior->SetContour(1, level);
814  hist_2dprior->Draw("CONT3");
815  level[0] = bch2d_2dposterior->GetLevel(0.32);
816  hist_2dposterior->SetContour(1, level);
817  hist_2dposterior->Draw("CONT3 SAME");
818 
819  std::vector<double> mode_prior = fPriorModel->GetBestFitParameters();
820  std::vector<double> mode_posterior = fModel->GetBestFitParameters();
821 
822  marker_prior->DrawMarker(mode_prior.at(j), mode_prior.at(i));
823  marker_posterior->DrawMarker(mode_posterior.at(j), mode_posterior.at(i));
824  arrow->DrawArrow(mode_prior.at(j), mode_prior.at(i), mode_posterior.at(j), mode_posterior.at(i));
825 
826  if (i==1 && j == 0) {
827  legend2d->AddEntry(hist_2dprior, "smallest 68% interval(s) of prior", "L");
828  legend2d->AddEntry(hist_2dposterior, "smallest 68% interval(s) of posterior", "L");
829  legend2d->AddEntry(marker_prior, "prior mode", "P");
830  legend2d->AddEntry(marker_posterior, "posterior mode", "P");
831  legend2d->AddEntry(arrow, "change in mode", "L");
832  }
833  gPad->SetTopMargin(0.02);
834 
835  double height = 0.03*legend2d->GetNRows();
836 
837  double xlegend1 = gPad->GetLeftMargin();
838  double xlegend2 = 1.0-gPad->GetRightMargin();
839  double ylegend1 = 1.-gPad->GetTopMargin()-height;
840  double ylegend2 = 1.-gPad->GetTopMargin();
841 
842  // place legend on top of histogram
843  legend2d->SetX1NDC(xlegend1);
844  legend2d->SetX2NDC(xlegend2);
845  legend2d->SetY1NDC(ylegend1);
846  legend2d->SetY2NDC(ylegend2);
847 
848  legend2d->Draw();
849 
850  gPad->SetTopMargin(1.-ylegend1+0.01);
851 
852  gPad->RedrawAxis();
853  c->Print(file.c_str());
854  }
855  }
856 
857  // close output
858  c->Print(std::string(file + "]").c_str());
859  c->Update();
860 
861  // free memory
862  delete legend2d;
863  delete marker_prior;
864  delete marker_posterior;
865  delete arrow;
866  delete c;
867 
868  // no error
869  return 1;
870 }
int BCSummaryTool::PrintParameterLatex ( const char *  filename)

Print a Latex table of the parameters.

Returns
An error flag.

Definition at line 873 of file BCSummaryTool.cxx.

874 {
875  // open file
876  std::ofstream ofi(filename);
877  ofi.precision(3);
878 
879  // check if file is open
880  if(!ofi.is_open()) {
881  std::cerr << "Couldn't open file " << filename <<std::endl;
882  return 0;
883  }
884 
885  // get number of parameters and quantiles
886  int npar = fModel->GetNParameters();
887 
888  // print table
889  ofi
890  << "\\documentclass[11pt, a4paper]{article}" << std::endl
891  << std::endl
892  << "\\begin{document}" << std::endl
893  << std::endl
894  << "\\begin{table}[ht!]" << std::endl
895  << "\\begin{center}" << std::endl
896  <<"\\begin{tabular}{llllllll}" << std::endl
897  << "\\hline" << std::endl
898  << "Parameter & Mean & RMS & Gl. mode & Mode & Median & 16\\% quant. & 84\\% quant. \\\\" << std::endl
899  << "\\hline" << std::endl;
900 
901  for (int i = 0; i < npar; ++i) {
902  const BCParameter * par = fModel->GetParameter(i);
903  BCH1D * bch1d = fModel->GetMarginalized(par);
904  ofi
905  << par->GetName() << " & "
906  << bch1d->GetMean() << " & "
907  << bch1d->GetRMS() << " & "
908  << fModel->GetBestFitParameters().at(i) << " & "
909  << bch1d->GetMode() << " & "
910  << bch1d->GetMedian() << " & "
911  << bch1d->GetQuantile(0.16) << " & "
912  << bch1d->GetQuantile(0.84) << " \\\\" << std::endl;
913  }
914  ofi
915  << "\\hline" << std::endl
916  << "\\end{tabular}" << std::endl
917  << "\\caption{Summary of the parameter estimates.}" << std::endl
918  << "\\end{center}" << std::endl
919  << "\\end{table}" << std::endl
920  << std::endl
921  << "\\end{document}" << std::endl;
922 
923  // close file
924  ofi.close();
925 
926  // no error
927  return 1;
928 }
int BCSummaryTool::PrintParameterPlot ( const char *  filename = "parameters.pdf")

Print a summary plot for the parameters.

Returns
An error flag.

Definition at line 166 of file BCSummaryTool.cxx.

167 {
168  // copy summary data
169  if (!CopySummaryData())
170  return 0;
171 
172  // get number of parameters and quantiles
173  int npar = fModel->GetNParameters();
174  int nquantiles = int( fSumProb.size() );
175 
176  // create histogram
177  TH1D * hist_axes = new TH1D(
178  TString::Format("hist_axes_par_%d",getNextIndex()),
179  ";;Scaled parameter range [a.u.]",npar, -0.5, npar-0.5);
180  hist_axes->SetStats(kFALSE);
181  for (int i = 0; i < npar; ++i)
182  hist_axes->GetXaxis()->SetBinLabel( i+1, fParName.at(i).c_str() );
183  hist_axes->GetXaxis()->SetLabelOffset(0.015);
184  hist_axes->GetXaxis()->SetLabelSize(0.06);
185  hist_axes->GetXaxis()->SetTickLength(0.0);
186  hist_axes->GetYaxis()->SetRangeUser(-0.1, 1.1);
187  hist_axes->GetYaxis()->SetTickLength(0.0);
188 
189  // create graphs
190  TGraphErrors * graph_quantiles = new TGraphErrors(npar*nquantiles);
191  graph_quantiles->SetMarkerSize(0);
192  graph_quantiles->SetLineColor(38);
193  graph_quantiles->SetLineStyle(2);
194 
195  TGraphErrors * graph_mean = new TGraphErrors(npar);
196  graph_mean->SetMarkerColor(kBlack);
197  graph_mean->SetMarkerStyle(20);
198 
199  TGraphErrors * graph_mode = new TGraphErrors(npar);
200  graph_mode->SetMarkerColor(kRed);
201  graph_mode->SetMarkerStyle(20);
202 
203  TGraphAsymmErrors * graph_intervals = new TGraphAsymmErrors(0);
204  graph_intervals->SetFillColor(kYellow);
205  graph_intervals->SetLineStyle(2);
206  graph_intervals->SetLineColor(kRed);
207  graph_intervals->SetMarkerSize(0);
208 
209  // fill graphs
210  int indexintervals = 0;
211 
212  // fill graph quantiles
213  if (fFlagInfoMarg) {
214  for (int i = 0; i < npar; ++i) {
215  for (int j = 0; j < nquantiles; ++j) {
216  graph_quantiles->SetPoint(i*nquantiles+j,double(i),
217  (fQuantiles.at(i*nquantiles+j) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
218  graph_quantiles->SetPointError(i*nquantiles+j, 0.5, 0.0);
219  }
220  }
221  }
222 
223  // fill graph mean and rms
224  if (fFlagInfoMarg) {
225  for (int i = 0; i < npar; ++i) {
226  // fill graph mean
227  graph_mean->SetPoint(i, double(i), (fMean.at(i) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
228  graph_mean->SetPointError(i, 0.0, fRMS.at(i)/(fParMax.at(i)-fParMin.at(i)));
229  }
230  }
231 
232  // fill graph mode
233  if (fFlagInfoOpt) {
234  for (int i = 0; i < npar; ++i)
235  graph_mode->SetPoint(i, double(i), (fGlobalMode.at(i) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
236  }
237 
238  // fill graph smallest intervals
239  if (fFlagInfoMarg) {
240  for (int i = 0; i < npar; ++i) {
241  int nintervals = int(fSmallInt.at(indexintervals++));
242  for (int j = 0; j < nintervals; ++j) {
243  double xmin = fSmallInt.at(indexintervals++);
244  double xmax = fSmallInt.at(indexintervals++);
245  indexintervals++;
246  double xlocalmaxpos = fSmallInt.at(indexintervals++);
247  indexintervals++;
248  int npoints = graph_intervals->GetN();
249  graph_intervals->SetPoint(npoints,double(i),
250  (xlocalmaxpos - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
251  graph_intervals->SetPointError(npoints,0.5, 0.5,
252  (xlocalmaxpos - xmin)/(fParMax.at(i)-fParMin.at(i)),
253  (xmax - xlocalmaxpos)/(fParMax.at(i)-fParMin.at(i)));
254  }
255  }
256  }
257 
258  // create legend
259  TLegend * legend = new TLegend(0.15, 0.88, 0.85, 0.99);
260  legend->SetBorderSize(0);
261  legend->SetFillColor(0);
262 
263  // create latex
264  TLatex * latex = new TLatex();
265  latex->SetTextSize(0.02);
266 
267  // create lines
268  TLine * line_top = new TLine(-0.5, 1.0, npar-0.5, 1.0);
269  line_top->SetLineColor(kBlack);
270  line_top->SetLineStyle(1);
271  line_top->SetLineWidth(2);
272 
273  TLine * line_bot = new TLine(-0.5, 0.0, npar-0.5, 0.0);
274  line_bot->SetLineColor(kBlack);
275  line_bot->SetLineStyle(1);
276  line_bot->SetLineWidth(2);
277 
278  // print to file
279  TCanvas * c_par = new TCanvas(TString::Format("c_par_%d",getNextIndex()));
280  c_par->cd();
281  hist_axes->Draw();
282  line_top->Draw();
283  line_bot->Draw();
284  if (fFlagInfoMarg) {
285  graph_intervals->DrawClone("SAME2");
286  for (int i = 0; i < graph_intervals->GetN(); ++i)
287  graph_intervals->SetPointError(i, 0.5, 0.5, 0.0, 0.0);
288  graph_intervals->Draw("SAMEPZ");
289  graph_quantiles->Draw("SAMEPZ");
290  graph_mean->Draw("SAMEP");
291  legend->AddEntry(graph_quantiles, "Quantiles (5%, 10%, 16%, 50%, 84%, 90%, 95%)", "L");
292  legend->AddEntry(graph_mean, "Mean and RMS", "LEP");
293  legend->AddEntry(graph_intervals, "Smallest 68% intervals and local modes", "FL");
294  }
295  if (fFlagInfoOpt) {
296  graph_mode->Draw("SAMEP");
297  legend->AddEntry(graph_mode, "Global mode", "P");
298  }
299  for (int i = 0; i < npar;++i) {
300  // latex->DrawLatex(double(i)-0.1, 0.010, Form("%+3.3f", fParMin.at(i)));
301  // latex->DrawLatex(double(i)-0.1, 0.965, Form("%+3.3f", fParMax.at(i)));
302  latex->DrawLatex(double(i)-0.1, 0.010-0.07, Form("%+3.3g", fParMin.at(i)));
303  latex->DrawLatex(double(i)-0.1, 0.965+0.07, Form("%+3.3g", fParMax.at(i)));
304  }
305  latex->SetNDC();
306  latex->DrawLatex(0.9, 0.175, "Par. min.");
307  latex->DrawLatex(0.9, 0.83, "Par. max.");
308  legend->Draw("SAME");
309  gPad->RedrawAxis();
310  c_par->Print(filename);
311 
312  // no error
313  return 1;
314 }
void BCSummaryTool::SetModel ( BCModel *  model)
inline

Set the model to be summarized.

Parameters
modelThe BCModel to be summarized.

Definition at line 69 of file BCSummaryTool.h.

70  { fModel = model; };

Member Data Documentation

std::vector<double> BCSummaryTool::fCorrCoeff
private

Correlation coefficients. Length of vector equals number of parameters * number of parameters.

Definition at line 156 of file BCSummaryTool.h.

bool BCSummaryTool::fFlagInfoMarg
private

A flag: check if marginalized information is present

Definition at line 202 of file BCSummaryTool.h.

bool BCSummaryTool::fFlagInfoOpt
private

A flag: check if optimization information is present

Definition at line 206 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fGlobalMode
private

Global modes.
Length of vector equals number of parameters.

Definition at line 171 of file BCSummaryTool.h.

unsigned int BCSummaryTool::fHCounter =0
staticprivate

helper variable to get an unique number to be used in histogram naming

Definition at line 135 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fMargMode
private

Marginalized modes.
Length of vector equals number of parameters.

Definition at line 161 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fMean
private

Mean values.
Length of vector equals number of parameters.

Definition at line 166 of file BCSummaryTool.h.

BCModel* BCSummaryTool::fModel
private

The model whose results are summarized

Definition at line 139 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fParMax
private

Parameter maxima

Definition at line 151 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fParMin
private

parameter minima

Definition at line 147 of file BCSummaryTool.h.

std::vector<std::string> BCSummaryTool::fParName
private

parameter names

Definition at line 143 of file BCSummaryTool.h.

BCSummaryPriorModel* BCSummaryTool::fPriorModel
private

A model for calculating the marginalized distributions for the prior probabilities.

Definition at line 198 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fQuantiles
private

Quantiles.
The following quantiles are stored: 0.05, 0.10, 0.16, 0.5, 0.84, 0.90, 0.95.
Length of vector equals number of parameters * number of quantiles.

Definition at line 177 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fRMS
private

RMS values.
Length of vector equals number of parameters.

Definition at line 189 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fSmallInt
private

Smallest intervals.
For each parameter a set of the smallest intervals is recorded.
Structure: number of intervals n + n * (start, stop, local max, local max pos, integral) Length of vector equals number of parameters * number of quantiles.

Definition at line 184 of file BCSummaryTool.h.

std::vector<double> BCSummaryTool::fSumProb
private

Sum of probabilities for quantiles

Definition at line 193 of file BCSummaryTool.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.4