I have been using the rs101_limitexample.C to calculate upper limits for a simple counting experiment. I was told that for small background (<10) that the MCMCCalculator was the most appropriate for this since it should give an UL of 3 events when 0 events are observed. However, I have noticed two things:
When I use obs=0, I don’t get 3 but 3.9 . If I look at the probability distribution of signal events plotted in the canvas, I see that the MCMCCalculator looks like it is underestimating near zero.
If I vary the upper limits on either s or obs I get different results. On the level of 5% which seems to much if the MC has really converged.
Does anyone have any ideas why I am seeing the 2 features above? Any solutions?
Just an update to this post. I have solved the first point. I didn’t realize that the script was calculating a central limit (equal sized lower tails and upper tails). I used the function, mc.SetLeftSideTailFraction(0.0), and I can now get 3.0 for the UL with zero expected background.
I am still seeing the UL change with the upper bound of both obs and s. I would be VERY interested to have ANY input about this.
setting very small errors on s and b makes very narrow Gaussians which is very difficult to integrate in a large range, since the phase space available is very small. This is a known effect and I have observed it also in the BayesianCalculator which uses numerical integration.
You should put sensible limits in s and b (something like 10-20 sigma, but not 1000’s or more ). Otherwise do not put any uncertainty, which avoids integrating on them
I tried taking out the uncertainty on the means so that the pdf is just a poisson. However, I still see that the UL is changing with the upper limit of obs and s. Likewise, the problem persists when I use reasonable uncertainties, e.g. 20%. I can understand that a monte carlo method would suffer when the bounds are too large though. Do you know of some good interval calculator that would be appropriate for a simple model like this (poisson(N;s,b) where s and b are smeared by a gaussian) but might also be more robust?
I have attached the script that I modified by taking out the uncertainties in case anyone is interested.
I tried the attached rs101_limitexampleNoUnc.C script with ROOT 5.28. The upper limit from MCMC was not 3.0 but 3.19171. After adding a line which changes the random seed
(RooRandom::randomGenerator())->SetSeed(100);
I’ve got a rather different MCMC upper limit on s = 2.89693. Is that a ~10% uncertainty at 20K iterations?
Then I’ve changed the number of iterations from 20K to 1M (I believe that’s a fairly large number for such a simple model) and got the limits of 3.13082 and 3.12864 respectively for the default and random seed of 100.
Next changes: set burn-in steps from 100 to 1000 and for the proposal function add a uniform fraction of 0.1. The corresponding UL results become 3.15157 and 3.10369. It’s still fluctuating and is far from 3.0.
We’ve seen similar instability problems with MCMCCalculator for more complex models. So either there’s something wrong with MCMCCalculator at low background expectations or with the way we use it.
