RHat.cpp File Reference

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

#include "Fitters/RHatCalculator.h"

Include dependency graph for RHat.cpp:

Go to the source code of this file.

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

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.

[36].

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.

Author

Kamil Skwarczynski

Michael Reh

Definition in file RHat.cpp.

Function Documentation

main()

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

Definition at line 29 of file RHat.cpp.

                                 {
// *******************
  SetMaCh3LoggerFormat();
  M3::Utils::MaCh3Welcome();
 
  std::vector<std::string> MCMCFile;
  int Nchains = 0;
  if (argc < 2)
  {
    MACH3LOG_ERROR("Wrong arguments");
    MACH3LOG_ERROR("./RHat NThin MCMCchain_1.root MCMCchain_2.root MCMCchain_3.root ... [how many you like]");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  int NThin = atoi(argv[1]);
 
  //KS Gelman suggests to diagnose on more than one chain
  for (int i = 2; i < argc; i++)
  {
    MCMCFile.push_back(std::string(argv[i]));
    MACH3LOG_INFO("Adding file: {}", MCMCFile.back());
    Nchains++;
  }
 
  if(Nchains == 1)
  {
    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.");
    MACH3LOG_WARN("Multiple chains are more likely to reveal multimodality and poor adaptation or mixing:");
  }
  MACH3LOG_INFO("Diagnosing {} chains", Nchains);
  auto RHatCalc = std::make_unique<RHatCalculator>(false, MCMCFile, NThin);
  RHatCalc->RunDiagnostic();
 
  return 0;
}