From 0538159b15e00a0664dc1cfa39fbc6e5149d7937 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 21 Apr 2026 13:19:51 +0200 Subject: [PATCH] [RF][RS] Unify `BayesianCalculator::GetInterval()` 1D and ND code paths For models with a single POI and no nuisance parameters, `BayesianCalculator::GetInterval()` used `ComputeIntervalUsingRooFit()`, which built the CDF via `RooAbsPdf::createCdf()` On a computation graph that contains the NLL. The `createCdf()` relies on cloning the integrand subtree and substituting the POI with a prime variable, but the POI dependency inside the NLL compiled for the new vectorizing CPU evaluation backend is not propagated through the clone. The CDF and therefore the calculated interval is then wrong. Users worked around this by calling `SetScanOfPosterior(N)` instead (see #17567). Fortunately, there was already a code path for the n-dimensional parameter case present that doesn't use `createCdf()`, but the ROOT Math integrator classes directly. This commit commit suggests to also do this for the 1D POI-only case to avoid the issue with using the compiled NLL in `createCdf()` and simplifying the code by only providing one code path for interval caluclation with integration. Add a gtest regression test that reproduces the issue scenario and checks the default-path upper limit matches the explicit scan reference. Closes #17567. --- .../inc/RooStats/BayesianCalculator.h | 2 - roofit/roostats/src/BayesianCalculator.cxx | 98 ++++++------------- roofit/roostats/test/CMakeLists.txt | 1 + .../roostats/test/testBayesianCalculator.cxx | 65 ++++++++++++ 4 files changed, 98 insertions(+), 68 deletions(-) create mode 100644 roofit/roostats/test/testBayesianCalculator.cxx diff --git a/roofit/roostats/inc/RooStats/BayesianCalculator.h b/roofit/roostats/inc/RooStats/BayesianCalculator.h index f0800cb866abf..5529a456d362a 100644 --- a/roofit/roostats/inc/RooStats/BayesianCalculator.h +++ b/roofit/roostats/inc/RooStats/BayesianCalculator.h @@ -151,8 +151,6 @@ namespace RooStats { void ComputeIntervalFromCdf(double c1, double c2) const; - void ComputeIntervalUsingRooFit(double c1, double c2) const; - void ComputeShortestInterval() const; private: diff --git a/roofit/roostats/src/BayesianCalculator.cxx b/roofit/roostats/src/BayesianCalculator.cxx index d836e39948c75..3d5d2e513231d 100644 --- a/roofit/roostats/src/BayesianCalculator.cxx +++ b/roofit/roostats/src/BayesianCalculator.cxx @@ -54,7 +54,6 @@ credible interval from the given function. #include "RooAbsReal.h" #include "RooRealVar.h" #include "RooArgSet.h" -#include "RooBrentRootFinder.h" #include "RooFormulaVar.h" #include "RooGenericPdf.h" #include "RooPlot.h" @@ -167,7 +166,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fFunctor(nll, bindParams, RooArgList() ), // functor fPriorFunc(nullptr), fLikelihood(fFunctor, nullptr, nllMinimum), // integral of exp(-nll) function - fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ), // integrator fXmin(bindParams.size()), // vector of parameters (min values) fXmax(bindParams.size()) { @@ -176,7 +174,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fLikelihood.SetPrior(fPriorFunc.get() ); } - fIntegrator.SetFunction(fLikelihood, bindParams.size() ); + // ROOT::Math::IntegratorMultiDim does not support dim = 1, so fall back to + // the 1D integrator when there are no nuisance parameters. + if (bindParams.size() == 1) { + fIntegratorOneDim = std::make_unique(ROOT::Math::IntegratorOneDim::GetType(integType)); + fIntegratorOneDim->SetFunction(fLikelihood); + } else { + fIntegrator = std::make_unique(ROOT::Math::IntegratorMultiDim::GetType(integType)); + fIntegrator->SetFunction(fLikelihood, bindParams.size()); + } ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. " << " nllMinimum is " << nllMinimum << std::endl; @@ -191,7 +197,10 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { << " in interval [ " << fXmin[i] << " , " << fXmax[i] << " ] " << std::endl; } - fIntegrator.Options().Print(ooccoutD(nullptr,NumericIntegration)); + if (fIntegrator) + fIntegrator->Options().Print(ooccoutD(nullptr,NumericIntegration)); + else + fIntegratorOneDim->Options().Print(ooccoutD(nullptr,NumericIntegration)); // store max POI value because it will be changed when evaluating the function fMaxPOI = fXmax[0]; @@ -214,7 +223,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { //fPriorFunc(std::shared_ptr((RooFunctor*)0)), fPriorFunc(rhs.fPriorFunc), fLikelihood(fFunctor, fPriorFunc.get(), rhs.fLikelihood.fOffset), - fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ), // integrator fXmin( rhs.fXmin), fXmax( rhs.fXmax), fNorm( rhs.fNorm), @@ -226,7 +234,13 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fError(rhs.fError), fNormCdfValues(rhs.fNormCdfValues) { - fIntegrator.SetFunction(fLikelihood, fXmin.size() ); + if (rhs.fIntegratorOneDim) { + fIntegratorOneDim = std::make_unique(ROOT::Math::IntegratorOneDim::GetType(rhs.fIntegratorOneDim->Name().c_str())); + fIntegratorOneDim->SetFunction(fLikelihood); + } else { + fIntegrator = std::make_unique(ROOT::Math::IntegratorMultiDim::GetType(rhs.fIntegrator->Name().c_str())); + fIntegrator->SetFunction(fLikelihood, fXmin.size()); + } // need special treatment for the smart pointer // if (rhs.fPriorFunc.get() ) { // fPriorFunc = std::shared_ptr(new RooFunctor(*(rhs.fPriorFunc) ) ); @@ -277,8 +291,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { fFunctor.binding().resetNumCall(); // reset number of calls for debug - double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]); - double error = fIntegrator.Error(); + double cdf; + double error; + if (fIntegratorOneDim) { + cdf = fIntegratorOneDim->Integral(fXmin[0], fXmax[0]); + error = fIntegratorOneDim->Error(); + } else { + cdf = fIntegrator->Integral(&fXmin[0], &fXmax[0]); + error = fIntegrator->Error(); + } double normcdf = cdf/fNorm; // normalize the cdf ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , " @@ -324,7 +345,8 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction { mutable RooFunctor fFunctor; // functor binding nll mutable std::shared_ptr fPriorFunc; // functor binding the prior LikelihoodFunction fLikelihood; // likelihood function - mutable ROOT::Math::IntegratorMultiDim fIntegrator; // integrator (mutable because Integral() is not const + mutable std::unique_ptr fIntegrator; // multi-dim integrator (for >=2 dims) + mutable std::unique_ptr fIntegratorOneDim; // 1D integrator (for POI-only case) mutable std::vector fXmin; // min value of parameters (poi+nuis) - mutable std::vector fXmax; // max value of parameters (poi+nuis) - max poi changes so it is mutable double fNorm = 1.0; // normalization value (computed in constructor) @@ -1127,14 +1149,7 @@ SimpleInterval* BayesianCalculator::GetInterval() const } else { - // use integration method if there are nuisance parameters - if (!fNuisanceParameters.empty()) { - ComputeIntervalFromCdf(lowerCutOff, upperCutOff); - } - else { - // case of no nuisance - just use createCdf from roofit - ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff); - } + ComputeIntervalFromCdf(lowerCutOff, upperCutOff); // case cdf failed (scan then the posterior) if (!fValidInterval) { fNScanBins = 100; @@ -1187,55 +1202,6 @@ double BayesianCalculator::GetMode() const { //return fApproxPosterior->GetMaximumX(); } -//////////////////////////////////////////////////////////////////////////////// -/// internal function compute the interval using RooFit - -void BayesianCalculator::ComputeIntervalUsingRooFit(double lowerCutOff, double upperCutOff ) const { - - coutI(Eval) << "BayesianCalculator: Compute interval using RooFit: posteriorPdf + createCdf + RooBrentRootFinder " << std::endl; - - RooRealVar* poi = dynamic_cast(fPOI.first() ); - assert(poi); - - fValidInterval = false; - if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf(); - if (!fPosteriorPdf) return; - - std::unique_ptr cdf{fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf())}; - if (!cdf) return; - - std::unique_ptr cdf_bind{cdf->bindVars(fPOI,&fPOI)}; - if (!cdf_bind) return; - - RooBrentRootFinder brf(*cdf_bind); - brf.setTol(fBrfPrecision); // set the brf precision - - double tmpVal = poi->getVal(); // patch used because findRoot changes the value of poi - - bool ret = true; - if (lowerCutOff > 0) { - double y = lowerCutOff; - ret &= brf.findRoot(fLower,poi->getMin(),poi->getMax(),y); - } - else - fLower = poi->getMin(); - - if (upperCutOff < 1.0) { - double y=upperCutOff; - ret &= brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y); - } - else - fUpper = poi->getMax(); - if (!ret) { - coutE(Eval) << "BayesianCalculator::GetInterval " - << "Error returned from Root finder, estimated interval is not fully correct" << std::endl; - } else { - fValidInterval = true; - } - - poi->setVal(tmpVal); // patch: restore the original value of poi -} - //////////////////////////////////////////////////////////////////////////////// /// internal function compute the interval using Cdf integration diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt index 75eff086a88d6..7770dc075e93e 100644 --- a/roofit/roostats/test/CMakeLists.txt +++ b/roofit/roostats/test/CMakeLists.txt @@ -1,4 +1,5 @@ ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats) +ROOT_ADD_GTEST(testBayesianCalculator testBayesianCalculator.cxx LIBRARIES RooStats) ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx LIBRARIES RooStats COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testHypoTestInvResult_1.root) diff --git a/roofit/roostats/test/testBayesianCalculator.cxx b/roofit/roostats/test/testBayesianCalculator.cxx new file mode 100644 index 0000000000000..164ed8a1357ae --- /dev/null +++ b/roofit/roostats/test/testBayesianCalculator.cxx @@ -0,0 +1,65 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +/// Regression test for https://github.com/root-project/root/issues/17567 + +TEST(BayesianCalculator, UpperLimitNoNuisance_Issue17567) +{ + RooRandom::randomGenerator()->SetSeed(12345); + + // Background-only generation model. + RooRealVar mass("mass", "mass", 6500, 7500); + RooRealVar slopeGen("slopeGen", "slopeGen", -8e-4); + RooExponential expoGen("expoGen", "expoGen", mass, slopeGen); + std::unique_ptr data{expoGen.generate(mass, 1000)}; + + // Fit model: Gaussian signal on top of exponential background. + // All parameters except the POI ngauss are fixed so there are no nuisances. + RooRealVar mean("mean", "mean", 7024); + RooRealVar sigma("sigma", "sigma", 44); + RooGaussian gauss("gauss", "gauss", mass, mean, sigma); + RooRealVar slope("slope", "slope", -8e-4); + RooExponential expo("expo", "expo", mass, slope); + RooRealVar ngauss("ngauss", "ngauss", 10, 0, 50); + RooRealVar nexpo("nexpo", "nexpo", 1e3); + RooAddPdf pdftot("pdftot", "pdftot", RooArgList(gauss, expo), RooArgList(ngauss, nexpo)); + + RooUniform prior("prior", "prior", ngauss); + + // Run with the default interval method (numerical integration instead of + // brute-force scan). + RooStats::BayesianCalculator bc(*data, pdftot, RooArgSet(ngauss), prior); + bc.SetConfidenceLevel(0.95); + bc.SetLeftSideTailFraction(0.); + std::unique_ptr bcInterval{bc.GetInterval()}; + ASSERT_NE(bcInterval, nullptr); + const double limitDefault = bcInterval->UpperLimit(); + EXPECT_EQ(bcInterval->LowerLimit(), 0.0); + + // Reference: explicit scan with many bins gives the same answer up to + // numerical integration and scan discretization errors. + RooStats::BayesianCalculator bcScan(*data, pdftot, RooArgSet(ngauss), prior); + bcScan.SetConfidenceLevel(0.95); + bcScan.SetLeftSideTailFraction(0.); + bcScan.SetScanOfPosterior(200); + std::unique_ptr scanInterval{bcScan.GetInterval()}; + ASSERT_NE(scanInterval, nullptr); + const double limitScan = scanInterval->UpperLimit(); + + // Both limits should sit well inside the POI range (if the old bug were + // present the limit would be at or near poi->getMax() = 50, or pulled down + // by a broken CDF root) and must agree within 1%. + EXPECT_GT(limitDefault, 20.0); + EXPECT_LT(limitDefault, 35.0); + EXPECT_NEAR(limitDefault, limitScan, 0.01 * limitScan); +}