I am experiencing some issues while extracting limits on a process using the AsympoticCalculator.
The code is similar to the one discussed and shown in Upper limit calculation with uncertainty using Asymptotic Calculator (with CLs).
The issue is that the extracted limits happen to be numerically unstable for certain mass ranges, in particular at the beginning and at the end of the overall scanned mass range.
I’m attaching example for the upper end of the range.
I’m showing 10 random extraction (I’m performing my extraction tests on Monte Carlo, and I sample it 10 times, randomly).
It is clear that up to 9.2 GeV everything is stable and works fine; above that mass, the extractions start to be unstable.
On the other hand, the fits are stable and they evolve smoothly, returning meaningful parameters and significance, so the problem is not in the fit itself but rather in the upper limit calculator.
I’m attaching also two pictures (with random seed 30, green line) of two adjacent fits: they are very similar to each other, as expected but the returned upper limit is very different, as seeable in the previous picture.
What I noticed is that in the second fit we enter an area in which the overall fit pdf would want to go negative (see the very last couple of bins), but of course RooFit forces it to be non negative.
I reckon that this may actually be the cause of these instabilities.
I spotted this error message during the limit calculation:
[#1] INFO:ObjectHandling -- RooWorkspace::saveSnaphot(w) replacing previous snapshot with name S+B_modelNM__snapshot
----> Doing a re-scan first
----> trying with improve
[#0] ERROR:Minization -- AsymptoticCalculator: failure in fitting for qmu or qmuA 3.59799 return a dummy result
[#0] WARNING:Eval -- HypoTestInverter - Skip invalid result for point xsec = 3.59799
where xsec is the parameter of interest of my scan.
I also tried to shrink the fit window to no more than 100 GeV, so that no empty bins are included and therefore the pdf does not want to go negative, and this indeed gets rid of all the instabilities.
Do you have any idea on how to solve this issue and make the limit extraction stable even in areas in which the pdf is force to be negative?
Hi Michael,
Our RooFit expert (@StephanH) is currently on the road. Once he is back he will reply! If you did not receive an answer within 10 days please ping him here!
Axel.
In the meanwhile, given that we haven’t been able to solve our issue, I attach a minimal working example (with instructions), so that Stephan may take a look in detail as soon as he returns ForRooFitForum.zip (366.9 KB)
I think you are right that the problem is the negative PDF (which is not a PDF, actually, as these cannot be negative). I would try to protect the fit model against going negative. Since I have to go through a lot of posts, could you write down the essence of the model, e.g. PDF = a * Gauss(signal) + b * Pol4(Bkg) or something like that, and could you also explain what exactly leads to negative bins? (Is it the signal that goes negative or the background? Is it the PDF or the scale factor?) The solution could be to constrain the parameters of the signal or background or to use a different distribution such that no negative predictions are encountered.
to be more precise: the pdf would like to go negative, in some areas, but RooFit is smart enough to forcefully set the pdf to zero in these areas. To my understanding, also the normalization should be computed properly.
The core of my model is:
PDF = A * CrystalBall(signal) + B * PolN(bkg)
where N = 2, 3, 4 depending on the mass hypothesis, and A>=0 and B>=0.
So what wants to go negative is the polynomial itself.
This happens at the edges of the distributions, where the terminal bins can be empty.
I’m attaching two example fits to show where and how this happens. You can see the PDF being forced to zero.
I reckon that this may bring numerical instabilities when it comes to extract the upper limits with the AsymptoticCalculator.
I know you’re very busy right now, but any help would be greatly appreciated, as we’ve not been able to solve this for a long time.
Cheers,
Michael
P.S.: the error that I quoted in my first message happens only when is set to zero, and for close windows, where, I guess, the wiggling operated by the AsymptoticCalculator to do its computations may make the pdf “negative”.
That’s correct, but RooFit uses a very brutal method to drive the minimiser away from the parameter regions, which lead to negative PDFs. At this point, there is no proper gradient such that the fitter doesn’t know what to do next. It basically just sees that it made a mistake, and reverses that decision. If the model would be written in a way that this doesn’t happen, it would make the job of the fitter much easier. This may very well have an impact on the limits, as uncertainties are probably only correct if you are “far” from the disallowed regions. If you make sure that the fitter doesn’t end up in those disallowed regions, everything should stabilise.
You could offset the polynomial. Check out this constructor: ROOT: RooPolynomial Class Reference
You can use it to switch on the 0th (i.e. constant component), and instead of one, you can set it to a constant 2 or 5 etc. In this way, the polynomial can be made positive across the whole range. This offset will be “normalised away” when the distribution is integrated.
You can also limit the allowed range of the components such that the highest even order is always positive or you can make sure that the highest uneven component is small, because it will eventually drive the polynomial negative on one side.
Another options is to e.g. use the Chebychev Polynomials. They might behave a bit better.
If all of those don’t work, one can introduce constaint terms in the likelihood that make sure that the polynomial doesn’t go negative.
many thanks for the answer and suggestions.
I will try these options.
First of all I would like to stress that the actual fit works perfectly fine and returns perfectly meaningful results; only the limits obtained with the AsymptoticCalculator are unstable.
I am actually already using Chebychev polynomials, the unstable results that I report, and that are reproducible using the minimal working example I attached some posts above, are indeed obtained with such functions.
I also tried adding manually a very small constant to the overall PDF, in such a way that it became:
PDF = ChebyZero(constant) * F +[ A * CrystalBall(signal) + B * ChebyN(bkg) ]
where if tried fractions F of 10^-15, 10^-5, and 10^-1, the latter just to check that the result is what I expected. I do end up having a “baseline” which makes the overall PDF non-negative, but the problem is that, apparently, RooFit complains even if just one of the components of the PDF goes (or would like to go) negative, and the Chebychev component still would like to go negative.
I will test constraining the components of the Chebychev, but I would hope for another solution, cause this will likely make the overall fit “worse” in terms of likelihood/chisquare.
I’d start modifying the offset; is there anything similar for the ChebyChev polynomials too? I cannot find anything similar in the constructors https://root.cern.ch/doc/master/classRooChebychev.html#a9914b5ac39fe62addb49511de3817c0c
Yes, that’s exactly how it was designed to behave. If you are summing PDFs, it is ensured that every component of the sum is non-negative.
No, unfortunately not. The constant polynomial is always assumed to be 1. What you can do, though is to make the coefficients of the Chebychevs small. Their value range is doesn’t exceed [-1 and 1], so if you e.g. take the first five Chebychevs, and limit the coefficients to e.g. [-0.2, 0.2], they cannot possibly go negative.
Some more ideas:
Bernstein polynomials are always positive if the Chebychevs are just not working.
When you change the mass range, maybe the fits for the asymptotic calculator are starting in a very unfavourable position. Try to set the mass range, run a manual fit to get reasonable starting values, and then run the asymptotic calculator. That might work because RooFit’s way of telling the fitter that there was a problem is to return a very low likelihood. The asymptotic calculator then compares the two likelihoods, and gets unreasonable results.
short answer: I tried basically everything you suggested (Bernstein, make the parameter interval smaller, feeding the AsymptoticCalculator with better starting value), but nothing was 100% functioning: reducing the parameter range was enough for low masses, but nothing helped for high masses.
Below you see an example in which the upper limit computed with the AsymptoticCalculator unstable (i.e. it depends heavily on a single entry difference and on the parameters’ ranges and starting values).
For the last handful of bins, you can clearly see how the overall pdf “would like” to go negative, cause it is following the trend of the descending slope.
Do you have any solution on how to solve this problem?
Best,
Michael
I have one more idea:
When a PDF goes negative, the minimiser is instructed to go away from this region by handing over an extremely large value to it. When that happens, there is no way to properly estimate a gradient etc, and that’s why you might get unreliable limits.
You can actually switch that off. The minimiser will not pull away from this region, and the PDF can continue to try to go negative (but the bins where it actually does go negative are zeroed). Given that the data are also zero, that can actually be the right thing to do. You do that by myPdf.fitTo(data, ..., ..., EvalErrorWall(false), ...);
Hi @StephanH,
This sounds very nice! But… how can I use that within the AsymptoticCalculator machinery?
The minimal lines of interest should be the following:
OK, I see the problem. That requires codes changes in the AsymptoticCalculator. I could add this functionality, and you could test with the master version of ROOT.
we tried setting up Wed’s nightlie.
Setting up the proper ROOT version was successful, but no python nor python3 are able to import ROOT.
Is there a quick way to make it work?
My script is in python, and converting that into C++ could take quite some times. If there is a smart, fast way of making python recognize ROOT it would be great.
------------------------------------------------------------------
| Welcome to ROOT 6.19/01 https://root.cern |
| (c) 1995-2019, The ROOT Team; conception: R. Brun, F. Rademakers |
| Built for linuxx8664gcc on Dec 18 2019, 04:50:00 |
| From heads/master@v6-19-01-2197-g13d4ee4 |
| Try '.help', '.demo', '.license', '.credits', '.quit'/'.q' |
------------------------------------------------------------------
Error message:
[denuccio@naf-belle12 LocalFolder]$ python3
Python 3.6.8 (default, Aug 7 2019, 17:28:10)
[GCC 4.8.5 20150623 (Red Hat 4.8.5-39)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import ROOT
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/cvmfs/sft.cern.ch/lcg/nightlies/dev3/Wed/ROOT/HEAD/x86_64-centos7-gcc8-opt/lib/ROOT.py", line 24, in <module>
import cppyy
File "/cvmfs/sft.cern.ch/lcg/nightlies/dev3/Wed/ROOT/HEAD/x86_64-centos7-gcc8-opt/lib/cppyy.py", line 61, in <module>
import libPyROOT as _backend
ImportError: libvdt.so: cannot open shared object file: No such file or directory
Otherwise I ran into problems that vdt wasn’t setup and that my TBB wasn’t in the correct version (and probably more later on).
@StephanH you think it might be possible to update the nightlies page you linked to be a bit more explicit, also on how to use a python3 version by using dev3python3 instead of dev3?
