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); +}