[RF][RS] Add SetSigned() option to the AsymptoticCalculator

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;        ///<! flag to check if calculator is initialized
       mutable int fUseQTilde;             ///< flag to indicate if using qtilde or not (-1 (default based on RooRealVar)), 0 false, 1 (true)

roofit/roostats/src/AsymptoticCalculator.cxx

@@ -701,20 +701,34 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
       useQTilde = fUseQTilde;
    }
 
-
-   //check for one side condition (remember this is valid only for one poi)
+   // check for one side condition (remember this is valid only for one poi)
+   // for a signed (uncapped) test statistic the sign of sqrt(qmu) is flipped instead of setting qmu to zero
+   bool flipSign = false;
    if (fOneSided ) {
       if ( muHat->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);

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 <memory>
+
 // 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<HypoTestResult> 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<HypoTestResult> 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<HypoTestResult> resRightSide{calc.GetHypoTest()};
+   EXPECT_NEAR(resRightSide->NullPValue(), ROOT::Math::normal_cdf_c(sqrtqmu), 1e-3);
+
+   calc.SetSigned(false);
+   std::unique_ptr<HypoTestResult> resRightSideCapped{calc.GetHypoTest()};
+   EXPECT_NEAR(resRightSideCapped->NullPValue(), resRightSide->NullPValue(), 1e-6);
+}