From 72d935cf09208668bc58a88d0f190e2aa81e6995 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Wed, 19 Aug 2026 10:33:42 +0000 Subject: [PATCH] [RF][RS] Add SetSigned() option to the AsymptoticCalculator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Support the signed (uncapped) profile likelihood test statistic in the AsymptoticCalculator, like the ProfileLikelihoodTestStat already does with its SetSigned() option. When enabled and the best fit value of the POI is beyond the tested value, the sign of sqrt(qmu) is flipped instead of setting qmu to zero. The Gaussian asymptotic formulae remain valid in that case, so the p-values continue above 0.5 instead of saturating there. The qtilde corrections are skipped when the sign was flipped, because they apply when the best fit value is at the POI boundary, which is on the other side of the tested value. The option has no effect for the two-sided test statistics, whose asymptotic p-value formulae are symmetric in the sign. Closes [ROOT-8257](https://its.cern.ch/jira/browse/ROOT-8257). 🤖 Done with the help of AI --- .../inc/RooStats/AsymptoticCalculator.h | 8 +++ roofit/roostats/src/AsymptoticCalculator.cxx | 38 +++++++--- .../test/testAsymptoticCalculator.cxx | 72 +++++++++++++++++++ 3 files changed, 109 insertions(+), 9 deletions(-) diff --git a/roofit/roostats/inc/RooStats/AsymptoticCalculator.h b/roofit/roostats/inc/RooStats/AsymptoticCalculator.h index fd0b5e0477479..7d5e76fc2a97c 100644 --- a/roofit/roostats/inc/RooStats/AsymptoticCalculator.h +++ b/roofit/roostats/inc/RooStats/AsymptoticCalculator.h @@ -64,6 +64,13 @@ namespace RooStats { /// set the test statistics for one-sided discovery void SetOneSidedDiscovery(bool on) { fOneSidedDiscovery = on; } + /// use the signed (uncapped) profile likelihood test statistic. + /// When using a one-sided test statistic, the value is not set to zero when the best fit value of the + /// POI is beyond the tested value, but the sign of the test statistic is flipped instead + /// (see also ProfileLikelihoodTestStat::SetSigned). + /// It has no effect for the two-sided test statistics. + void SetSigned(bool on = true) { fSigned = on; } + /// re-implementation of setters since they are needed to re-initialize the calculator void SetNullModel(const ModelConfig &nullModel) override { HypoTestCalculatorGeneric::SetNullModel(nullModel); @@ -98,6 +105,7 @@ namespace RooStats { private: bool fOneSided; ///< for one sided PL test statistic (upper limits) mutable bool fOneSidedDiscovery; ///< for one sided PL test statistic (for discovery) + bool fSigned = false; ///< use signed (uncapped) PL test statistic bool fNominalAsimov; ///< make Asimov at nominal parameter values mutable bool fIsInitialized; ///getVal() > muTest->getVal() ) { - oocoutI(nullptr,Eval) << "Using one-sided qmu - setting qmu to zero muHat = " << muHat->getVal() + if (fSigned) { + oocoutI(nullptr, Eval) << "Using signed one-sided qmu - flipping the sign of the test statistic muHat = " + << muHat->getVal() << " muTest = " << muTest->getVal() << std::endl; + flipSign = true; + } else { + oocoutI(nullptr, Eval) << "Using one-sided qmu - setting qmu to zero muHat = " << muHat->getVal() << " muTest = " << muTest->getVal() << std::endl; - qmu = 0; + qmu = 0; + } } } if (fOneSidedDiscovery ) { if ( muHat->getVal() < muTest->getVal() ) { - oocoutI(nullptr,Eval) << "Using one-sided discovery qmu - setting qmu to zero muHat = " << muHat->getVal() + if (fSigned) { + oocoutI(nullptr, Eval) + << "Using signed one-sided discovery qmu - flipping the sign of the test statistic muHat = " + << muHat->getVal() << " muTest = " << muTest->getVal() << std::endl; + flipSign = true; + } else { + oocoutI(nullptr, Eval) << "Using one-sided discovery qmu - setting qmu to zero muHat = " << muHat->getVal() << " muTest = " << muTest->getVal() << std::endl; - qmu = 0; + qmu = 0; + } } } @@ -733,7 +747,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { // asymptotic formula for pnull (for only one POI) // From fact that qmu is a chi2 with ndf=1 + // for the signed test statistic, sqrtqmu becomes negative when the best fit + // value is beyond the tested value; the Gaussian asymptotic formulae below + // remain valid also in that case double sqrtqmu = (qmu > 0) ? std::sqrt(qmu) : 0; + if (flipSign) + sqrtqmu = -sqrtqmu; double sqrtqmu_A = (qmu_A > 0) ? std::sqrt(qmu_A) : 0; @@ -758,7 +777,10 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { } - if (useQTilde ) { + // the qtilde corrections apply when the best fit value is at the boundary + // (qmu > qmu_A); they cannot apply when the sign was flipped, because then + // the best fit value is on the other side of the tested value + if (useQTilde && !flipSign) { if (fOneSided) { // for bounded one-sided (q_mu_tilde: equations 64,65) if ( qmu > qmu_A && (qmu_A > 0 || qmu > tol) ) { // to avoid case 0/0 @@ -780,8 +802,6 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const { } } - - // create an HypoTest result but where the sampling distributions are set to zero string resultname = "HypoTestAsymptotic_result"; HypoTestResult* res = new HypoTestResult(resultname.c_str(), pnull, palt); diff --git a/roofit/roostats/test/testAsymptoticCalculator.cxx b/roofit/roostats/test/testAsymptoticCalculator.cxx index 6d97b7daa5c3b..f49d8548f125a 100644 --- a/roofit/roostats/test/testAsymptoticCalculator.cxx +++ b/roofit/roostats/test/testAsymptoticCalculator.cxx @@ -1,10 +1,19 @@ // Author: Jonas Rembser, CERN 01/2025 +#include "RooDataSet.h" #include "RooMultiVarGaussian.h" +#include "RooRealVar.h" +#include "RooWorkspace.h" #include "RooStats/AsymptoticCalculator.h" +#include "RooStats/HypoTestResult.h" +#include "RooStats/ModelConfig.h" + +#include "Math/ProbFuncMathCore.h" #include "gtest/gtest.h" +#include + // Check if asymptotic datasets for counting experiments can also be generated // from the RooMultiVarGaussian. TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian) @@ -37,3 +46,66 @@ TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian) EXPECT_EQ(dataX.getVal(), mu.getVal()); } } + +// Check the signed (uncapped) one-sided profile likelihood test statistic +// enabled with AsymptoticCalculator::SetSigned() (ROOT-8257). +TEST(AsymptoticCalculator, SignedTestStatistic) +{ + using namespace RooStats; + + RooWorkspace ws; + ws.factory("Gaussian::model(x[1.5, -5, 5], mu[1.0, -5, 5], 1.0)"); + + RooRealVar &x = *ws.var("x"); + RooRealVar &mu = *ws.var("mu"); + + // A single observation at x = 1.5, so muHat = 1.5 is beyond the tested value mu = 1. + RooDataSet data{"data", "data", {x}}; + data.add({x}); + + ModelConfig nullModel{"null_model", &ws}; + nullModel.SetPdf(*ws.pdf("model")); + nullModel.SetObservables({x}); + nullModel.SetParametersOfInterest({mu}); + mu.setVal(1.0); + nullModel.SetSnapshot({mu}); + + ModelConfig altModel{nullModel}; + altModel.SetName("alt_model"); + mu.setVal(0.0); + altModel.SetSnapshot({mu}); + + AsymptoticCalculator calc{data, altModel, nullModel}; + calc.SetOneSided(true); + + // For a single Gaussian observation, muHat = x = 1.5 with sigma(muHat) = 1, + // so qmu = (x - mu)^2 = 0.25, and on the alt-hypothesis Asimov data (x = 0) + // qmu_A = mu^2 = 1. + const double sqrtqmu = 0.5; + const double sqrtqmuA = 1.0; + + // With the capped one-sided statistic, qmu is set to zero because muHat is + // beyond the tested value. + std::unique_ptr resCapped{calc.GetHypoTest()}; + EXPECT_NEAR(resCapped->NullPValue(), ROOT::Math::normal_cdf_c(0.0), 1e-3); + EXPECT_NEAR(resCapped->AlternatePValue(), ROOT::Math::normal_cdf(sqrtqmuA), 1e-3); + + // With the signed statistic, sqrt(qmu) enters the asymptotic formulae with + // a negative sign instead. + calc.SetSigned(); + std::unique_ptr resSigned{calc.GetHypoTest()}; + EXPECT_NEAR(resSigned->NullPValue(), ROOT::Math::normal_cdf_c(-sqrtqmu), 1e-3); + EXPECT_NEAR(resSigned->AlternatePValue(), ROOT::Math::normal_cdf(sqrtqmuA + sqrtqmu), 1e-3); + + // On the "right side" (muHat below the tested value) the signed and capped + // statistics must agree: test mu = 2, where sqrt(qmu) = 0.5. + mu.setVal(2.0); + nullModel.SetSnapshot({mu}); + calc.SetNullModel(nullModel); + std::unique_ptr resRightSide{calc.GetHypoTest()}; + EXPECT_NEAR(resRightSide->NullPValue(), ROOT::Math::normal_cdf_c(sqrtqmu), 1e-3); + + calc.SetSigned(false); + std::unique_ptr resRightSideCapped{calc.GetHypoTest()}; + EXPECT_NEAR(resRightSideCapped->NullPValue(), resRightSide->NullPValue(), 1e-6); +}