Some questions about "AsymptoticLimits.cc"

[desycms107]

Dear Combine experts,

I really appreciate the scripts that you provide that make limits for so many CMS analysis. But while I read through the codes, I have some questions that I hope can be clarified so that we might reduce misunderstanding from users:

1. May I ask why we set rMin directly to 0 when we use $\widetilde{q}$? I wonder if MIGRAD will get the gradient calculation wrong if $r$ isn't allowed to get slightly toward negative values, causing biased minimisation. Could we follow the definition in eq. (16) of the Asymptotic formulae paper, re-define the NLL if $\widehat{r}$ is either less than 0 or larger than $r$?
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 161 in 62005b0

   r->setMin(qtilde_ ? 0 : -r->getMax());

2. I think it is good that rMax is set to be 3 times the uncertainty of r (rErr). So the while loop has much smaller chance to become an infinite loop (than the one mentioned in point 8). But I think we should put a warning whenever the rMax is doubled to let the user know the uncertainty of r increases fast with r, which could be due to incorrect model configuration or large systematics applying to signals.
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 225 in 62005b0

   double rMin = std::max<double>(0, r->getVal()), rMax = rMin + 3 * rErr;

3. I wonder why we do a logarithmic interpolation for the observed limit but bissection method for the expected limit as default? When a user creates an Asimov dataset using "-t -1", the user would expect the observed limit to match the 50% quantile expected limit. But one wouldn't because the limit-searching methods are not same.
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 247 in 62005b0

   double rCross = rMin + (rMax-rMin)*log(clsMax/clsTarget)/log(clsMax/clsMin);

   
   I also wonder what is the mathematical basis that the r should be lienar to the natural log of CLs. According to equations like (69) in the Asymptotic formulae paper, the relation is an inverse CDF of the normal distribution, not natural log.

4. May I ask the mathematical basis fo this 80-20 interpolation ?
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 249 in 62005b0

   limit = 0.8*rCross + 0.2*rMax;

5. Is it intentional to use the function "minimize" for unconstrained fit:
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 184 in 62005b0

   minim.minimize(verbose-2);

   
   But then we use the function "improve" for the constrained fit"
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 294 in 62005b0

   if (!minimD.improve(verbose-2) && picky_) return -999;

6. I thought qmu can never become negative according to the definition of qmu in the Asymptotic formulae paper's eq (16).
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 299 in 62005b0

   double qmu = 2*(nllD_->getVal() - minNllD_); if (qmu < 0) qmu = 0;

   
   Is this to prevent numerical errors? The setting of qmu to 0 when $\widehat{r}>r$ is done here already:
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 304 in 62005b0

   qmu = 0.;

7. According to page 20 of Minuit documentation, the error level (the value to 1 $\sigma$) affects the precision of the minimiser: estimated distance to minimum (EDM) requirement is 0.001*[tolerance = 0.1]*[error-level]. The RooFit default value is 1. But in the following snippet, it changed it to 1.92
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 403 in 62005b0

   minim.setErrorLevel(0.5*pow(ROOT::Math::normal_quantile(1-0.5*(1-cl),1.0), 2)); // the 0.5 is because qmu is -2*NLL

   
   So the precision becomes worse by this change. But in the following lines, it used this less precise fit to determine whether the Asimov dataset is well-behaved.
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 415 in 62005b0

   CombineLogger::instance().log("AsymptoticLimits.cc",__LINE__,std::string(Form("[WARNING] Best fit of asimov dataset is at %s = %f (%f times %sMax), while it should be at zero",r->GetName(), r->getVal(), r->getVal()/r->getMax(), r->GetName())),__func__);

   
   What's stranger is that the reference is 0.001 times rMax, which is rather arbitrary. I suggested we turn off this warning as it might mislead careful users that turn up the verbose level, or we don't change the error-level for higher precision and could use the ratio $\widehat{r}$ over $\sigma_{\text{fit}}(r)$ as a indicator of the fit healthiness.

8. The current default way of calculating expected limit is via the "bisection", which contains a while loop:
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 517 in 62005b0

   while (rErr > std::max(rRelAccuracy_*rCross, rAbsAccuracy_)) {

   
   This implied that user might be stuck in an infinite loop if the rMax and rMin were set too narrow. In contrast, the method "stepping" has a warning when it couldn't find the required NLL shift and won't get stuck in an infinite loop:
   

   HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

   Line 557 in 62005b0

   if (verbose > 1) CombineLogger::instance().log("AsymptoticLimits.cc",__LINE__,std::string(Form("reached hard bound at %s = %f\n", r->GetName(), rCross)),__func__);

Best regards,

Kuan
B2G Combine contact

[gpetruc]

Hi,

Much of that code is rather old, so the precise motivations may have been lost.

For (3) and (4) my guess from what I committed 13 years ago ( 22a7319 ) was that it was just an empyrical approach to get the fit converge more quickly to the right value when CLs was in the right ballpark by jumping to something closer to the extrapolated intercept instead of blindly bisecting the interval.

For (5), I believe the rationale to use improve instead of minimize was that after the unconstrained fit the values and errors should already be in a reasonable ballpark and so minuit shouldn't need to re-guess them. I'm not sure that we ever tested whether it brings and improvement or not, though.

For (6) has been there since ever, I'd guess it could be due to numerical precision or stablitiy issues that could cause the NLL to go lower than the previously-found minimum.

For (7), the use of 1.92 instead of 1 is so that when Minuit computes the uncertainties it gets the 95% instead of the 68%. Yes, that reduces the default uncertainty on the search, but that shouldn't matter too much for the later check about whether the best fit value is close to zero. The 0.001 is indeed arbitrary and could be increased a bit if there's many false positives. It is unrelated to the 0.001 used by Minuit (the one in the warning is for horizontal distances, i.e. in r, and relative, while Minuit's one is on nll i.e. vertical, and absolute). The condition uses 0.001 times rMax just because anyway in minuit rMax the scale of the relevant values of the parameter r; one could have used also the sigma (e.g. r < 0.01 * sigma(r)), but at the risk of getting bogus results when the fit converges poorly and so gives a wrong estimate of sigma (e.g. sigma = zero, which could happen since this is just the first fit from minuit before any scanning of the NLL to get precise uncertainties). If the starting ranges of the parameters are sensible then you would expect sigma(r) to be smaller than rMax but not by orders of magnitude, and so the two should give reasonably similar results (e.g. if sigma(r) is of order 0.1 times rMax).
Still, that number 0.001 was put in long ago when datacards were simpler, so if nowadays on the average datacard the asimov best fit is less accurate than the tolerance of the warning message could be changed.

For (8), back in the days I had experimented a bit with different strategies for finding the NLL crossing, but didn't find one that was universally better, i.e. with low rate of failures due to numerical issues and also converging quickly, across many different datacards. That said, from a quick reading for how it's implemented the bisection should never give an infinite loop since the interval is always halved at each iteration; if the initial range is bad and so no crossing can be found, provided that rMin < rMax, it will just converge to a wrong value close to the boundary, which it not a very good behaviour but it is not an infinite loop.

[desycms107]

Hi @gpetruc,

Thanks for your reply. For the crossing search, may I ask if you know why the default option is not minos in here?

HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

Line 469 in 62005b0

if (minosAlgo_ == "minos") {

I believe that you have optimised the current method. But I believe it would be faster to do a Newton's method to take advantage of the gradient. And then I recalled that TMinuit could do fast optimisation using Davidon–Fletcher–Powell formula; TMinuit2 can even adopt the more advanced BFGS method. With these fast converging algorithm, we should be able approach to ''r'' with the desired ''qmu'' faster than a line-search.

And I am still worried that the "while" loop on L517 could be unstoppable -- it doesn't seem to have the ability to increase rMax if the NLL is always less than the target. If the user's rMax set too low, then rCross will just keep approaching rMax. But because the NLL-level at rMax is still below the target, it will never success. L518 has a tiny adjustment for rMax. But rCross is 0.5*(rMin+rMax), so I wonder how can it ever get larger than (r->getMax())? L483 only increases rMax by two times if the original rMax is less than 90% of the up-error. But what if 2 times of that is still not enough or if the up-error was wrongly calculated by minos because rMax was set too small? The default rMax in Combine seems to be 20, but it is quite common in B2G with r-limit at 100 before they optimise the analysis. I think we need more Poka-yoke here.

Best,

Kuan

[nucleosynthesis]

Hi, I don't want to interrupt the discussion but we now have a few unit tests built into our GitHub repo and they should include some limit calculations using this method. Kuan, if you would be willing to test out a few different methods to find the crossing here, and compare the performance in accuracy/speed and even robustness against certain cases (user not setting large enough limits on r etc) then we would definitely appreciate it from the development side. Just wanted to inject that into this thread, Best, Nick

…

On Wed, 12 Feb 2025, 09:52 desycms107, ***@***.***> wrote: Hi @gpetruc <https://github.com/gpetruc>, Thanks for your reply. For the crossing search, may I ask if you know why the default option is not minos in here? https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/62005b0b00ef5d284e3a4d4732c82d443f9123d2/src/AsymptoticLimits.cc#L469 I believe that you have optimised the current method. But I believe it would be faster to do a Newton's method to take advantage of the gradient. And then I recalled that TMinuit could do fast optimisation using Davidon–Fletcher–Powell formula <https://en.wikipedia.org/wiki/Davidon%E2%80%93Fletcher%E2%80%93Powell_formula>; TMinuit2 can even adopt the more advanced BFGS method <https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm>. With these fast converging algorithm, we should be able approach to ''r'' with the desired ''qmu'' faster than a line-search. And I am still worried that the "while" loop on L517 could be unstoppable -- it doesn't seem to have the ability to increase rMax if the NLL is always less than the target. If the user's rMax set too low, then rCross will just keep approaching rMax. But because the NLL-level at rMax is still below the target, it will never success. L518 has a tiny adjustment for rMax. But rCross is 0.5*(rMin+rMax), so I wonder how can it ever get larger than (r->getMax())? L483 only increases rMax by two times if the original rMax is less than 90% of the up-error. But what if 2 times of that is still not enough or if the up-error was wrongly calculated by minos because rMax was set too small? The default rMax in Combine seems to be 20, but it is quite common in B2G with r-limit at 100 before they optimise the analysis. I think we need more Poka-yoke <https://en.wikipedia.org/wiki/Poka-yoke> here. Best, Kuan — Reply to this email directly, view it on GitHub <#1037 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMEVW2VXJJ35TQZVU6VAP32PMKXVAVCNFSM6AAAAABWANUVECVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJTGE4DMMJUGE> . You are receiving this because you are subscribed to this thread.Message ID: <cms-analysis/HiggsAnalysis-CombinedLimit/issues/1037/2653186141@ github.com>

[desycms107]

Hi Nick,

Sorry I may have just pointed out things without any detailed studies. These are mostly due to some weird Combine behaviours I heard from other analysers. But I haven't checked if my observation is correct or not.

Could you point me to the unit tests? I am curious which methods are tried and the output. I will also try to make extreme examples to check if the difference between "stepping" and "bissection".

Finally, may I ask if you know the concept behind this "new" method?

HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

Line 567 in 62005b0

} else if (minosAlgo_ == "new") {

I am not sure what "k" is doing here and what is the capital X.

HiggsAnalysis-CombinedLimit/src/AsymptoticLimits.cc

Line 581 in 62005b0

// X_approx(x) = X(x, y(x1)) - k * (x-x1)^2

Best,
Kuan

[anigamova]

Hi Kuan, I think Nick is referring to the checks implemented here for example. We do not check AsymptoticLimits explicitly, probably would be good to add it there

[nucleosynthesis]

Yep - thats what I meant, I think we could add an AsymptoticLimits run there and Kuan, you could try out some comparison of the algorithms, Also for the line you pointed to, The k is just supposed to mean a step size and X is the function being approximated. Cheers, Nick

…

On Thu, Feb 13, 2025 at 11:16 AM Aliya Nigamova ***@***.***> wrote: Hi Kuan, I think Nick is referring to the checks implemented here <https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/main/.github/workflows/cvmfs-ci.yml> for example. We do not check AsymptoticLimits explicitly, probably would be good to add it there — Reply to this email directly, view it on GitHub <#1037 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMEVW4P4HRWBQJXSFS3YRT2PR5KBAVCNFSM6AAAAABWANUVECVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJWGI4DCNJVGE> . You are receiving this because you commented.Message ID: ***@***.*** com>