MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
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.

[29].

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.

[29].

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:23
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;
655 delete[] BetweenChainVariance;
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.;
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;
94
95 BetweenChainVariance = 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
127
129
130 //KS: Main function
132
133 SaveResults();
134
136
137 return 0;
138}
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:30
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
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:11

◆ 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:22
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:212

◆ 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:26

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.