MaCh3  2.2.3
Reference Guide
Functions | Variables
RHat.cpp File Reference

This executable calculates the \( \hat{R} \) estimator for Markov Chain Monte Carlo (MCMC) convergence. More...

#include "Manager/Manager.h"
#include "Samples/SampleStructs.h"
#include "TObjArray.h"
#include "TChain.h"
#include "TFile.h"
#include "TBranch.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TLegend.h"
#include "TString.h"
#include "TH1.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TColor.h"
#include "TStyle.h"
#include "TROOT.h"
Include dependency graph for RHat.cpp:

Go to the source code of this file.

Functions

void PrepareChains ()
 
void InitialiseArrays ()
 
void RunDiagnostic ()
 
void CalcRhat ()
 
void SaveResults ()
 
void DestroyArrays ()
 
double CalcMedian (double arr[], int size)
 
void CapVariable (double var, double cap)
 
int main (int argc, char *argv[])
 

Variables

int * Ntoys_requested
 
int * Ntoys_filled
 
int TotToys
 
int NThin
 
int Nchains
 
int nDraw
 
std::vector< TString > BranchNames
 
std::vector< std::string > MCMCFile
 
std::vector< bool > ValidPar
 
double * S1_global
 
double * S2_global
 
double ** S1_chain
 
double ** S2_chain
 
double ** Mean
 
double ** StandardDeviation
 
double * MeanGlobal
 
double * StandardDeviationGlobal
 
double * BetweenChainVariance
 
double * MarginalPosteriorVariance
 
double * RHat
 
double * EffectiveSampleSize
 

Detailed Description

This executable calculates the \( \hat{R} \) estimator for Markov Chain Monte Carlo (MCMC) convergence.

KS: This exe is meant to calculate the \( \hat{R} \) estimator. For a well-converged chain, this distribution should be centered at one. The \( \hat{R} \) statistic is used to assess the convergence of MCMC simulations and helps determine whether the chains have reached a stable distribution.

[32].

MJR: Update – Improved memory usage so that whole chains can be quickly loaded without requiring copious amounts of RAM. This comes at the cost of not being able to calculate the Folded RHat since finding the median requires the loading of full chains at a time. The method has been validated to give identical results to the "High Memory" (original) version at a fraction of the runtime and resources.

The input format is also slightly altered; since we can now load entire chains, there's less need to specify how many toys are desired for a sub-sample, so the Ntoys input has been removed.

KS: This exe is meant to calculate the \( \hat{R} \) estimator. For a well-converged chain, this distribution should be centered at one. The \( \hat{R} \) statistic is used to assess the convergence of MCMC simulations and helps determine whether the chains have reached a stable distribution.

[32].

Definition in file RHat.cpp.

Function Documentation

◆ CalcMedian()

double CalcMedian ( double  arr[],
int  size 
)

Definition at line 680 of file RHat.cpp.

680  {
681 // *******************
682  std::sort(arr, arr+size);
683  if (size % 2 != 0)
684  return arr[size/2];
685  return (arr[(size-1)/2] + arr[size/2])/2.0;
686 }
int size

◆ CalcRhat()

void CalcRhat ( )

Definition at line 383 of file RHat.cpp.

383  {
384 // *******************
385 
386  TStopwatch clock;
387  clock.Start();
388 
389  //KS: Start parallel region
390  // If we would like to do this for thousands of chains we might consider using GPU for this
391  #ifdef MULTITHREAD
392  #pragma omp parallel
393  {
394  #endif
395 
396  #ifdef MULTITHREAD
397  #pragma omp for
398  #endif
399  //KS: loop over chains and draws are independent so might as well collapse for sweet cache hits
400  //Calculate the mean for each parameter within each considered chain
401  // MJR: Calculate using running totals to massively save on time and memory
402  for (int m = 0; m < Nchains; ++m)
403  {
404  for (int j = 0; j < nDraw; ++j)
405  {
406  Mean[m][j] = S1_chain[m][j] / static_cast<double>(Ntoys_filled[m]);
407  StandardDeviation[m][j] = S2_chain[m][j]/static_cast<double>(Ntoys_filled[m]) - Mean[m][j]*Mean[m][j];
408  }
409  }
410 
411  #ifdef MULTITHREAD
412  #pragma omp for
413  #endif
414  //Calculate the mean for each parameter global means we include information from several chains
415  for (int j = 0; j < nDraw; ++j)
416  {
417  for (int m = 0; m < Nchains; ++m)
418  {
420  }
421  MeanGlobal[j] = S1_global[j] / static_cast<double>(TotToys);
422  StandardDeviationGlobal[j] = StandardDeviationGlobal[j] / static_cast<double>(Nchains);
423  }
424 
425  #ifdef MULTITHREAD
426  #pragma omp for
427  #endif
428  for (int j = 0; j < nDraw; ++j)
429  {
430  //KS: This term only makes sense if we have at least 2 chains
431  if(Nchains == 1)
432  {
433  BetweenChainVariance[j] = 0.;
434  }
435  else
436  {
437  for (int m = 0; m < Nchains; ++m)
438  {
439  BetweenChainVariance[j] += ( Mean[m][j] - MeanGlobal[j])*( Mean[m][j] - MeanGlobal[j]) * Ntoys_filled[m];
440  }
442  }
443  }
444 
445  int avgNtoys = TotToys/Nchains;
446  #ifdef MULTITHREAD
447  #pragma omp for
448  #endif
449  for (int j = 0; j < nDraw; ++j)
450  {
451  MarginalPosteriorVariance[j] = (avgNtoys-1) * StandardDeviationGlobal[j] / (avgNtoys) + BetweenChainVariance[j]/avgNtoys;
452  }
453 
454  #ifdef MULTITHREAD
455  #pragma omp for
456  #endif
457  //Finally calculate our estimator
458  for (int j = 0; j < nDraw; ++j)
459  {
461 
462  //KS: For flat params values can be crazy so cap at 0
463  CapVariable(RHat[j], 0);
464  }
465 
466  #ifdef MULTITHREAD
467  #pragma omp for
468  #endif
469  //KS: Additionally calculates effective step size which is an estimate of the sample size required to achieve the same level of precision if that sample was a simple random sample.
470  for (int j = 0; j < nDraw; ++j)
471  {
473 
474  //KS: For flat params values can be crazy so cap at 0
476  }
477  #ifdef MULTITHREAD
478  } //End parallel region
479  #endif
480 
481  clock.Stop();
482  MACH3LOG_INFO("Finished calculating RHat, it took {:.2f}s to finish", clock.RealTime());
483 }
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:25
double * S1_global
Definition: RHat.cpp:53
double ** S1_chain
Definition: RHat.cpp:55
int Nchains
Definition: RHat.cpp:46
void CapVariable(double var, double cap)
Definition: RHat.cpp:690
int * Ntoys_filled
Definition: RHat.cpp:43
double * MeanGlobal
Definition: RHat.cpp:61
double * BetweenChainVariance
Definition: RHat.cpp:64
double ** StandardDeviation
Definition: RHat.cpp:59
double * StandardDeviationGlobal
Definition: RHat.cpp:62
double * MarginalPosteriorVariance
Definition: RHat.cpp:65
int TotToys
Definition: RHat.cpp:44
double * RHat
Definition: RHat.cpp:66
double ** Mean
Definition: RHat.cpp:58
double ** S2_chain
Definition: RHat.cpp:56
int nDraw
Definition: RHat.cpp:47
double * EffectiveSampleSize
Definition: RHat.cpp:67

◆ CapVariable()

void CapVariable ( double  var,
double  cap 
)

Definition at line 690 of file RHat.cpp.

690  {
691 // *******************
692  if(std::isnan(var) || !std::isfinite(var)) var = cap;
693 }

◆ DestroyArrays()

void DestroyArrays ( )

Definition at line 649 of file RHat.cpp.

649  {
650 // *******************
651 
652  MACH3LOG_INFO("Killing all arrays");
653  delete[] MeanGlobal;
654  delete[] StandardDeviationGlobal;
655  delete[] BetweenChainVariance;
656  delete[] MarginalPosteriorVariance;
657  delete[] RHat;
658  delete[] EffectiveSampleSize;
659 
660  for(int m = 0; m < Nchains; m++)
661  {
662  delete[] Mean[m];
663  delete[] StandardDeviation[m];
664  delete[] S1_chain[m];
665  delete[] S2_chain[m];
666  }
667  delete[] Mean;
668  delete[] StandardDeviation;
669  delete[] S1_chain;
670  delete[] S2_chain;
671  delete[] S1_global;
672  delete[] S2_global;
673 
674  delete[] Ntoys_requested;
675  delete[] Ntoys_filled;
676 }
double * S2_global
Definition: RHat.cpp:54
int * Ntoys_requested
Definition: RHat.cpp:42

◆ InitialiseArrays()

void InitialiseArrays ( )

Definition at line 335 of file RHat.cpp.

335  {
336 // *******************
337 
338  MACH3LOG_INFO("Initialising arrays");
339  Mean = new double*[Nchains]();
340  StandardDeviation = new double*[Nchains]();
341 
342  MeanGlobal = new double[nDraw]();
343  StandardDeviationGlobal = new double[nDraw]();
344  BetweenChainVariance = new double[nDraw]();
345 
346  MarginalPosteriorVariance = new double[nDraw]();
347  RHat = new double[nDraw]();
348  EffectiveSampleSize = new double[nDraw]();
349 
350  for (int m = 0; m < Nchains; ++m)
351  {
352  Mean[m] = new double[nDraw]();
353  StandardDeviation[m] = new double[nDraw]();
354 
355  for (int j = 0; j < nDraw; ++j)
356  {
357  Mean[m][j] = 0.;
358  StandardDeviation[m][j] = 0.;
359 
360  if(m == 0)
361  {
362  MeanGlobal[j] = 0.;
363  StandardDeviationGlobal[j] = 0.;
364  BetweenChainVariance[j] = 0.;
366  RHat[j] = 0.;
367  EffectiveSampleSize[j] = 0.;
368  }
369  }
370  }
371 }

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 83 of file RHat.cpp.

83  {
84 // *******************
85 
88 
89  Mean = nullptr;
90  StandardDeviation = nullptr;
91 
92  MeanGlobal = nullptr;
93  StandardDeviationGlobal = nullptr;
94 
95  BetweenChainVariance = nullptr;
96  MarginalPosteriorVariance = nullptr;
97  RHat = nullptr;
98  EffectiveSampleSize = nullptr;
99 
100  Nchains = 0;
101 
102  if (argc < 2)
103  {
104  MACH3LOG_ERROR("Wrong arguments");
105  MACH3LOG_ERROR("./RHat NThin MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]");
106  throw MaCh3Exception(__FILE__ , __LINE__ );
107  }
108 
109  NThin = atoi(argv[1]);
110 
111  //KS Gelman suggests to diagnose on more than one chain
112  for (int i = 2; i < argc; i++)
113  {
114  MCMCFile.push_back(std::string(argv[i]));
115  MACH3LOG_INFO("Adding file: {}", MCMCFile.back());
116  Nchains++;
117  }
118 
119  if(Nchains == 1)
120  {
121  MACH3LOG_WARN("Gelman is going to be sad :(. He suggested you should use more than one chain (at least 4). Code works fine for one chain, however, estimator might be biased.");
122  MACH3LOG_WARN("Multiple chains are more likely to reveal multimodality and poor adaptation or mixing:");
123  }
124  MACH3LOG_INFO("Diagnosing {} chains", Nchains);
125 
126  PrepareChains();
127 
129 
130  //KS: Main function
131  RunDiagnostic();
132 
133  SaveResults();
134 
135  DestroyArrays();
136 
137  return 0;
138 }
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:27
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:51
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:26
void SaveResults()
Definition: RHat.cpp:487
void InitialiseArrays()
Definition: RHat.cpp:335
void RunDiagnostic()
Definition: RHat.cpp:374
int NThin
Definition: RHat.cpp:45
std::vector< std::string > MCMCFile
Definition: RHat.cpp:50
void DestroyArrays()
Definition: RHat.cpp:649
void PrepareChains()
Definition: RHat.cpp:142
Custom exception class for MaCh3 errors.
void MaCh3Welcome()
KS: Prints welcome message with MaCh3 logo.
Definition: Monitor.cpp:12

◆ PrepareChains()

void PrepareChains ( )

Definition at line 142 of file RHat.cpp.

142  {
143 // *******************
144 
145  TStopwatch clock;
146  clock.Start();
147 
148  Ntoys_requested = new int[Nchains]();
149  Ntoys_filled = new int[Nchains]();
150  TotToys = 0;
151  std::vector<int> BurnIn(Nchains);
152  std::vector<int> nEntries(Nchains);
153  std::vector<int> nBranches(Nchains);
154  std::vector<int> step(Nchains);
155 
156  S1_chain = new double*[Nchains]();
157  S2_chain = new double*[Nchains]();
158 
159  // KS: This can reduce time necessary for caching even by half
160  #ifdef MULTITHREAD
161  //ROOT::EnableImplicitMT();
162  #endif
163 
164  // Open the Chain
165  //It is tempting to multithread here but unfortunately, ROOT files are not thread safe :(
166  for (int m = 0; m < Nchains; m++)
167  {
168  TChain* Chain = new TChain("posteriors");
169  Chain->Add(MCMCFile[m].c_str());
170 
171  nEntries[m] = int(Chain->GetEntries());
172  Ntoys_requested[m] = nEntries[m]/NThin;
173  Ntoys_filled[m] = 0;
174 
175  MACH3LOG_INFO("On file: {}", MCMCFile[m].c_str());
176  MACH3LOG_INFO("Generating {} Toys", Ntoys_requested[m]);
177 
178  // Set the step cut to be 20%
179  BurnIn[m] = nEntries[m]/5;
180 
181  // Get the list of branches
182  TObjArray* brlis = Chain->GetListOfBranches();
183 
184  // Get the number of branches
185  nBranches[m] = brlis->GetEntries();
186 
187  if(m == 0) BranchNames.reserve(nBranches[m]);
188 
189  // Set all the branches to off
190  Chain->SetBranchStatus("*", false);
191 
192  // Loop over the number of branches
193  // Find the name and how many of each systematic we have
194  for (int i = 0; i < nBranches[m]; i++)
195  {
196  // Get the TBranch and its name
197  TBranch* br = static_cast<TBranch *>(brlis->At(i));
198  if(!br){
199  MACH3LOG_ERROR("Invalid branch at position {}", i);
200  throw MaCh3Exception(__FILE__,__LINE__);
201  }
202  TString bname = br->GetName();
203 
204  // Read in the step
205  if (bname == "step") {
206  Chain->SetBranchStatus(bname, true);
207  Chain->SetBranchAddress(bname, &step[m]);
208  }
209  //Count all branches
210  else if (bname.BeginsWith("PCA_") || bname.BeginsWith("accProb") || bname.BeginsWith("stepTime") )
211  {
212  continue;
213  }
214  else
215  {
216  //KS: Save branch name only for one chain, we assume all chains have the same branches, otherwise this doesn't make sense either way
217  if(m == 0)
218  {
219  BranchNames.push_back(bname);
220  //KS: We calculate R Hat also for LogL, just in case, however we plot them separately
221  if(bname.BeginsWith("LogL"))
222  {
223  ValidPar.push_back(false);
224  }
225  else
226  {
227  ValidPar.push_back(true);
228  }
229  }
230  Chain->SetBranchStatus(bname, true);
231  MACH3LOG_DEBUG("{}", bname);
232  }
233  }
234 
235  if(m == 0) nDraw = int(BranchNames.size());
236 
237  // MJR: Initialize quantities needed for calculating RHat
238  S1_chain[m] = new double[nDraw]();
239  S2_chain[m] = new double[nDraw]();
240  if (m == 0)
241  {
242  S1_global = new double[nDraw]();
243  S2_global = new double[nDraw]();
244  }
245  for (int id = 0; id < nDraw; ++id)
246  {
247  S1_chain[m][id] = 0.0;
248  S2_chain[m][id] = 0.0;
249  if (m == 0)
250  {
251  S1_global[id] = 0.0;
252  S2_global[id] = 0.0;
253  }
254  }
255 
256  //TN: Qualitatively faster sanity check, with the very same outcome (all chains have the same #branches)
257  if(m > 0)
258  {
259  if(nBranches[m] != nBranches[0])
260  {
261  MACH3LOG_ERROR("Ups, something went wrong, chain {} called {} has {} branches, while 0 called {} has {} branches", m, MCMCFile[m], nBranches[m], MCMCFile[0], nBranches[0]);
262  MACH3LOG_ERROR("All chains should have the same number of branches");
263  throw MaCh3Exception(__FILE__ , __LINE__ );
264  }
265  }
266 
267  // MJR: Create an array to hold branch values. Resetting branch addresses
268  // for every step is very expensive.
269  double* branch_values = new double[nDraw]();
270  for (int id = 0; id < nDraw; ++id)
271  {
272  Chain->SetBranchAddress(BranchNames[id].Data(), &branch_values[id]);
273  }
274 
275  //TN: move looping over toys here, so we don't need to loop over chains more than once
276  if(BurnIn[m] >= nEntries[m])
277  {
278  MACH3LOG_ERROR("You are running on a chain shorter than BurnIn cut");
279  MACH3LOG_ERROR("Number of entries {} BurnIn cut {}", nEntries[m], BurnIn[m]);
280  MACH3LOG_ERROR("You will run into the infinite loop");
281  MACH3LOG_ERROR("You can make a new chain or modify BurnIn cut");
282  throw MaCh3Exception(__FILE__ , __LINE__ );
283  }
284 
285  MACH3LOG_INFO("Loading chain {} / {}...", m, Nchains);
286  for (int i = 0; i < Ntoys_requested[m]; i++)
287  {
288  // This is here as a placeholder in case we want to do some thinning later
289  int entry = i*NThin;
290 
291  Chain->GetEntry(entry);
292 
293  // If we have combined chains by hadd need to check the step in the chain
294  // Note, entry is not necessarily the same as the step due to merged ROOT files, so can't choose an entry in the range BurnIn - nEntries :(
295  if (step[m] < BurnIn[m])
296  {
297  continue;
298  }
299 
300  // Output some info for the user
301  if (Ntoys_requested[m] > 10 && i % (Ntoys_requested[m]/10) == 0) {
302  MaCh3Utils::PrintProgressBar(i+m*Ntoys_requested[m], static_cast<Long64_t>(Ntoys_requested[m])*Nchains);
303  MACH3LOG_DEBUG("Getting random entry {}", entry);
304  }
305 
306  // MJR: Fill running quantities instead of loading everything into RAM.
307  // This is where we save big on both memory and time (resetting
308  // branch addresses and calling GetEntry() again here is very slow).
309  for (int j = 0; j < nDraw; ++j)
310  {
311  S1_global[j] += branch_values[j];
312  S2_global[j] += branch_values[j]*branch_values[j];
313  S1_chain[m][j] += branch_values[j];
314  S2_chain[m][j] += branch_values[j]*branch_values[j];
315  }
316 
317  // Increment counters
318  Ntoys_filled[m]++;
319  TotToys++;
320 
321  }//end loop over toys
322 
323  //TN: There, we now don't need to keep the chain in memory anymore
324  delete Chain;
325  delete[] branch_values;
326  MACH3LOG_INFO("Finished loading chain {}!", m);
327  }
328 
329  clock.Stop();
330  MACH3LOG_INFO("Finished calculating Toys, it took {:.2f}s to finish", clock.RealTime());
331 }
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:24
std::vector< bool > ValidPar
Definition: RHat.cpp:51
std::vector< TString > BranchNames
Definition: RHat.cpp:49
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:213

◆ RunDiagnostic()

void RunDiagnostic ( )

Definition at line 374 of file RHat.cpp.

374  {
375 // *******************
376  CalcRhat();
377  //In case in future we expand this
378 }
void CalcRhat()
Definition: RHat.cpp:383

◆ SaveResults()

void SaveResults ( )

Definition at line 487 of file RHat.cpp.

487  {
488 // *******************
489  #pragma GCC diagnostic ignored "-Wfloat-conversion"
490 
491  std::string NameTemp = "";
492  //KS: If we run over many many chains there is danger that name will be so absurdly long we run over system limit and job will be killed :(
493  if(Nchains < 5)
494  {
495  for (int i = 0; i < Nchains; i++)
496  {
497  std::string temp = MCMCFile[i];
498 
499  while (temp.find(".root") != std::string::npos) {
500  temp = temp.substr(0, temp.find(".root"));
501  }
502 
503  NameTemp = NameTemp + temp + "_";
504  }
505  }
506  else {
507  NameTemp = std::to_string(Nchains) + "Chains" + "_";
508  }
509  NameTemp += "diag.root";
510 
511  TFile* DiagFile = new TFile(NameTemp.c_str(), "recreate");
512 
513  DiagFile->cd();
514 
515  TH1D *StandardDeviationGlobalPlot = new TH1D("StandardDeviationGlobalPlot", "StandardDeviationGlobalPlot", nDraw, 0, nDraw);
516  TH1D *BetweenChainVariancePlot = new TH1D("BetweenChainVariancePlot", "BetweenChainVariancePlot", nDraw, 0, nDraw);
517  TH1D *MarginalPosteriorVariancePlot = new TH1D("MarginalPosteriorVariancePlot", "MarginalPosteriorVariancePlot", nDraw, 0, nDraw);
518  TH1D *RhatPlot = new TH1D("RhatPlot", "RhatPlot", 200, 0, 2);
519  TH1D *EffectiveSampleSizePlot = new TH1D("EffectiveSampleSizePlot", "EffectiveSampleSizePlot", 400, 0, 10000);
520 
521  TH1D *RhatLogPlot = new TH1D("RhatLogPlot", "RhatLogPlot", 200, 0, 2);
522 
523  int Criterium = 0;
524  for(int j = 0; j < nDraw; j++)
525  {
526  //KS: Fill only valid parameters
527  if(ValidPar[j])
528  {
529  StandardDeviationGlobalPlot->Fill(j,StandardDeviationGlobal[j]);
530  BetweenChainVariancePlot->Fill(j,BetweenChainVariance[j]);
531  MarginalPosteriorVariancePlot->Fill(j,MarginalPosteriorVariance[j]);
532  RhatPlot->Fill(RHat[j]);
533  EffectiveSampleSizePlot->Fill(EffectiveSampleSize[j]);
534  if(RHat[j] > 1.1) Criterium++;
535  }
536  else
537  {
538  RhatLogPlot->Fill(RHat[j]);
539  }
540  }
541  //KS: We set criterium of 1.1 based on Gelman et al. (2003) Bayesian Data Analysis
542  MACH3LOG_WARN("Number of parameters which has R hat greater than 1.1 is {}({:.2f}%)", Criterium, 100*double(Criterium)/double(nDraw));
543  for(int j = 0; j < nDraw; j++)
544  {
545  if( (RHat[j] > 1.1) && ValidPar[j])
546  {
547  MACH3LOG_CRITICAL("Parameter {} has R hat higher than 1.1", BranchNames[j]);
548  }
549  }
550  StandardDeviationGlobalPlot->Write();
551  BetweenChainVariancePlot->Write();
552  MarginalPosteriorVariancePlot->Write();
553  RhatPlot->Write();
554  EffectiveSampleSizePlot->Write();
555 
556  RhatLogPlot->Write();
557 
558  //KS: Now we make fancy canvases, consider some function to have less copy pasting
559  auto TempCanvas = std::make_unique<TCanvas>("Canvas", "Canvas", 1024, 1024);
560  gStyle->SetOptStat(0);
561  TempCanvas->SetGridx();
562  TempCanvas->SetGridy();
563 
564  // Random line to write useful information to TLegend
565  auto TempLine = std::make_unique<TLine>(0, 0, 0, 0);
566  TempLine->SetLineColor(kBlack);
567 
568  RhatPlot->GetXaxis()->SetTitle("R hat");
569  RhatPlot->SetLineColor(kRed);
570  RhatPlot->SetFillColor(kRed);
571 
572  TLegend *Legend = new TLegend(0.55, 0.6, 0.9, 0.9);
573  Legend->SetTextSize(0.04);
574  Legend->SetFillColor(0);
575  Legend->SetFillStyle(0);
576  Legend->SetLineWidth(0);
577  Legend->SetLineColor(0);
578 
579  Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), "");
580  Legend->AddEntry(RhatPlot, "Rhat Gelman 2013", "l");
581 
582  RhatPlot->Draw();
583  Legend->Draw("same");
584  TempCanvas->Write("Rhat");
585  delete Legend;
586  Legend = nullptr;
587 
588  //Now R hat for log L
589  RhatLogPlot->GetXaxis()->SetTitle("R hat for LogL");
590  RhatLogPlot->SetLineColor(kRed);
591  RhatLogPlot->SetFillColor(kRed);
592 
593  Legend = new TLegend(0.55, 0.6, 0.9, 0.9);
594  Legend->SetTextSize(0.04);
595  Legend->SetFillColor(0);
596  Legend->SetFillStyle(0);
597  Legend->SetLineWidth(0);
598  Legend->SetLineColor(0);
599 
600  Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), "");
601  Legend->AddEntry(RhatLogPlot, "Rhat Gelman 2013", "l");
602 
603  RhatLogPlot->Draw();
604  Legend->Draw("same");
605  TempCanvas->Write("RhatLog");
606  delete Legend;
607  Legend = nullptr;
608 
609  //Now canvas for effective sample size
610  EffectiveSampleSizePlot->GetXaxis()->SetTitle("S_{eff, BDA2}");
611  EffectiveSampleSizePlot->SetLineColor(kRed);
612 
613  Legend = new TLegend(0.45, 0.6, 0.9, 0.9);
614  Legend->SetTextSize(0.03);
615  Legend->SetFillColor(0);
616  Legend->SetFillStyle(0);
617  Legend->SetLineWidth(0);
618  Legend->SetLineColor(0);
619 
620  const double Mean1 = EffectiveSampleSizePlot->GetMean();
621  const double RMS1 = EffectiveSampleSizePlot->GetRMS();
622 
623  Legend->AddEntry(TempLine.get(), Form("Number of throws=%.0i, Number of chains=%.1i", TotToys, Nchains), "");
624  Legend->AddEntry(EffectiveSampleSizePlot, Form("S_{eff, BDA2} #mu = %.2f, #sigma = %.2f",Mean1 ,RMS1), "l");
625 
626  EffectiveSampleSizePlot->Draw();
627  Legend->Draw("same");
628  TempCanvas->Write("EffectiveSampleSize");
629 
630  //Fancy memory cleaning
631  delete StandardDeviationGlobalPlot;
632  delete BetweenChainVariancePlot;
633  delete MarginalPosteriorVariancePlot;
634  delete RhatPlot;
635  delete EffectiveSampleSizePlot;
636 
637  delete Legend;
638 
639  delete RhatLogPlot;
640 
641  DiagFile->Close();
642  delete DiagFile;
643 
644  MACH3LOG_INFO("Finished and wrote results to {}", NameTemp);
645 }
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:28

Variable Documentation

◆ BetweenChainVariance

double* BetweenChainVariance

Definition at line 64 of file RHat.cpp.

◆ BranchNames

std::vector<TString> BranchNames

Definition at line 49 of file RHat.cpp.

◆ EffectiveSampleSize

double* EffectiveSampleSize

Definition at line 67 of file RHat.cpp.

◆ MarginalPosteriorVariance

double* MarginalPosteriorVariance

Definition at line 65 of file RHat.cpp.

◆ MCMCFile

std::vector<std::string> MCMCFile

Definition at line 50 of file RHat.cpp.

◆ Mean

double** Mean

Definition at line 58 of file RHat.cpp.

◆ MeanGlobal

double* MeanGlobal

Definition at line 61 of file RHat.cpp.

◆ Nchains

int Nchains

Definition at line 46 of file RHat.cpp.

◆ nDraw

int nDraw

Definition at line 47 of file RHat.cpp.

◆ NThin

int NThin

Definition at line 45 of file RHat.cpp.

◆ Ntoys_filled

int* Ntoys_filled

Definition at line 43 of file RHat.cpp.

◆ Ntoys_requested

int* Ntoys_requested

Definition at line 42 of file RHat.cpp.

◆ RHat

double* RHat

Definition at line 66 of file RHat.cpp.

◆ S1_chain

double** S1_chain

Definition at line 55 of file RHat.cpp.

◆ S1_global

double* S1_global

Definition at line 53 of file RHat.cpp.

◆ S2_chain

double** S2_chain

Definition at line 56 of file RHat.cpp.

◆ S2_global

double* S2_global

Definition at line 54 of file RHat.cpp.

◆ StandardDeviation

double** StandardDeviation

Definition at line 59 of file RHat.cpp.

◆ StandardDeviationGlobal

double* StandardDeviationGlobal

Definition at line 62 of file RHat.cpp.

◆ TotToys

int TotToys

Definition at line 44 of file RHat.cpp.

◆ ValidPar

std::vector<bool> ValidPar

Definition at line 51 of file RHat.cpp.