Skip to content

feat: posterior to Chi2 conversion #911

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
KSkwarczynski merged 4 commits into from
May 20, 2026
Merged

feat: posterior to Chi2 conversion #911

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b138889
Implement posterior to Chi2 conversion
KSkwarczynski May 19, 2026
92414b3
adress comment by Henry
KSkwarczynski May 19, 2026
b4e3440
more tweaks after Henry
KSkwarczynski May 19, 2026
13cea0e
missed header
KSkwarczynski May 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Diagnostics/ProcessMCMC.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,7 @@ void ProcessMCMC(const std::string& inputFile)
//KS: When creating covariance matrix longest time is spend on caching every step, since we already cached we can run some fancy covariance stability diagnostic
if(GetFromManager<bool>(Settings["DiagnoseCovarianceMatrix"], false)) DiagnoseCovarianceMatrix(Processor.get(), inputFile);
}
Processor->ProduceChi2(GetFromManager<std::string>(Settings["Chi2Group"], "Osc"));
if(GetFromManager<bool>(Settings["JarlskogAnalysis"], true)) Processor->PerformJarlskogAnalysis();
if(GetFromManager<bool>(Settings["ProducePMNSElements"], true)) Processor->ProducePMNSElements();
if(GetFromManager<bool>(Settings["ProduceUnitarityTriangles"], true)) Processor->ProduceUnitarityTriangles();
Expand Down
52 changes: 42 additions & 10 deletions Fitters/MCMCProcessor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,7 @@ void MCMCProcessor::MakePostfit(const std::map<std::string, std::pair<double, do
if (hpost[i]->GetMaximum() == hpost[i]->Integral()*DrawRange)
{
MACH3LOG_WARN("Found fixed parameter: {} ({}), moving on", Title, i);
IamVaried[i] = false;
ParamVaried[i] = false;
//KS:Set mean and error to prior for fixed parameters, it looks much better when fixed parameter has mean on prior rather than on 0 with 0 error.
(*Means_HPD)(i) = Prior;
(*Errors_HPD)(i) = PriorError;
Expand All @@ -360,7 +360,7 @@ void MCMCProcessor::MakePostfit(const std::map<std::string, std::pair<double, do
}

// Store that this parameter is indeed being varied
IamVaried[i] = true;
ParamVaried[i] = true;

// Draw onto the TCanvas
hpost[i]->Draw();
Expand Down Expand Up @@ -746,7 +746,7 @@ void MCMCProcessor::MakeCredibleIntervals(const std::vector<double>& CredibleInt

for (int i = 0; i < nDraw; ++i)
{
if(!IamVaried[i]) continue;
if(!ParamVaried[i]) continue;

// Now make the TLine for the Asimov
TString Title = "";
Expand Down Expand Up @@ -854,7 +854,7 @@ void MCMCProcessor::MakeViolin() {
for (int x = 0; x < nDraw; ++x)
{
//KS: Consider another treatment for fixed params
//if (IamVaried[x] == false) continue;
//if (ParamVaried[x] == false) continue;
for (int k = 0; k < nEntries; ++k)
{
//KS: Burn in cut
Expand Down Expand Up @@ -985,7 +985,7 @@ void MCMCProcessor::MakeCovariance() {
if (j == i) continue;

// If this parameter isn't varied
if (IamVaried[j] == false) {
if (ParamVaried[j] == false) {
(*Covariance)(i,j) = 0.0;
(*Covariance)(j,i) = (*Covariance)(i,j);
(*Correlation)(i,j) = 0.0;
Expand Down Expand Up @@ -1228,7 +1228,7 @@ void MCMCProcessor::MakeCovariance_MP(const bool Mute) {
if (j == i) continue;

// If this parameter isn't varied
if (IamVaried[j] == false) {
if (ParamVaried[j] == false) {
(*Covariance)(i,j) = 0.0;
(*Covariance)(j,i) = (*Covariance)(i,j);
(*Correlation)(i,j) = 0.0;
Expand Down Expand Up @@ -1271,7 +1271,7 @@ void MCMCProcessor::MakeCovariance_MP(const bool Mute) {
{
// Skip the diagonal elements which we've already done above
if (j == i) continue;
if (IamVaried[j] == false) continue;
if (ParamVaried[j] == false) continue;

if(ParamType[i] == kXSecPar && ParamType[j] == kXSecPar)
{
Expand Down Expand Up @@ -1707,7 +1707,7 @@ void MCMCProcessor::DrawCorrelations1D() {

for(int i = 0; i < nDraw; i++)
{
if (IamVaried[i] == false) continue;
if (ParamVaried[i] == false) continue;

Corr1DHist[i][0]->GetXaxis()->LabelsOption("v");
Corr1DHist[i][0]->SetMaximum(+1.);
Expand Down Expand Up @@ -1764,6 +1764,38 @@ void MCMCProcessor::DrawCorrelations1D() {
gStyle->SetOptTitle(OptTitle);
}


// *********************
// Convert posterior likelihood to Delta Chi2 used for comparison with frequentists fitter
void MCMCProcessor::ProduceChi2(const std::string& GroupName) const {
// *********************
if(GroupName == "") return;
MACH3LOG_INFO("Starting {}", __func__);
TDirectory* Chi2Folder = OutputFile->mkdir("DeltaChi2");
Comment thread
henry-wallace-phys marked this conversation as resolved.

Chi2Folder->cd();
for (int iPar = 0; iPar < nDraw; iPar++)
{
std::string GroupNameCurr;
if(ParamType[iPar] == kXSecPar){
const int InternalNumeration = iPar - ParamTypeStartPos[kXSecPar];
GroupNameCurr = ParameterGroup[InternalNumeration];
} else {
GroupNameCurr = "Other"; // Use Other for all legacy params
}
if (ParamVaried[iPar] == false) continue;
if (GroupName != "All" && GroupNameCurr != GroupName) continue;

auto Chi2 = GetDeltaChi2(hpost[iPar]);
RemoveFitter(Chi2.get(), "Gauss");

Chi2->Write();
}
Chi2Folder->Close();
delete Chi2Folder;
OutputFile->cd();
}

// *********************
// Make fancy Credible Intervals plots
void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegions,
Expand Down Expand Up @@ -1822,7 +1854,7 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
{
// Skip the diagonal elements which we've already done above
if (j == i) continue;
if (IamVaried[j] == false) continue;
if (ParamVaried[j] == false) continue;

auto legend = std::make_unique<TLegend>(0.20, 0.7, 0.4, 0.92);
legend->SetTextColor(kRed);
Expand Down Expand Up @@ -2280,7 +2312,7 @@ void MCMCProcessor::ScanInput() {
// Check order of parameter types
ScanParameterOrder();

IamVaried.resize(nDraw, true);
ParamVaried.resize(nDraw, true);

// Print useful Info
PrintInfo();
Expand Down
5 changes: 4 additions & 1 deletion Fitters/MCMCProcessor.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,9 @@ class MCMCProcessor {
const std::string& Title,
const double EvaluationPoint) const ;

/// @brief Convert posterior likelihood to Delta Chi2 used for comparison with frequentists fitter
void ProduceChi2(const std::string& GroupName) const;

/// @brief Reweight Prior by giving new central value and new error
/// @param Names Parameter names for which we do reweighting
/// @param NewCentral New central value for which we reweight
Expand Down Expand Up @@ -460,7 +463,7 @@ class MCMCProcessor {
std::vector<std::string> ExcludedGroups;

/// Is the ith parameter varied
std::vector<bool> IamVaried;
std::vector<bool> ParamVaried;
/// Name of parameters which we are going to analyse
std::vector<std::vector<TString>> ParamNames;
/// Parameters central values which we are going to analyse
Expand Down
26 changes: 26 additions & 0 deletions Fitters/StatisticalUtils.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -672,7 +672,33 @@ void Get2DBayesianpValue(TH2D *Histogram) {
TempCanvas->Write();
}

// ****************
// Converts posterior likelihood to dchi2
std::unique_ptr<TH1D> GetDeltaChi2(TH1D* posterior_probability_hist) {
// ****************
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);
Comment thread
henry-wallace-phys marked this conversation as resolved.
delta_chi2->SetBinContent(iBin, chi2_likelihood);
NewMaximum = std::max(NewMaximum, chi2_likelihood);
}
delta_chi2->SetMaximum(NewMaximum*1.1);
Comment thread
henry-wallace-phys marked this conversation as resolved.
return delta_chi2;
}

// ****************
void PassErrorToRatioPlot(TH1D* RatioHist, TH1D* Hist1, TH1D* DataHist) {
Expand Down
11 changes: 11 additions & 0 deletions Fitters/StatisticalUtils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,17 @@ double ComputeKLDivergence(TH2Poly* DataPoly, TH2Poly* PolyMC);
/// @return The combined p-value, representing the overall significance.
double FisherCombinedPValue(const std::vector<double>& pvalues);

/// @brief Convert a posterior probability histogram into a \f$\Delta\chi^2\f$ distribution. Using
/// the likelihood-ratio definition:
///
/// \f[
/// \Delta\chi^2 = -2 \ln\left(\frac{L}{L_{\max}}\right)
/// \f]
///
/// @param posterior_probability_hist Pointer to a TH1D histogram containing posterior probabilities
/// @note based on CompareMaCh3PThetaDeltaChi2.C
std::unique_ptr<TH1D> GetDeltaChi2(TH1D* posterior_probability_hist);

/// @brief Thin MCMC Chain, to save space and maintain low autocorrelations.
///
/// @param FilePath Path to MCMC chain you want to thin
Expand Down
Loading