Dear experts,
In a couple of ATLAS analyses we are using a particular setup based RooStats and HistFactory to compute significances of potential excesses. In one of the analyses, we have for a while been having numerical problems with failing fits. Some investigation and comparison to an independent code led to the suspicion that the problems were caused by very large log-likelihood values (order 10^7) associated with bins with large contents (we have bin contents under the background hypothesis ranging from order 10^5-10^6 down to essentially 0). This suspicion was further strengthened when the other analysis also started to include bins with larger content and saw the same kind of behavior.
While we do have a working alternative in the before mentioned independent code, there is still interest to keep the HistFactory setup going, so we therefore kindly request your assistance to get to the bottom of this and hopefully solve the problem. Being not an expert in RooStats, the information below may be incomplete, so please request any additional information as needed.
---- The code setup -----
The setup is based on HistFactory, and the central ingredient in the setting up of the workspace seems to be the lines
#include <RooStats/HistFactory/MakeModelAndMeasurementsFast.h>
RooWorkspace* ws = MakeModelAndMeasurementFast(meas);Here, the measurement “meas” is created along the lines of
Measurement meas(name.data(), name.data());
meas.SetPOI("mu");
meas.SetLumi(lumi);
Channel channel_ee("channel_ee");
channel_ee.SetStatErrorConfig(relStatErrorThreshold, "Poisson");
Sample sample_ee_bkg_dy("sample_ee_bkg_dy", "h_ee_bkg_dy", inputFile);
HistoSys h_elecSys(sysString);
sample_ee_bkg_dy.AddHistoSys(h_elecSys);
channel_ee.AddSample(sample_ee_bkg_dy);
channel_ee.SetData("h_ee_data", inputFile);
meas.AddChannel(channel_ee);----- Observations ------
Comparisons to the independent code clearly show that the log-likelihood resulting from this setup is not the obvious sum over Poisson log-likelihood values for each bin (which is used in the independent code). The order of magnitude suggests that the Poisson factorial normalization is missing, but further observations (see next paragraph) indicate that it may not be as simple as this.
Having received the tip of looking into the RooNLLVar::evaluatePartition function, we printed some debugging information, and found the following surprising structure:
-
The boolean “_binnedPdf” is false, even though the inputs we provide are all in histogram form.
-
The code appears to “pretend” to be doing an unbinned fit where it is looping over the “events” (which appear to be the bins in our histogram), adding some “logPdf” value to the log-likelihood which is multiplied by the “event weight” (which appears to be the bin contents in data). We did not manage to fully understand what is the interpretation of “logPdf” in this context.
-
Finally, a Poisson term is added in the log-likelihood representing the total normalization of the whole histogram.
It is not clear to us at present whether this calculation is meant to be equivalent to the obvious sum over Poisson log-likelihood values for all bins, or whether the statistical model/interpretation of the data is somehow different in this setup. (In the cases where the fits converge in RooStats, the results seem to be identical to those of the independent code using simple Poisson probabilities).
Furthermore, we would be interested to know how the numerical stability can be improved. It is indeed the case that the factorial normalization is missing from the Poisson term 3), but adding it back in does not help, as the numbers from 2) are already large when bin contents are large.
Best regards,
Magnar K. Bugge