Revert "[RF] Replace asymptotic Poisson interval with exact Garwood intervals"

This reverts commit 770f648.

Author: guitargeek

README/ReleaseNotes/v640/index.md

@@ -205,28 +205,6 @@ As part of this migration, the following build options are deprecated. From ROOT
   overhead in `RooAbsRealLValue::inRange()`, which is visible in many-parameter
   fits. Therefore, these functions are removed.
 
-### New implementation of `RooHistError::getPoissonInterval`
-
-**RooHistError::getPoissonInterval** was reimplemented to use an exact chi-square–based construction (Garwood interval) instead of asymptotic approximations and lookup tables.
-
-Previously:
-
-  * For `n > 100`, a hard-coded asymptotic formula was used, which was not explained in the documentation.
-  * That approximation corresponds to inverting the Poisson score test and was only correct for `nSigma = 1`.
-  * The behavior for `n > 100` was not statistically consistent because of the hard transition between exact formula and approximation, resulting in a discrete and unexpected jump in approximation bias.
-  * For small `n`, numerical root finding and a lookup table were used.
-
-Now:
-
-  * The interval is computed directly using chi-square quantiles for all `n`.
-  * The construction provides exact frequentist coverage.
-  * The implementation works consistently for arbitrary `nSigma`.
-  * The hard-coded `n > 100` threshold, lookup table, and numerical root finding were removed.
-  * For `n = 0`, a one-sided upper limit is used (with lower bound fixed at 0), consistent with the physical constraint `mu ≥ 0`.
-
-Results may differ from previous ROOT versions for `n > 100` or `nSigma != 1`.
-The new implementation is statistically consistent and recommended.
-
 ## RDataFrame
 
 - The message shown in ROOT 6.38 to inform users about change of default compression setting used by Snapshot (was 101 before 6.38, became 505 in 6.38) is now removed.

roofit/roofitcore/inc/RooHistError.h

@@ -37,6 +37,13 @@ class RooHistError {
   static RooAbsFunc *createBinomialSum(Int_t n, Int_t m, bool eff) ;
 
 private:
+
+
+  bool getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma= 1) const;
+  double _poissonLoLUT[1000] ;
+  double _poissonHiLUT[1000] ;
+
+  RooHistError();
   double seek(const RooAbsFunc &f, double startAt, double step, double value) const;
 
   // -----------------------------------------------------------

roofit/roofitcore/src/RooHistError.cxx

@@ -27,13 +27,10 @@ a specified area of a Poisson or Binomail error distribution.
 #include "RooHistError.h"
 #include "RooBrentRootFinder.h"
 #include "RooMsgService.h"
-
-#include "Math/QuantFuncMathCore.h" // ROOT::Math::chisquared_quantile, chisquared_quantile_c
+#include "TMath.h"
 
 #include "Riostream.h"
-
 #include <cassert>
-#include <cmath>
 
 using std::endl;
 
@@ -50,48 +47,76 @@ const RooHistError &RooHistError::instance()
 }
 
 
+////////////////////////////////////////////////////////////////////////////////
+/// Construct our singleton object.
+
+RooHistError::RooHistError()
+{
+  // Initialize lookup table ;
+  Int_t i ;
+  for (i=0 ; i<1000 ; i++) {
+    getPoissonIntervalCalc(i,_poissonLoLUT[i],_poissonHiLUT[i],1.) ;
+  }
+
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+/// Return a confidence interval for the expected number of events given n
+/// observed (unweighted) events. The interval will contain the same probability
+/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
+/// the point estimate n (ie, for small n) in which case the interval is adjusted
+/// to start at n. This method uses a lookup table to return precalculated results
+/// for n<1000
+
+bool RooHistError::getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma) const
+{
+  // Use lookup table for most common cases
+  if (n<1000 && nSigma==1.) {
+    mu1=_poissonLoLUT[n] ;
+    mu2=_poissonHiLUT[n] ;
+    return true ;
+  }
+
+  // Forward to calculation method
+  bool ret =  getPoissonIntervalCalc(n,mu1,mu2,nSigma) ;
+  return ret ;
+}
+
+
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Calculate a confidence interval for the expected number of events given n
 /// observed (unweighted) events. The interval will contain the same probability
 /// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
 /// the point estimate n (ie, for small n) in which case the interval is adjusted
 /// to start at n.
 
-bool RooHistError::getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma) const
+bool RooHistError::getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma) const
 {
-   // sanity checks
-   if (n < 0) {
-      oocoutE(nullptr, Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n
-                                 << std::endl;
-      return false;
-   }
-   if (!(nSigma > 0.0)) {
-      oocoutE(nullptr, Plotting) << "RooHistError::getPoissonInterval: nSigma must be > 0, got " << nSigma << std::endl;
-      return false;
-   }
-
-   // Convert "number of sigmas" to central Gaussian probability beta, and
-   // corresponding two-sided tail probability alpha.
-   const double beta = std::erf(nSigma / std::sqrt(2.0));
-   const double alpha = 1.0 - beta;
-
-   // Special case n = 0 (boundary at mu >= 0).
-   // Use a one-sided interval including the MLE mu=0:
-   //   mu1 = 0
-   //   P(N <= 0 | mu2) = alpha,
-   // with alpha = 1 - erf(nSigma / sqrt(2)).
-   if (n == 0) {
-      mu1 = 0.0;
-      mu2 = 0.5 * ROOT::Math::chisquared_quantile(1.0 - alpha, 2.0 * (n + 1));
-      return true;
-   }
-
-   // For n>0, use the central (equal-tailed) Garwood interval, which
-   // corresponds to allocating alpha/2 in each tail.
-   const double a2 = 0.5 * alpha;
-   mu1 = 0.5 * ROOT::Math::chisquared_quantile(a2, 2.0 * n);
-   mu2 = 0.5 * ROOT::Math::chisquared_quantile_c(a2, 2.0 * (n + 1));
-   return true;
+  // sanity checks
+  if(n < 0) {
+    oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << std::endl;
+    return false;
+  }
+
+  // use asymptotic error if possible
+  if(n > 100) {
+    mu1= n - sqrt(n+0.25) + 0.5;
+    mu2= n + sqrt(n+0.25) + 0.5;
+    return true;
+  }
+
+  // create a function object to use
+  PoissonSum upper(n);
+  if(n > 0) {
+    PoissonSum lower(n-1);
+    return getInterval(&upper,&lower,(double)n,1.0,mu1,mu2,nSigma);
+  }
+
+  // Backup solution for negative numbers
+  return getInterval(&upper,nullptr,(double)n,1.0,mu1,mu2,nSigma);
 }