HypoTestInverter Interval does not contain a root

[beojan]

I’m using HypoTestInverter to set a 95% CL upper limit. Sometimes after running a full scan, I get errors of the form:

Error in <ROOT::Math::BrentRootFinder>: Interval does not contain a root
HypoTestInverterResult - interpolation failed for interval [-12.5978,4436.69 ]  g(xmin,xmax) =0.0117351,0 target=0.05 return inf

Is there any way I can reliably avoid this issue (e.g. by running additional points)?

EDIT: I should add, I know the limit should lie within the interval, because I determined the interval by running a coarse scan, and expanding the ±2σ bands returned by that scan.

[beojan]

It appears the underlying issue is that CLs is set to -1 if CLb is 0.
Is it possible to exclude such points when calculating the interval?

[StephanH]

Hello,

@beojan

It appears the underlying issue is that CLs is set to -1 if CLb is 0.

This may be possible, but recently, there was a also report about a change in behaviour of the Brent root finder:

Program using GetX() working in root 5.34 but not root 6.18

Dear Colleagues, I need to invert a very complicated function (i.e. given y = f(x), I wish to find x = g(y)) and neither the inverse nor the function can be given analytically. I have written a simple test program (attached) that uses GetX() to invert an equation of the form y = exp(x). It appears to work in root version 5.34 but not version 6.18. Can anyone kindly explain what I have to do to make it work in version 6.18? shine.c (1.1 KB) Thanks in advance, R Please read tips for efficie…

@moneta knows RooStats much better than me, could you have a look?

Also:
Do you happen to know where the CLs is set to -1? Maybe it’s easy to fix.

[beojan]

Do you happen to know where the CLs is set to -1? Maybe it’s easy to fix.

Here

[StephanH]

Thanks.

I started a pull request. Hopefully, automatic skipping of such points will stabilise the root finder. It should be in the nightly builds by tomorrow, which can be found at https://root.cern/nightlies.
If you are able to test one of those, I would appreciate some feedback if it works.

[RS] Make HypoTestInverterPlot more robust against failed fits.

[beojan]

It still fails to find a root with the new nightly:

HypoTestInverterResult - interpolation failed for interval [0,4436.69 ]  g(xmin,xmax) =-nan,0 target=0.05 return inf
Info in <ROOT::Math::BrentMethods::MinimStep>: Grid search failed to find a root in the  interval
Info in <ROOT::Math::BrentMethods::MinimStep>: xmin = 0 xmax = 4436.69 npts = 100
Error in <ROOT::Math::BrentRootFinder>: Interval does not contain a root
HypoTestInverterResult - interpolation failed for interval [0,4436.69 ]  g(xmin,xmax) =-nan,0 target=0.05 return inf
Info in <ROOT::Math::BrentMethods::MinimStep>: Grid search failed to find a root in the  interval
Info in <ROOT::Math::BrentMethods::MinimStep>: xmin = 0 xmax = 4436.69 npts = 100
Error in <ROOT::Math::BrentRootFinder>: Interval does not contain a root
HypoTestInverterResult - interpolation failed for interval [0,4436.69 ]  g(xmin,xmax) =-nan,0 target=0.05 return inf

The plot I get from HypoTestInverterPlot is attached.

The (failed fit?). Skipping this point. messages you added appear when I make the plot, but not when I run result->GetInterval().

[StephanH]

It makes sense that it doesn’t find a root if CLs at mu=0 is NaN, but I don’t know why that happens.

@moneta, do you have an idea?

[moneta]

If the fit failed the value can then be a NaN. You should exclude then the point with mu=0.
I would try to scan only between values slightly larger than zero until mu ~ 50-100, looking at the plot above

Lorenzo

[beojan]

Can’t these failed fits be excluded automatically? If the failed fit occurs in the middle of the mu range, it can’t be excluded by changing the bounds. It would also be helpful not to have to tweak the bounds manually.

[moneta]

Yes, I agree, we should exclude those points, as we are doing now excluding points where CL_s = -1

[StephanH]

@beojan

Can’t these failed fits be excluded automatically?

That’s what the PR was meant to do, but the problem apparently is not only the CLs values of -1. I suspect that infinities or nans have to be removed as well, but the question is where to do that best.

Should I just blindly try to skip those values at the location where I skipped the -1s, or can you give us an example that we can run?

[beojan]

I’m not sure the PR even worked for CLs=-1. The file changes was HypoTestInverterPlot.cxx. I’m not sure, but I don’t think that would affect HypoTestInverterResult'.

I think this code (https://github.com/root-project/root/blob/master/roofit/roostats/src/HypoTestInverterResult.cxx#L193) is where this skipping should be done already, and I’m probably running into a corner case that isn’t accounted for.
EDIT 2: It also looks like the above code only looks at the observed CLs, not the expected.

EDIT: Almost forgot. The example I have is just one example. With slightly different data the specific problematic points may be different.

[StephanH]

It worked partially I think. Since you didn’t include an example, we can only guess what’s happening in your case, though. I dug a bit more through the code, and found that the inverter uses a two-step procedure in RunLimit() (I hope that’s what you are using.)

1. It tries to bisect the scan range to quickly make it to the limit. Here, any funny CLs values will throw it off.
2. If it’s not close enough to the limit, it creates a HypoTestInverterPlot (== TGraph), and uses TGraph’s interpolation functions to estimate the limit. Here, the “-1 fix” will have an effect, but non-finite numbers will still throw it off. I also fixed the non-finite numbers now, but I’ll have to look more at the bisection step.

BTW: The function you posted isn’t called anywhere. I asked one of the authors why that is the case.

[beojan]

I’m not using RunLimit(). The problem with that function is that it cares only about the observed limit, not the expected limit or the uncertainty bands on this. Even with background only Asimov data (so the observed and expected limits are roughly identical) the bands can be rather inaccurate because it hasn’t run enough points in the vicinity of the band boundaries.

I use RunOnePoint to run the mu = 0 case (this was a recent change, while attempting to fix this issue) and RunFixedScan to run a logarithmic scan between 1e-2 and 1e4. I then take the ±2σ bands this returns and run a fine grained logarithmic scan (with a fresh HypoTestInverter) between -2σ/(coarse step size) and +2σ * (coarse step size) to get my final limits.

If you’re in atlas-current-physicists, you should be able to see the code here.

[StephanH]

Here’s a second try. Still flying blindly, but I hope that the fixed scan now simply steps over all CL values that are not within [0, 1]. Please let me know if it works.

[RS] Skip NaNs and infinities in HypoTestInverter.

[beojan]

I think I’ve determined what the problem is. It seems to arise only when I get CLs=-1 at the start of the series of test points (where CLs is actually asymptotically approaching 1).

ExclusionCleanup then gets confused because the order the checks are done means CLsobsprev gets updated to reflect these -1 values, leading to all points being removed. If -1 points (or better, points < 0) are removed in ExclusionCleanup then I think that function could be run (by the user if not automatically) and would correctly remove the problematic points.

[StephanH]

Third attempt:

[RS] Update HypoTestInverterResult::ExclusionCleanup.