Notebook

rf614_binned_fit_problems¶

A tutorial that explains you how to solve problems with binning effects and numerical stability in binned fits.

Introduction¶

In this tutorial, you will learn three new things:

How to reduce the bias in binned fits by changing the definition of the normalization integral
How to completely get rid of binning effects by integrating the pdf over each bin
How to improve the numeric stability of fits with a greatly different number of events per bin, using a constant per-bin counterterm

Author: Jonas Rembser
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Thursday, March 23, 2023 at 10:53 AM.

Generate binned Asimov dataset for a continuous pdf. One should in principle be able to use pdf.generateBinned(x, nEvents, RooFit::ExpectedData()). Unfortunately it has a problem: it also has the bin bias that this tutorial demostrates, to if we would use it, the biases would cancel out.

In [1]:

%%cpp -d
std::unique_ptr<RooDataHist> generateBinnedAsimov(RooAbsPdf const &pdf, RooRealVar &x, int nEvents)
{
   auto dataH = std::make_unique<RooDataHist>("dataH", "dataH", RooArgSet{x});
   RooAbsBinning &xBinning = x.getBinning();
   for (int iBin = 0; iBin < x.numBins(); ++iBin) {
      x.setRange("bin", xBinning.binLow(iBin), xBinning.binHigh(iBin));
      std::unique_ptr<RooAbsReal> integ{pdf.createIntegral(x, RooFit::NormSet(x), RooFit::Range("bin"))};
      integ->getVal();
      dataH->set(iBin, nEvents * integ->getVal(), -1);
   }
   return dataH;
}

Force numeric integration and do this numeric integration with the RooBinIntegrator, which sums the function values at the bin centers.

In [2]:

%%cpp -d
void enableBinIntegrator(RooAbsReal &func, int numBins)
{
   RooNumIntConfig customConfig(*func.getIntegratorConfig());
   customConfig.method1D().setLabel("RooBinIntegrator");
   customConfig.getConfigSection("RooBinIntegrator").setRealValue("numBins", numBins);
   func.setIntegratorConfig(customConfig);
   func.forceNumInt(true);
}

Reset the integrator config to disable the RooBinIntegrator.

In [3]:

%%cpp -d
void disableBinIntegrator(RooAbsReal &func)
{
   func.setIntegratorConfig();
   func.forceNumInt(false);
}

In [4]:

using namespace RooFit;

Silence info output for this tutorial

In [5]:

RooMsgService::instance().getStream(1).removeTopic(Minimization);
RooMsgService::instance().getStream(1).removeTopic(Fitting);
RooMsgService::instance().getStream(1).removeTopic(Generation);

Exponential example¶

Set up the observable

In [6]:

RooRealVar x{"x", "x", 0.1, 5.1};
x.setBins(10); // fewer bins so we have larger binning effects for this demo

Let's first look at the example of an exponential function

In [7]:

RooRealVar c{"c", "c", -1.8, -5, 5};
RooExponential expo{"expo", "expo", x, c};

Generate an Asimov dataset such that the only difference between the fit result and the true parameters comes from binning effects.

In [8]:

std::unique_ptr<RooAbsData> expoData{generateBinnedAsimov(expo, x, 10000)};

[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'bin' created with bounds [0.1,0.6]

If you do the fit the usual was in RooFit, you will get a bias in the result. This is because the continuous, normalized pdf is evaluated only at the bin centers.

In [9]:

std::unique_ptr<RooFitResult> fit1{expo.fitTo(*expoData, Save(), PrintLevel(-1), SumW2Error(false))};
fit1->Print();

  RooFitResult: minimized FCN value: 4754.37, estimated distance to minimum: 6.17653e-09
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     c   -1.6862e+00 +/-  1.70e-02

In the case of an exponential function, the bias that you get by evaluating the pdf only at the bin centers is a constant scale factor in each bin. Here, we can do a trick to get rid of the bias: we also evaluate the normalization integral for the pdf the same way, i.e., summing the values of the unnormalized pdf at the bin centers. Like this the bias cancels out. You can achieve this by customizing the way how the pdf is integrated (see also the rf901_numintconfig tutorial).

In [10]:

enableBinIntegrator(expo, x.numBins());
std::unique_ptr<RooFitResult> fit2{expo.fitTo(*expoData, Save(), PrintLevel(-1), SumW2Error(false))};
fit2->Print();
disableBinIntegrator(expo);

[#0] WARNING:Integration -- RooBinIntegrator::RooBinIntegrator WARNING: integrand provide no binning definition observable #0 substituting default binning of 10 bins
[#1] INFO:NumericIntegration -- RooRealIntegral::init(expo_Int[x]) using numeric integrator RooBinIntegrator to calculate Int(x)

  RooFitResult: minimized FCN value: 4440.6, estimated distance to minimum: 1.1198e-06
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     c   -1.8000e+00 +/-  1.87e-02

Power law example¶

Let's not look at another example: a power law $$x^a$$.

In [11]:

RooRealVar a{"a", "a", -0.3, -5.0, 5.0};
RooPower powerlaw{"powerlaw", "powerlaw", x, RooConst(1.0), a};
std::unique_ptr<RooAbsData> powerlawData{generateBinnedAsimov(powerlaw, x, 10000)};

Again, if you do a vanilla fit, you'll get a bias

In [12]:

std::unique_ptr<RooFitResult> fit3{powerlaw.fitTo(*powerlawData, Save(), PrintLevel(-1), SumW2Error(false))};
fit3->Print();

  RooFitResult: minimized FCN value: 15816.4, estimated distance to minimum: 9.94074e-07
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     a   -2.6106e-01 +/-  1.06e-02

This time, the bias is not the same factor in each bin! This means our trick by sampling the integral in the same way doesn't cancel out the bias completely. The average bias is canceled, but there are per-bin biases that remain. Still, this method has some value: it is cheaper than rigurously correcting the bias by integrating the pdf in each bin. So if you know your per-bin bias variations are small or performance is an issue, this approach can be sufficient.

In [13]:

enableBinIntegrator(powerlaw, x.numBins());
std::unique_ptr<RooFitResult> fit4{powerlaw.fitTo(*powerlawData, Save(), PrintLevel(-1), SumW2Error(false))};
fit4->Print();
disableBinIntegrator(powerlaw);

[#0] WARNING:Integration -- RooBinIntegrator::RooBinIntegrator WARNING: integrand provide no binning definition observable #0 substituting default binning of 10 bins
[#1] INFO:NumericIntegration -- RooRealIntegral::init(powerlaw_Int[x]) using numeric integrator RooBinIntegrator to calculate Int(x)

  RooFitResult: minimized FCN value: 15739.9, estimated distance to minimum: 9.98948e-07
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     a   -3.1481e-01 +/-  1.15e-02

To get rid of the binning effects in the general case, one can use the IntegrateBins() command argument. Now, the pdf is not evaluated at the bin centers, but numerically integrated over each bin and divided by the bin width. The parameter for IntegrateBins() is the required precision for the numeric integrals. This is computationally expensive, but the bias is now not a problem anymore.

In [14]:

std::unique_ptr<RooFitResult> fit5{
   powerlaw.fitTo(*powerlawData, IntegrateBins(1e-3), Save(), PrintLevel(-1), SumW2Error(false))};
fit5->Print();

  RooFitResult: minimized FCN value: 15739.6, estimated distance to minimum: 2.09956e-09
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     a   -3.0000e-01 +/-  1.07e-02

Improving numerical stability¶

There is one more problem with binned fits that is related to the binning effects because often, a binned fit is affected by both problems.

The issue is numerical stability for fits with a greatly different number of events in each bin. For each bin, you have a term $$n\log(p)$$ in the NLL, where $$n$$ is the number of observations in the bin, and $$p$$ the predicted probability to have an event in that bin. The difference in the logarithms for each bin is small, but the difference in $$n$$ can be orders of magnitudes! Therefore, when summing these terms, lots of numerical precision is lost for the bins with less events.

We can study this with the example of an exponential plus a Gaussian. The Gaussian is only a faint signal in the tail of the exponential where there are not so many events. And we can't afford any precision loss for these bins, otherwise we can't fit the Gaussian.

In [15]:

x.setBins(100); // It's not about binning effects anymore, so reset the number of bins.

RooRealVar mu{"mu", "mu", 3.0, 0.1, 5.1};
RooRealVar sigma{"sigma", "sigma", 0.5, 0.01, 5.0};
RooGaussian gauss{"gauss", "gauss", x, mu, sigma};

RooRealVar nsig{"nsig", "nsig", 10000, 0, 1e9};
RooRealVar nbkg{"nbkg", "nbkg", 10000000, 0, 1e9};
RooRealVar frac{"frac", "frac", nsig.getVal() / (nsig.getVal() + nbkg.getVal()), 0.0, 1.0};

RooAddPdf model{"model", "model", {gauss, expo}, {nsig, nbkg}};

std::unique_ptr<RooAbsData> modelData{model.generateBinned(x)};

[#0] PROGRESS:Generation -- RooAbsPdf::generateBinned(model) Performing costly accept/reject sampling. If this takes too long, use extended mode to speed up the process.

Set the starting values for the Gaussian parameters away from the true value such that the fit is not trivial.

In [16]:

mu.setVal(2.0);
sigma.setVal(1.0);

std::unique_ptr<RooFitResult> fit6{model.fitTo(*modelData, Save(), PrintLevel(-1), SumW2Error(false))};
fit6->Print();

  RooFitResult: minimized FCN value: -1.47174e+08, estimated distance to minimum: 0.324117
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=-1 HESSE=4 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     c   -1.7972e+00 +/-  7.39e-04
                    mu    2.9756e+00 +/-  3.90e-02
                  nbkg    1.0001e+07 +/-  3.25e+03
                  nsig    9.4264e+03 +/-  7.36e+02
                 sigma    4.6849e-01 +/-  2.75e-02

You should see in the previous fit result that the fit did not converge: the MINIMIZE return code should be -1 (a successful fit has status code zero).

To improve the situation, we can apply a numeric trick: if we subtract in each bin a constant counterterm $$n\log(n/N)$$, we get terms for each bin that are closer to each other in order of magnitude as long as the initial model is not extremely off. Proving this mathematically is left as an excercise to the reader.

This counterterms can be enabled in RooFit if you use a binned RooDataHist to do your fit and pass the Offset("bin") option to RooAbsPdf::fitTo() or RooAbsPdf::createNLL().

In [17]:

std::unique_ptr<RooFitResult> fit7{
   model.fitTo(*modelData, Offset("bin"), Save(), PrintLevel(-1), SumW2Error(false))};
fit7->Print();

  RooFitResult: minimized FCN value: -0.151899, estimated distance to minimum: 3.26001e-06
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     c   -1.7971e+00 +/-  7.26e-04
                    mu    2.9942e+00 +/-  3.64e-02
                  nbkg    1.0001e+07 +/-  3.24e+03
                  nsig    9.2412e+03 +/-  6.94e+02
                 sigma    4.5770e-01 +/-  2.59e-02

You should now see in the last fit result that the fit has converged.

Draw all canvases

In [18]:

%jsroot on
gROOT->GetListOfCanvases()->Draw()