From 5d26660fb3aeaef1fc5e5efdb6383ae7602c9f0d Mon Sep 17 00:00:00 2001 From: "y.maslov" Date: Sun, 6 Sep 2026 23:07:06 +0700 Subject: [PATCH] [math] Remove TRolke model 1 Numerical solution implemented in model 1 contains errors (wrong signs in lines 1068-1069 of LikeGradMod1) that may lead to wrong limits estimation. These errors can't be fixed without massive changes in algorithm. Detailed report about errors can be found here: https://github.com/y-maslov/TRolke-Model-1-Implementation-Errors AI disclosure: ChatGPT was used to assist with translating the technical report linked above from Russian into English. The technical analysis, identified issues, conclusions, and code changes were produced and reviewed by the authors. --- math/physics/inc/TRolke.h | 9 +-- math/physics/src/TRolke.cxx | 140 ++-------------------------------- tutorials/legacy/math/Rolke.C | 54 ------------- 3 files changed, 7 insertions(+), 196 deletions(-) diff --git a/math/physics/inc/TRolke.h b/math/physics/inc/TRolke.h index 801d68e6f34bf..78ce61b24cfb4 100644 --- a/math/physics/inc/TRolke.h +++ b/math/physics/inc/TRolke.h @@ -65,11 +65,7 @@ class TRolke : public TObject // LIKELIHOOD ROUTINE Double_t Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what); - //MODEL 1 - Double_t EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what); - Double_t LikeMod1(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); - void ProfLikeMod1(Double_t mu, Double_t &b, Double_t &e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); - Double_t LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); + //MODEL 1 is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm //MODEL 2 Double_t EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t sde, Double_t tau, Int_t what); @@ -131,8 +127,7 @@ class TRolke : public TObject } // The Set methods for the different models are described in Rolke.cxx - // model 1 - void SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); + // model 1 is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm // model 2 void SetPoissonBkgGaussEff(Int_t x, Int_t y, Double_t em, Double_t tau, Double_t sde); diff --git a/math/physics/src/TRolke.cxx b/math/physics/src/TRolke.cxx index 39535830ed29e..b673e39d4b1e6 100644 --- a/math/physics/src/TRolke.cxx +++ b/math/physics/src/TRolke.cxx @@ -25,6 +25,7 @@ If unsure, first consider model 3, 4 or 5. 1: SetPoissonBkgBinomEff(x,y,z,tau,m) +\warning This model is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm ~~~ Background: Poisson Efficiency: Binomial @@ -190,30 +191,6 @@ TRolke::~TRolke() { } -//////////////////////////////////////////////////////////////////////////////// -/// Model 1: Background - Poisson, Efficiency - Binomial -/// - x : number of observed events in the experiment -/// - y : number of observed events in background region -/// - z : number of MC events observed -/// - tau : ratio parameter (read TRolke.cxx for details) -/// - m : number of MC events generated - -void TRolke::SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - SetModelParameters( - x , // Int_t x, - y , // Int_t y, - z , // Int_t z, - 0 , // Double_t bm, - 0 , // Double_t em, - 0 , // Double_t e, - 1 , // Int_t mid, - 0 , // Double_t sde, - 0 , // Double_t sdb, - tau, // Double_t tau, - 0 , // Double_t b, - m); // Int_t m -} //////////////////////////////////////////////////////////////////////////////// /// Model 2: Background - Poisson, Efficiency - Gaussian @@ -934,7 +911,10 @@ Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, { switch (mid) { case 1: - return EvalLikeMod1(mu, x, y, z, tau, m, what); + std::cerr << "TRolke is removed due to an implementation error " << + "that may lead to wrong limits estimation and can't be " << + "fixed without massive changes in the algorithm" << std::endl; + return 0; case 2: return EvalLikeMod2(mu, x, y, em, sde, tau, what); case 3: @@ -956,120 +936,10 @@ Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, return 0; } -//////////////////////////////////////////////////////////////////////////////// -/// Calculates the Profile Likelihood for MODEL 1: -/// Poisson background/ Binomial Efficiency -/// - what = 1: Maximum likelihood estimate is returned -/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned. -/// - what = 3: Profile Likelihood of Test hypothesis is returned -/// otherwise parameters as described in the beginning of the class) - -Double_t TRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what) -{ - Double_t f = 0; - Double_t zm = Double_t(z) / m; - - if (what == 1) { - f = (x - y / tau) / zm; - } - - if (what == 2) { - mu = (x - y / tau) / zm; - Double_t b = y / tau; - Double_t e = zm; - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } - - if (what == 3) { - if (mu == 0) { - Double_t b = (x + y) / (1.0 + tau); - Double_t e = zm; - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } else { - Double_t e = 0; - Double_t b = 0; - ProfLikeMod1(mu, b, e, x, y, z, tau, m); - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } - } - - return f; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Profile Likelihood function for MODEL 1: -/// Poisson background/ Binomial Efficiency - -Double_t TRolke::LikeMod1(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - double s = e*mu+b; - double lls = - s; - if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x); - double bg = tau*b; - double llb = -bg; - if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y); - - double lle = 0; // binomial log-like - if (z == 0) lle = m * TMath::Log(1-e); - else if ( z == m) lle = m * TMath::Log(e); - else lle = z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z); - - double f = 2*( lls + llb + lle); - return f; -} - - // this code is non-sense - // need to solve using Minuit struct LikeFunction1 { }; -//////////////////////////////////////////////////////////////////////////////// -/// Helper for calculation of estimates of efficiency and background for model 1 - -void TRolke::ProfLikeMod1(Double_t mu, Double_t &b, Double_t &e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - Double_t med = 0.0, fmid; - Int_t maxiter = 1000; - Double_t acc = 0.00001; - Double_t emin = ((m + mu * tau) - TMath::Sqrt((m + mu * tau) * (m + mu * tau) - 4 * mu * tau * z)) / 2 / mu / tau; - - Double_t low = TMath::Max(1e-10, emin + 1e-10); - Double_t high = 1 - 1e-10; - - for (Int_t i = 0; i < maxiter; i++) { - med = (low + high) / 2.; - - fmid = LikeGradMod1(med, mu, x, y, z, tau, m); - - if (high < 0.5) acc = 0.00001 * high; - else acc = 0.00001 * (1 - high); - - if ((high - low) < acc*high) break; - - if (fmid > 0) low = med; - else high = med; - } - - e = med; - Double_t eta = Double_t(z) / e - Double_t(m - z) / (1 - e); - - b = Double_t(y) / (tau - eta / mu); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Gradient model likelihood - -Double_t TRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - Double_t eta, etaprime, bprime, f; - eta = static_cast(z) / e - static_cast(m - z) / (1.0 - e); - etaprime = (-1) * (static_cast(m - z) / ((1.0 - e) * (1.0 - e)) + static_cast(z) / (e * e)); - Double_t b = y / (tau - eta / mu); - bprime = (b * b * etaprime) / mu / y; - f = (mu + bprime) * (x / (e * mu + b) - 1) + (y / b - tau) * bprime + eta; - return f; -} - //////////////////////////////////////////////////////////////////////////////// /// Calculates the Profile Likelihood for MODEL 2: /// Poisson background/ Gauss Efficiency diff --git a/tutorials/legacy/math/Rolke.C b/tutorials/legacy/math/Rolke.C index eb78e0a322d06..2a42ac11bcf79 100644 --- a/tutorials/legacy/math/Rolke.C +++ b/tutorials/legacy/math/Rolke.C @@ -41,31 +41,6 @@ void Rolke() Double_t ul ; // upper limit Double_t ll ; // lower limit -//----------------------------------------------- -// Model 1 assumes: -// -// Poisson uncertainty in the background estimate -// Binomial uncertainty in the efficiency estimate -// - cout << endl<<" ======================================================== " <