[RF][RS] Unify BayesianCalculator::GetInterval() 1D and ND code paths

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:

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::Integrator>(ROOT::Math::IntegratorOneDim::GetType(integType));
+         fIntegratorOneDim->SetFunction(fLikelihood);
+      } else {
+         fIntegrator = std::make_unique<ROOT::Math::IntegratorMultiDim>(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>((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::Integrator>(ROOT::Math::IntegratorOneDim::GetType(rhs.fIntegratorOneDim->Name().c_str()));
+         fIntegratorOneDim->SetFunction(fLikelihood);
+      } else {
+         fIntegrator = std::make_unique<ROOT::Math::IntegratorMultiDim>(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<RooFunctor>(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<RooFunctor> 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<ROOT::Math::IntegratorMultiDim> fIntegrator; // multi-dim integrator (for >=2 dims)
+   mutable std::unique_ptr<ROOT::Math::Integrator> fIntegratorOneDim;   // 1D integrator (for POI-only case)
    mutable std::vector<double> fXmin;    // min value of parameters (poi+nuis) -
    mutable std::vector<double> 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<RooRealVar*>(fPOI.first() );
-   assert(poi);
-
-   fValidInterval = false;
-   if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf();
-   if (!fPosteriorPdf) return;
-
-   std::unique_ptr<RooAbsReal> cdf{fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf())};
-   if (!cdf) return;
-
-   std::unique_ptr<RooAbsFunc> 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
 

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)

roofit/roostats/test/testBayesianCalculator.cxx

@@ -0,0 +1,65 @@
+#include <RooAddPdf.h>
+#include <RooArgSet.h>
+#include <RooDataSet.h>
+#include <RooExponential.h>
+#include <RooGaussian.h>
+#include <RooRandom.h>
+#include <RooRealVar.h>
+#include <RooStats/BayesianCalculator.h>
+#include <RooStats/SimpleInterval.h>
+#include <RooUniform.h>
+
+#include <gtest/gtest.h>
+
+/// 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<RooDataSet> 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<RooStats::SimpleInterval> 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<RooStats::SimpleInterval> 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);
+}