StatisticalUtils.h File Reference

Utility functions for statistical interpretations in MaCh3. More...

#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <cmath>
#include "Manager/Manager.h"
#include "Samples/HistogramUtils.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TROOT.h"
#include "TStyle.h"

Include dependency graph for StatisticalUtils.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Namespaces
	M3
	Main namespace for MaCh3 software.

Enumerations
enum	M3::kInfCrit { M3::kBIC , M3::kDIC , M3::kWAIC , M3::kInfCrits }
	KS: Different Information Criterion tests mostly based Gelman paper. More...

Functions
std::string	GetJeffreysScale (const double BayesFactor)
	KS: Following H. Jeffreys [22]. More...
std::string	GetDunneKaboth (const double BayesFactor)
	Convert a Bayes factor into an approximate particle-physics significance level using the Dunne–Kaboth scale. More...
double	GetSigmaValue (const int sigma)
	KS: Convert sigma from normal distribution into percentage. More...
double	GetBIC (const double llh, const int data, const int nPars)
	Get the Bayesian Information Criterion (BIC) or Schwarz information criterion (also SIC, SBC, SBIC) More...
double	GetNeffective (const int N1, const int N2)
	KS: See 14.3.10 in Numerical Recipes in C [29]. More...
void	CheckBonferoniCorrectedpValue (const std::vector< std::string > &SampleNameVec, const std::vector< double > &PValVec, const double Threshold=0.05)
	KS: For more see https://www.t2k.org/docs/technotes/429/TN429_v8#page=63. More...
double	GetAndersonDarlingTestStat (const double CumulativeData, const double CumulativeMC, const double CumulativeJoint)
int	GetNumberOfRuns (const std::vector< int > &GroupClasifier)
	KS: https://esjeevanand.uccollege.edu.in/wp-content/uploads/sites/114/2020/08/NON-PARAMTERIC-TEST-6.pdf. More...
double	GetBetaParameter (const double data, const double mc, const double w2, TestStatistic TestStat)
	KS: Calculate Beta parameter which will be different based on specified test statistic. More...
double	GetSubOptimality (const std::vector< double > &EigenValues, const int TotalTarameters)
	Based on [30]. More...
void	GetArithmetic (TH1D *const hist, double &Mean, double &Error)
	CW: Get Arithmetic mean from posterior. More...
void	GetGaussian (TH1D *&hist, TF1 *gauss, double &Mean, double &Error)
	CW: Fit Gaussian to posterior. More...
void	GetHPD (TH1D *const hist, double &Mean, double &Error, double &Error_p, double &Error_m, const double coverage=0.6827)
	Get Highest Posterior Density (HPD) More...
void	GetCredibleInterval (const std::unique_ptr< TH1D > &hist, std::unique_ptr< TH1D > &hpost_copy, const double coverage=0.6827)
	KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading. More...
void	GetCredibleIntervalSig (const std::unique_ptr< TH1D > &hist, std::unique_ptr< TH1D > &hpost_copy, const bool CredibleInSigmas, const double coverage=0.6827)
	KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading. More...
void	GetCredibleRegion (std::unique_ptr< TH2D > &hist2D, const double coverage=0.6827)
	KS: Set 2D contour within some coverage. More...
void	GetCredibleRegionSig (std::unique_ptr< TH2D > &hist2D, const bool CredibleInSigmas, const double coverage=0.6827)
	KS: Set 2D contour within some coverage. More...
double	GetIQR (TH1D *Hist)
	Interquartile Range (IQR) More...
double	ComputeKLDivergence (const std::vector< double > &Data, const std::vector< double > &MC)
	Compute the Kullback-Leibler divergence between two TH2Poly histograms. More...
double	ComputeKLDivergence (TH2Poly *DataPoly, TH2Poly *PolyMC)
	Compute the Kullback-Leibler divergence between two TH2Poly histograms. More...
double	FisherCombinedPValue (const std::vector< double > &pvalues)
	KS: Combine p-values using Fisher's method. More...
std::unique_ptr< TH1D >	GetDeltaChi2 (TH1D *posterior_probability_hist)
	Convert a posterior probability histogram into a \(\Delta\chi^2\) distribution. Using the likelihood-ratio definition: More...
void	ThinningMCMC (const std::string &FilePath, const int ThinningCut)
	Thin MCMC Chain, to save space and maintain low autocorrelations. More...
double	GetZScore (const double value, const double mean, const double stddev)
	Compute the Z-score for a given value. More...
double	GetPValueFromZScore (const double zScore)
	Compute the P-value from a given Z-score. More...
double	GetModeError (TH1D *hpost)
	Get the mode error from a TH1D. More...
void	Get2DBayesianpValue (TH2D *Histogram)
	Calculates the 2D Bayesian p-value and generates a visualization. More...
void	PassErrorToRatioPlot (TH1D *RatioHist, TH1D *Hist1, TH1D *DataHist)

Detailed Description

Utility functions for statistical interpretations in MaCh3.

Author

Kamil Skwarczynski

Definition in file StatisticalUtils.h.

Function Documentation

CheckBonferoniCorrectedpValue()

void CheckBonferoniCorrectedpValue	(	const std::vector< std::string > &	SampleNameVec,
		const std::vector< double > &	PValVec,
		const double	Threshold = 0.05
	)

KS: For more see https://www.t2k.org/docs/technotes/429/TN429_v8#page=63.

Parameters

SampleNameVec	vector of sample names
PValVec	pvalue for each sample
Threshold	pvalue accepted threshold, usually 5%

Definition at line 91 of file StatisticalUtils.cpp.

                                                                  {
// ****************
  MACH3LOG_INFO("");
  if(SampleNameVec.size() != PValVec.size())
  {
    MACH3LOG_ERROR("Size of vectors do not match");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  const size_t NumberOfStatisticalTests = SampleNameVec.size();
  //KS: 0.05 or 5% is value used by T2K.
  const double StatisticalSignificanceDown = Threshold / double(NumberOfStatisticalTests);
  const double StatisticalSignificanceUp = 1 - StatisticalSignificanceDown;
  MACH3LOG_INFO("Bonferroni-corrected statistical significance level: {:.2f}", StatisticalSignificanceDown);
 
  int Counter = 0;
  for(unsigned int i = 0; i < SampleNameVec.size(); i++)
  {
    if(  (PValVec[i] < 0.5 && PValVec[i] < StatisticalSignificanceDown) ) {
      MACH3LOG_INFO("Sample {} indicates disagreement between the model predictions and the data", SampleNameVec[i]);
      MACH3LOG_INFO("Bonferroni-corrected statistical significance level: {:.2f} p-value: {:.2f}", StatisticalSignificanceDown, PValVec[i]);
      Counter++;
    } else if( (PValVec[i] > 0.5 && PValVec[i] > StatisticalSignificanceUp) ) {
      MACH3LOG_INFO("Sample {} indicates disagreement between the model predictions and the data", SampleNameVec[i]);
      MACH3LOG_INFO("Bonferroni-corrected statistical significance level: {:.2f} p-value: {:.2f}", StatisticalSignificanceUp, PValVec[i]);
      Counter++;
    }
  }
  if(Counter == 0) {
    MACH3LOG_INFO("Every sample passed Bonferroni-corrected statistical significance level test");
  } else {
    MACH3LOG_WARN("{} samples didn't pass Bonferroni-corrected statistical significance level test", Counter);
  }
  MACH3LOG_INFO("");
}

ComputeKLDivergence() [1/2]

double ComputeKLDivergence	(	const std::vector< double > &	Data,
		const std::vector< double > &	MC
	)

Compute the Kullback-Leibler divergence between two TH2Poly histograms.

Parameters

Data	Vector of data entries.
MC	Vector of MC entries.

Returns

The Kullback-Leibler divergence value. Returns 0 if the data or MC integral is zero.

Definition at line 464 of file StatisticalUtils.cpp.

                                                        {
// ********************
  double klDivergence = 0.0;
  double DataIntegral = std::accumulate(Data.begin(), Data.end(), 0.0);
  double MCIntegral   = std::accumulate(MC.begin(), MC.end(), 0.0);
  for (size_t i = 0; i < Data.size(); ++i)
  {
    if (Data[i] > 0 && MC[i] > 0) {
      klDivergence += Data[i] / DataIntegral *
                      std::log((Data[i] / DataIntegral) / ( MC[i] / MCIntegral));
    }
  }
  return klDivergence;
}

ComputeKLDivergence() [2/2]

double ComputeKLDivergence	(	TH2Poly *	DataPoly,
		TH2Poly *	PolyMC
	)

Compute the Kullback-Leibler divergence between two TH2Poly histograms.

Parameters

DataPoly	Pointer to the data histogram (TH2Poly).
PolyMC	Pointer to the Monte Carlo histogram (TH2Poly).

Returns

The Kullback-Leibler divergence value. Returns 0 if the data or MC integral is zero.

Definition at line 481 of file StatisticalUtils.cpp.

                                                               {
// *********************
  int nBins = DataPoly->GetNumberOfBins();
  std::vector<double> Data(nBins);
  std::vector<double> MC(nBins);
 
  for (int i = 0; i < nBins; ++i) {
    Data[i] = DataPoly->GetBinContent(i+1);
    MC[i] = PolyMC->GetBinContent(i+1);
  }
 
  return ComputeKLDivergence(Data, MC);
}

FisherCombinedPValue()

double FisherCombinedPValue

(

const std::vector< double > &

pvalues

)

KS: Combine p-values using Fisher's method.

Parameters

pvalues

A vector of individual p-values to combine.

Returns

The combined p-value, representing the overall significance.

Definition at line 496 of file StatisticalUtils.cpp.

                                                              {
// ********************
  double testStatistic = 0;
  for(size_t i = 0; i < pvalues.size(); i++)
  {
    const double pval = std::max(0.00001, pvalues[i]);
    testStatistic += -2.0 * std::log(pval);
  }
  // Degrees of freedom is twice the number of p-values
  int degreesOfFreedom = int(2 * pvalues.size());
  double pValue = TMath::Prob(testStatistic, degreesOfFreedom);
 
  return pValue;
}

Get2DBayesianpValue()

void Get2DBayesianpValue

(

TH2D *

Histogram

)

Calculates the 2D Bayesian p-value and generates a visualization.

Parameters

Histogram

A pointer to a TH2D histogram object. The function modifies the histogram's title to include the p-value information.

Warning

The canvas is saved to the current ROOT file using TempCanvas->Write().

Definition at line 624 of file StatisticalUtils.cpp.

                                          {
// ****************
  const double TotalIntegral = Histogram->Integral();
  // Count how many fills are above y=x axis
  // This is the 2D p-value
  double Above = 0.0;
  for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) {
    const double xvalue = Histogram->GetXaxis()->GetBinCenter(i+1);
    for (int j = 0; j < Histogram->GetYaxis()->GetNbins(); ++j) {
      const double yvalue = Histogram->GetYaxis()->GetBinCenter(j+1);
      // We're only interested in being _ABOVE_ the y=x axis
      if (xvalue >= yvalue) {
        Above += Histogram->GetBinContent(i+1, j+1);
      }
    }
  }
  const double pvalue = Above/TotalIntegral;
  std::stringstream ss;
  ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue;
  Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str());
 
  // Now add the TLine going diagonally
  double minimum = Histogram->GetXaxis()->GetBinLowEdge(1);
  if (Histogram->GetYaxis()->GetBinLowEdge(1) > minimum) {
    minimum = Histogram->GetYaxis()->GetBinLowEdge(1);
  }
  double maximum = Histogram->GetXaxis()->GetBinLowEdge(Histogram->GetXaxis()->GetNbins());
  if (Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins()) < maximum) {
    maximum = Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins());
    //KS: Extend by bin width to perfectly fit canvas
    maximum += Histogram->GetYaxis()->GetBinWidth(Histogram->GetYaxis()->GetNbins());
  }
  else maximum += Histogram->GetXaxis()->GetBinWidth(Histogram->GetXaxis()->GetNbins());
  auto TempLine = std::make_unique<TLine>(minimum, minimum, maximum, maximum);
  TempLine->SetLineColor(kRed);
  TempLine->SetLineWidth(2);
 
  std::string Title = Histogram->GetName();
  Title += "_canv";
  auto TempCanvas = std::make_unique<TCanvas>(Title.c_str(), Title.c_str(), 1024, 1024);
  TempCanvas->SetGridx();
  TempCanvas->SetGridy();
  TempCanvas->cd();
  gStyle->SetPalette(81);
  Histogram->SetMinimum(-0.01);
  Histogram->Draw("colz");
  TempLine->Draw("same");
 
  TempCanvas->Write();
}

GetAndersonDarlingTestStat()

double GetAndersonDarlingTestStat	(	const double	CumulativeData,
		const double	CumulativeMC,
		const double	CumulativeJoint
	)

Parameters

CumulativeData	Value of CDF for data
CumulativeMC	Value of CDF for MC
CumulativeJoint	Value of CDF for joint data and MC distribution

Definition at line 129 of file StatisticalUtils.cpp.

                                                                                                                        {
// ****************
  double ADstat = std::fabs(CumulativeData - CumulativeMC)/ std::sqrt(CumulativeJoint*(1 - CumulativeJoint));
 
  if( std::isinf(ADstat) || std::isnan(ADstat)) return 0;
  return ADstat;
}

GetArithmetic()

void GetArithmetic	(	TH1D *const	hist,
		double &	Mean,
		double &	Error
	)

CW: Get Arithmetic mean from posterior.

Parameters

hist	histograms from which we extract arithmetic mean
Mean	Arithmetic Mean value
Error	Arithmetic Error value

Definition at line 202 of file StatisticalUtils.cpp.

                                                                   {
// **************************
  Mean = hist->GetMean();
  Error = hist->GetRMS();
}

GetBetaParameter()

double GetBetaParameter	(	const double	data,
		const double	mc,
		const double	w2,
		TestStatistic	TestStat
	)

KS: Calculate Beta parameter which will be different based on specified test statistic.

Parameters

data	Number of data events in a bin
mc	Number of MC events in a bin
w2	Value of weight squared in a bin
TestStat	Test statistic based on which we calculate beta

Definition at line 155 of file StatisticalUtils.cpp.

                                                                                                           {
// ****************
  double Beta = 0.0;
 
  if (TestStat == kDembinskiAbdelmotteleb) {
    //the so-called effective count
    const double k = mc*mc / w2;
    //Calculate beta which is scaling factor between true and generated MC
    Beta = (data + k) / (mc + k);
  }
  //KS: Below is technically only true for Cowan's BB, which will not be true for Poisson or IceCube, because why not...
  else {
    // CW: Barlow-Beeston uses fractional uncertainty on MC, so sqrt(sum[w^2])/mc
    const double fractional = std::sqrt(w2)/mc;
    // CW: -b/2a in quadratic equation
    const double temp = mc*fractional*fractional-1;
    // CW: b^2 - 4ac in quadratic equation
    const double temp2 = temp*temp + 4*data*fractional*fractional;
    if (temp2 < 0) {
      MACH3LOG_ERROR("Negative square root in Barlow Beeston coefficient calculation!");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    // CW: Solve for the positive beta
    Beta = (-1*temp+std::sqrt(temp2))/2.;
  }
  return Beta;
}

GetBIC()

double GetBIC	(	const double	llh,
		const int	data,
		const int	nPars
	)

Get the Bayesian Information Criterion (BIC) or Schwarz information criterion (also SIC, SBC, SBIC)

Definition at line 69 of file StatisticalUtils.cpp.

                                                                {
// ****************
  if(nPars == 0)
  {
    MACH3LOG_ERROR("You haven't passed number of model parameters as it is still zero");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  const double BIC = double(nPars * logl(data) + llh);
 
  return BIC;
}

GetCredibleInterval()

void GetCredibleInterval	(	const std::unique_ptr< TH1D > &	hist,
		std::unique_ptr< TH1D > &	hpost_copy,
		const double	coverage = 0.6827
	)

KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading.

Parameters

hist	histograms based on which we calculate credible interval
hpost_copy	To make code thread safe we use copy of histograms which user has to pass
coverage	What is defined coverage, by default 0.6827 (1 sigma)

Loop over histogram bins with highest number of entries until covered 90 or 68.3%

Get bin of highest content and save the number of entries reached so far

Replace bin value by -1 so it is not looped over as being maximum bin again

Definition at line 324 of file StatisticalUtils.cpp.

                                                                                                                  {
// ***************
  if(coverage > 1)
  {
    MACH3LOG_ERROR("Specified Credible Interval is greater that 1 and equal to {} Should be between 0 and 1", coverage);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  //KS: Reset first copy of histogram
  hpost_copy->Reset("");
  hpost_copy->Fill(0.0, 0.0);
 
  //KS: Temporary structure to be thread save
  std::vector<double> hist_copy(hist->GetXaxis()->GetNbins()+1);
  std::vector<bool> hist_copy_fill(hist->GetXaxis()->GetNbins()+1);
  for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i)
  {
    hist_copy[i] = hist->GetBinContent(i);
    hist_copy_fill[i] = false;
  }
 
  const long double Integral = hist->Integral();
  long double sum = 0;
 
  while ((sum / Integral) < coverage)
  {
    int max_entry_bin = 0;
    double max_entries = 0.;
    for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i)
    {
      if (hist_copy[i] > max_entries)
      {
        max_entries = hist_copy[i];
        max_entry_bin = i;
      }
    }
    hist_copy[max_entry_bin] = -1.;
    hist_copy_fill[max_entry_bin] = true;
 
    sum += max_entries;
  }
  //KS: Now fill our copy only for bins which got included in coverage region
  for(int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i)
  {
    if(hist_copy_fill[i]) hpost_copy->SetBinContent(i, hist->GetBinContent(i));
  }
}

GetCredibleIntervalSig()

void GetCredibleIntervalSig	(	const std::unique_ptr< TH1D > &	hist,
		std::unique_ptr< TH1D > &	hpost_copy,
		const bool	CredibleInSigmas,
		const double	coverage = 0.6827
	)

KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading.

Parameters

hist	histograms based on which we calculate credible interval
hpost_copy	To make code thread safe we use copy of histograms which user has to pass
CredibleInSigmas	Whether interval is in sigmas or percentage
coverage	What is defined coverage, by default 0.6827 (1 sigma)

Definition at line 375 of file StatisticalUtils.cpp.

                                                                                                                                                  {
// ***************
  //KS: Slightly different approach depending if intervals are in percentage or sigmas
  if(CredibleInSigmas) {
    //KS: Convert sigmas into percentage
    const double CredReg = GetSigmaValue(int(std::round(coverage)));
    GetCredibleInterval(hist, hpost_copy, CredReg);
  } else {
    GetCredibleInterval(hist, hpost_copy, coverage);
  }
}

GetCredibleRegion()

void GetCredibleRegion	(	std::unique_ptr< TH2D > &	hist2D,
		const double	coverage = 0.6827
	)

KS: Set 2D contour within some coverage.

Parameters

hist2D	histograms based on which we calculate credible regions
coverage	What is defined coverage, by default 0.6827 (1 sigma)

Loop over histogram bins with highest number of entries until covered 90 or 68.3%

Get bin of highest content and save the number of entries reached so far

Replace bin value by -1 so it is not looped over as being maximum bin again

Definition at line 388 of file StatisticalUtils.cpp.

                                                                           {
// ***************
  if(coverage > 1)
  {
    MACH3LOG_ERROR("Specified Credible Region is greater than 1 and equal to {:.2f} Should be between 0 and 1", coverage);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  //KS: Temporary structure to be thread save
  std::vector<std::vector<double>> hist_copy(hist2D->GetXaxis()->GetNbins()+1,
                                             std::vector<double>(hist2D->GetYaxis()->GetNbins()+1));
  for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) {
    for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j) {
      hist_copy[i][j] = hist2D->GetBinContent(i, j);
    }
  }
 
  const long double Integral = hist2D->Integral();
  long double sum = 0;
 
  //We need to as ROOT requires array to set to contour
  double Contour[1];
  while ((sum / Integral) < coverage)
  {
    int max_entry_bin_x = 0;
    int max_entry_bin_y = 0;
    double max_entries = 0.;
    for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i)
    {
      for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j)
      {
        if (hist_copy[i][j] > max_entries)
        {
          max_entries = hist_copy[i][j];
          max_entry_bin_x = i;
          max_entry_bin_y = j;
        }
      }
    }
    hist_copy[max_entry_bin_x][max_entry_bin_y] = -1.;
 
    sum += max_entries;
    Contour[0] = max_entries;
  }
  hist2D->SetContour(1, Contour);
}

GetCredibleRegionSig()

void GetCredibleRegionSig	(	std::unique_ptr< TH2D > &	hist2D,
		const bool	CredibleInSigmas,
		const double	coverage = 0.6827
	)

KS: Set 2D contour within some coverage.

Parameters

hist2D	histograms based on which we calculate credible regions
CredibleInSigmas	Whether interval is in sigmas or percentage
coverage	What is defined coverage, by default 0.6827 (1 sigma)

Definition at line 439 of file StatisticalUtils.cpp.

                                                                                                           {
// ***************
  if(CredibleInSigmas) {
    //KS: Convert sigmas into percentage
    const double CredReg = GetSigmaValue(int(std::round(coverage)));
    GetCredibleRegion(hist2D, CredReg);
  } else {
    GetCredibleRegion(hist2D, coverage);
  }
}

GetDeltaChi2()

std::unique_ptr<TH1D> GetDeltaChi2

(

TH1D *

posterior_probability_hist

)

Convert a posterior probability histogram into a \(\Delta\chi^2\) distribution. Using the likelihood-ratio definition:

\[ \Delta\chi^2 = -2 \ln\left(\frac{L}{L_{\max}}\right) \]

Parameters

posterior_probability_hist

Pointer to a TH1D histogram containing posterior probabilities

Note

based on CompareMaCh3PThetaDeltaChi2.C

Definition at line 677 of file StatisticalUtils.cpp.

                                                                   {
// ****************
  auto delta_chi2 = M3::Clone(posterior_probability_hist);
  delta_chi2->GetYaxis()->SetTitle("#Delta#chi^{2}");
 
  int max_bin = delta_chi2->GetMaximumBin();
  double max_content = delta_chi2->GetBinContent(max_bin);
  if (max_content == 0) {
    MACH3LOG_ERROR("Histogram {}, has larges bin with 0", delta_chi2->GetTitle());
    MACH3LOG_ERROR("This suggest you skewed binning for posterior probability or something else");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  double NewMaximum = M3::_BAD_DOUBLE_;
  for(int iBin = 1; iBin < delta_chi2->GetNbinsX()+1; iBin++) {
    double bin_content = delta_chi2->GetBinContent(iBin);
    if(bin_content == 0) bin_content = 1.0 ;
 
    double chi2_likelihood = -2*std::log(bin_content/max_content);
    delta_chi2->SetBinContent(iBin, chi2_likelihood);
    NewMaximum = std::max(NewMaximum, chi2_likelihood);
  }
  delta_chi2->SetMaximum(NewMaximum*1.1);
  return delta_chi2;
}

GetDunneKaboth()

std::string GetDunneKaboth

(

const double

BayesFactor

)

Convert a Bayes factor into an approximate particle-physics significance level using the Dunne–Kaboth scale.

This function maps a Bayes factor B(θ1,θ2) to an equivalent Gaussian significance ("n σ") assuming equal priors for the two hypotheses.

The thresholds are based on commonly used particle-physics probability levels and their corresponding Bayes factors:

Bayes factor (B) Approximate significance

B < 2.125 < 1 σ 2.125 ≤ B < 20.74 > 1 σ 20.74 ≤ B < 369.4 > 2 σ 369.4 ≤ B < 15800 > 3 σ 15800 ≤ B < 1745000 > 4 σ B ≥ 1745000 > 5 σ

For example, a Bayes factor of 369.4 corresponds approximately to a 3 σ effect in traditional particle-physics terminology.

Note

For T2K users see Table 1 in https://www.t2k.org/docs/technotes/435

Parameters

BayesFactor

The Bayes factor B(θ1,θ2).

Returns

A string describing the equivalent significance level.

Definition at line 21 of file StatisticalUtils.cpp.

                                                  {
// **************************
  std::string DunneKaboth = "";
  //KS: Get fancy DunneKaboth Scale as I am to lazy to look into table every time
 
  if(2.125 > BayesFactor)         DunneKaboth = "< 1 #sigma";
  else if( 20.74 > BayesFactor)   DunneKaboth = "> 1 #sigma";
  else if( 369.4 > BayesFactor)   DunneKaboth = "> 2 #sigma";
  else if( 15800 > BayesFactor)   DunneKaboth = "> 3 #sigma";
  else if( 1745000 > BayesFactor) DunneKaboth = "> 4 #sigma";
  else DunneKaboth = "> 5 #sigma";
 
  return DunneKaboth;
}

GetGaussian()

void GetGaussian	(	TH1D *&	hist,
		TF1 *	gauss,
		double &	Mean,
		double &	Error
	)

CW: Fit Gaussian to posterior.

Parameters

hist	histograms to which we fit gaussian
gauss	tf1 with gaussian, we pass pointer to make things faster
Mean	Gaussian Mean value
Error	Gaussian Error value

Definition at line 209 of file StatisticalUtils.cpp.

                                                                       {
// **************************
  // Supress spammy ROOT messages
  int originalErrorLevel = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kFatal;
 
  const double meanval = hist->GetMean();
  const double err = hist->GetRMS();
  const double peakval = hist->GetBinCenter(hist->GetMaximumBin());
 
  // Set the range for the Gaussian fit
  gauss->SetRange(meanval - 1.5*err , meanval + 1.5*err);
  // Set the starting parameters close to RMS and peaks of the histograms
  gauss->SetParameters(hist->GetMaximum()*err*std::sqrt(2*3.14), peakval, err);
 
  // Perform the fit
  hist->Fit(gauss->GetName(),"Rq");
  hist->SetStats(0);
 
  Mean = gauss->GetParameter(1);
  Error = gauss->GetParameter(2);
 
  // restore original warning setting
  gErrorIgnoreLevel = originalErrorLevel;
}

GetHPD()

void GetHPD	(	TH1D *const	hist,
		double &	Mean,
		double &	Error,
		double &	Error_p,
		double &	Error_m,
		const double	coverage = 0.6827
	)

Get Highest Posterior Density (HPD)

Parameters

hist	histograms from which we HPD
Mean	HPD Mean value
Error	HPD Error value
Error_p	HPD Negative (left hand side) Error value
Error_m	HPD Positive (right hand side) Error value
coverage	What is defined coverage, by default 0.6827 (1 sigma)

Definition at line 236 of file StatisticalUtils.cpp.

                                                                                                                    {
// ****************
  // Get the bin which has the largest posterior density
  const int MaxBin = hist->GetMaximumBin();
  // And it's value
  const double peakval = hist->GetBinCenter(MaxBin);
 
  // The total integral of the posterior
  const long double Integral = hist->Integral();
  //KS: and integral of left handed and right handed parts
  const long double LowIntegral = hist->Integral(1, MaxBin-1) + hist->GetBinContent(MaxBin)/2.0;
  const long double HighIntegral = hist->Integral(MaxBin+1, hist->GetNbinsX()) + hist->GetBinContent(MaxBin)/2.0;
 
  // Keep count of how much area we're covering
  //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error
  long double sum = hist->GetBinContent(MaxBin)/2.0;
 
  // Counter for current bin
  int CurrBin = MaxBin;
  while (sum/HighIntegral < coverage && CurrBin < hist->GetNbinsX()) {
    CurrBin++;
    sum += hist->GetBinContent(CurrBin);
  }
  const double sigma_p = std::fabs(hist->GetBinCenter(MaxBin)-hist->GetXaxis()->GetBinUpEdge(CurrBin));
  // Reset the sum
  //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error
  sum = hist->GetBinContent(MaxBin)/2.0;
 
  // Reset the bin counter
  CurrBin = MaxBin;
  // Counter for current bin
  while (sum/LowIntegral < coverage && CurrBin > 1) {
    CurrBin--;
    sum += hist->GetBinContent(CurrBin);
  }
  const double sigma_m = std::fabs(hist->GetBinCenter(CurrBin)-hist->GetBinLowEdge(MaxBin));
 
  // Now do the double sided HPD
  //KS: Start sum from the HPD
  sum = hist->GetBinContent(MaxBin);
  int LowBin = MaxBin;
  int HighBin = MaxBin;
  long double LowCon = 0.0;
  long double HighCon = 0.0;
 
  while (sum/Integral < coverage && (LowBin > 0 || HighBin < hist->GetNbinsX()+1))
  {
    LowCon = 0.0;
    HighCon = 0.0;
    //KS:: Move further only if you haven't reached histogram end
    if(LowBin > 1)
    {
      LowBin--;
      LowCon = hist->GetBinContent(LowBin);
    }
    if(HighBin < hist->GetNbinsX())
    {
      HighBin++;
      HighCon = hist->GetBinContent(HighBin);
    }
 
    // If we're on the last slice and the lower contour is larger than the upper
    if ((sum+LowCon+HighCon)/Integral > coverage && LowCon > HighCon) {
      sum += LowCon;
      break;
      // If we're on the last slice and the upper contour is larger than the lower
    } else if ((sum+LowCon+HighCon)/Integral > coverage && HighCon >= LowCon) {
      sum += HighCon;
      break;
    } else {
      sum += LowCon + HighCon;
    }
  }
 
  double sigma_hpd = 0.0;
  if (LowCon > HighCon) {
    sigma_hpd = std::fabs(hist->GetBinLowEdge(LowBin)-hist->GetBinCenter(MaxBin));
  } else {
    sigma_hpd = std::fabs(hist->GetXaxis()->GetBinUpEdge(HighBin)-hist->GetBinCenter(MaxBin));
  }
 
  Mean = peakval;
  Error = sigma_hpd;
  Error_p = sigma_p;
  Error_m = sigma_m;
}

GetIQR()

double GetIQR

(

TH1D *

Hist

)

Interquartile Range (IQR)

Parameters

Hist	histograms from which we IQR

Definition at line 451 of file StatisticalUtils.cpp.

                          {
// *********************
  if(Hist->Integral() == 0) return 0.0;
 
  constexpr double quartiles_x[3] = {0.25, 0.5, 0.75};
  double quartiles[3];
 
  Hist->GetQuantiles(3, quartiles, quartiles_x);
 
  return quartiles[2] - quartiles[0];
}

GetJeffreysScale()

std::string GetJeffreysScale

(

const double

BayesFactor

)

KS: Following H. Jeffreys [22].

Parameters

BayesFactor

Obtained value of Bayes factor

Definition at line 6 of file StatisticalUtils.cpp.

                                                    {
// **************************
  std::string JeffreysScale = "";
  //KS: Get fancy Jeffreys Scale as I am to lazy to look into table every time
  if(BayesFactor < 0)        JeffreysScale = "Negative";
  else if( 5 > BayesFactor)  JeffreysScale = "Barely worth mentioning";
  else if( 10 > BayesFactor) JeffreysScale = "Substantial";
  else if( 15 > BayesFactor) JeffreysScale = "Strong";
  else if( 20 > BayesFactor) JeffreysScale = "Very strong";
  else JeffreysScale = "Decisive";
 
  return JeffreysScale;
}

GetModeError()

double GetModeError

(

TH1D *

hpost

)

Get the mode error from a TH1D.

Parameters

hpost

hist from which we extract mode error

Definition at line 570 of file StatisticalUtils.cpp.

                                 {
// ****************
  // Get the bin which has the largest posterior density
  int MaxBin = hpost->GetMaximumBin();
 
  // The total integral of the posterior
  const double Integral = hpost->Integral();
 
  int LowBin = MaxBin;
  int HighBin = MaxBin;
  double sum = hpost->GetBinContent(MaxBin);;
  double LowCon = 0.0;
  double HighCon = 0.0;
  while (sum/Integral < 0.6827 && (LowBin > 0 || HighBin < hpost->GetNbinsX()+1) )
  {
    LowCon = 0.0;
    HighCon = 0.0;
    //KS:: Move further only if you haven't reached histogram end
    if(LowBin > 1)
    {
      LowBin--;
      LowCon = hpost->GetBinContent(LowBin);
    }
    if(HighBin < hpost->GetNbinsX())
    {
      HighBin++;
      HighCon = hpost->GetBinContent(HighBin);
    }
 
    // If we're on the last slice and the lower contour is larger than the upper
    if ((sum+LowCon+HighCon)/Integral > 0.6827 && LowCon > HighCon) {
      sum += LowCon;
      break;
      // If we're on the last slice and the upper contour is larger than the lower
    } else if ((sum+LowCon+HighCon)/Integral > 0.6827 && HighCon >= LowCon) {
      sum += HighCon;
      break;
    } else {
      sum += LowCon + HighCon;
    }
  }
 
  double Mode_Error = 0.0;
  if (LowCon > HighCon) {
    Mode_Error = std::fabs(hpost->GetBinCenter(LowBin)-hpost->GetBinCenter(MaxBin));
  } else {
    Mode_Error = std::fabs(hpost->GetBinCenter(HighBin)-hpost->GetBinCenter(MaxBin));
  }
 
  return Mode_Error;
}

GetNeffective()

double GetNeffective	(	const int	N1,
		const int	N2
	)

KS: See 14.3.10 in Numerical Recipes in C [29].

Definition at line 82 of file StatisticalUtils.cpp.

                                                 {
// ****************
  const double Nominator = (N1+N2);
  const double Denominator = (N1*N2);
  const double N_e = Nominator/Denominator;
  return N_e;
}

GetNumberOfRuns()

int GetNumberOfRuns

(

const std::vector< int > &

GroupClasifier

)

KS: https://esjeevanand.uccollege.edu.in/wp-content/uploads/sites/114/2020/08/NON-PARAMTERIC-TEST-6.pdf.

Definition at line 138 of file StatisticalUtils.cpp.

                                                          {
// ****************
  int NumberOfRuns = 0;
  int PreviousGroup = -999;
 
  //KS: If group changed increment run
  for (unsigned int i = 0; i < GroupClasifier.size(); i++)
  {
    if(GroupClasifier[i] != PreviousGroup)
      NumberOfRuns++;
    PreviousGroup = GroupClasifier[i];
  }
 
  return NumberOfRuns;
}

GetPValueFromZScore()

double GetPValueFromZScore

(

const double

zScore

)

Compute the P-value from a given Z-score.

The P-value represents the probability of observing a value as extreme as the given Z-score under the null hypothesis of a standard normal distribution.

Parameters

zScore

The Z-score for which to compute the P-value.

Returns

The P-value corresponding to the given Z-score.

Definition at line 563 of file StatisticalUtils.cpp.

                                                {
// ********************
  return 0.5 * std::erfc(-zScore / std::sqrt(2));
}

GetSigmaValue()

double GetSigmaValue

(

const int

sigma

)

KS: Convert sigma from normal distribution into percentage.

Definition at line 37 of file StatisticalUtils.cpp.

                                      {
// *********************
  double width = 0;
  switch (std::abs(sigma))
  {
    case 1:
      width = 0.682689492137;
      break;
    case 2:
      width = 0.954499736104;
      break;
    case 3:
      width = 0.997300203937;
      break;
    case 4:
      width = 0.999936657516;
      break;
    case 5:
      width = 0.999999426697;
      break;
    case 6:
      width = 0.999999998027;
      break;
    default:
      MACH3LOG_ERROR("{}  is unsupported value of sigma", sigma);
      throw MaCh3Exception(__FILE__ , __LINE__ );
      break;
  }
  return width;
}

GetSubOptimality()

double GetSubOptimality	(	const std::vector< double > &	EigenValues,
		const int	TotalTarameters
	)

Based on [30].

Parameters

EigenValues	Eigen values of covariance matrix
TotalTarameters	Total number of parameters needed to correctly calculate suboptimality

Definition at line 185 of file StatisticalUtils.cpp.

                                                                                         {
// *********************
  double sum_eigenvalues_squared_inv = 0.0;
  double sum_eigenvalues_inv = 0.0;
  for (unsigned int j = 0; j < EigenValues.size(); j++)
  {
    //KS: IF Eigen values are super small skip them
    //if(EigenValues[j] < 0.0000001) continue;
    sum_eigenvalues_squared_inv += std::pow(EigenValues[j], -2);
    sum_eigenvalues_inv += 1.0 / EigenValues[j];
  }
  const double SubOptimality = TotalTarameters * sum_eigenvalues_squared_inv / std::pow(sum_eigenvalues_inv, 2);
  return SubOptimality;
}

GetZScore()

double GetZScore	(	const double	value,
		const double	mean,
		const double	stddev
	)

Compute the Z-score for a given value.

The Z-score indicates how many standard deviations a value is from the mean. A positive Z-score means the value is above the mean, while a negative Z-score means it is below the mean.

Parameters

value	The data point for which to compute the Z-score.
mean	The mean of the data set.
stddev	The standard deviation of the data set. Must be non-zero.

Returns

The Z-score of the given value.

Warning

Ensure that stddev is not zero to avoid division by zero.

Definition at line 557 of file StatisticalUtils.cpp.

                                                                             {
// ********************
  return (value - mean) / stddev;
}

PassErrorToRatioPlot()

void PassErrorToRatioPlot	(	TH1D *	RatioHist,
		TH1D *	Hist1,
		TH1D *	DataHist
	)

Definition at line 704 of file StatisticalUtils.cpp.

                                                                        {
// ****************
  for (int j = 0; j <= RatioHist->GetNbinsX(); ++j)
  {
    if(DataHist->GetBinContent(j) > 0)
    {
      double dx = Hist1->GetBinError(j) / DataHist->GetBinContent(j);
      RatioHist->SetBinError(j, dx);
    }
  }
}

ThinningMCMC()

void ThinningMCMC	(	const std::string &	FilePath,
		const int	ThinningCut
	)

Thin MCMC Chain, to save space and maintain low autocorrelations.

Parameters

FilePath	Path to MCMC chain you want to thin
ThinningCut	every which entry you want to thin [24]

Warning

Thinning is done over entry not steps, it may now work very well for merged chains

Definition at line 512 of file StatisticalUtils.cpp.

                                                                    {
// ********************
  std::string FilePathNowRoot = FilePath;
  if (FilePath.size() >= 5 && FilePath.substr(FilePath.size() - 5) == ".root") {
    FilePathNowRoot = FilePath.substr(0, FilePath.size() - 5);
  }
  std::string TempFilePath = FilePathNowRoot + "_thinned.root";
  int ret = system(("cp " + FilePath + " " + TempFilePath).c_str());
  if (ret != 0) {
    MACH3LOG_WARN("System call to copy file failed with code {}", ret);
  }
 
  TFile *inFile = M3::Open(TempFilePath, "RECREATE", __FILE__, __LINE__);
  TTree *inTree = inFile->Get<TTree>("posteriors");
  if (!inTree) {
    MACH3LOG_ERROR("TTree 'posteriors' not found in file.");
    inFile->ls();
    inFile->Close();
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  // Clone the structure without data
  TTree *outTree = inTree->CloneTree(0);
 
  // Loop over entries and apply thinning
  Long64_t nEntries = inTree->GetEntries();
  double retainedPercentage = (double(nEntries) / ThinningCut) / double(nEntries) * 100;
  MACH3LOG_INFO("Thinning will retain {:.2f}% of chains", retainedPercentage);
  for (Long64_t i = 0; i < nEntries; i++) {
    if (i % (nEntries/10) == 0) {
      M3::Utils::PrintProgressBar(i, nEntries);
    }
    if (i % ThinningCut == 0) {
      inTree->GetEntry(i);
      outTree->Fill();
    }
  }
  inFile->WriteTObject(outTree, "posteriors", "kOverwrite");
  inFile->Close();
  delete inFile;
 
  MACH3LOG_INFO("Thinned TTree saved and overwrote original in: {}", TempFilePath);
}