ROOT   Reference Guide
Searching...
No Matches
rf614_binned_fit_problems.C File Reference

## Detailed Description

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:

1. How to reduce the bias in binned fits by changing the definition of the normalization integral
2. How to completely get rid of binning effects by integrating the pdf over each bin
3. How to improve the numeric stability of fits with a greatly different number of events per bin, using a constant per-bin counterterm
// 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
// demonstrates, to if we would use it, the biases would cancel out.
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.
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.
void disableBinIntegrator(RooAbsReal &func)
{
func.forceNumInt(false);
}
{
using namespace RooFit;
// Silence info output for this tutorial
// Exponential example
// -------------------
// Set up the observable
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
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.
std::unique_ptr<RooAbsData> expoData{generateBinnedAsimov(expo, x, 10000)};
// 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.
std::unique_ptr<RooFitResult> fit1{expo.fitTo(*expoData, Save(), PrintLevel(-1), SumW2Error(false))};
fit1->Print();
// 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
enableBinIntegrator(expo, x.numBins());
std::unique_ptr<RooFitResult> fit2{expo.fitTo(*expoData, Save(), PrintLevel(-1), SumW2Error(false))};
fit2->Print();
disableBinIntegrator(expo);
// Power law example
// -----------------
// Let's not look at another example: a power law \f[x^a\f].
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
std::unique_ptr<RooFitResult> fit3{powerlaw.fitTo(*powerlawData, Save(), PrintLevel(-1), SumW2Error(false))};
fit3->Print();
// 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.
enableBinIntegrator(powerlaw, x.numBins());
std::unique_ptr<RooFitResult> fit4{powerlaw.fitTo(*powerlawData, Save(), PrintLevel(-1), SumW2Error(false))};
fit4->Print();
disableBinIntegrator(powerlaw);
// 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.
std::unique_ptr<RooFitResult> fit5{
powerlaw.fitTo(*powerlawData, IntegrateBins(1e-3), Save(), PrintLevel(-1), SumW2Error(false))};
fit5->Print();
// 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 \f[n\log(p)\f] in
// the NLL, where \f[n\f] is the number of observations in the bin, and
// \f[p\f] the predicted probability to have an event in that bin. The
// difference in the logarithms for each bin is small, but the difference in
// \f[n\f] 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.
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)};
// Set the starting values for the Gaussian parameters away from the true
// value such that the fit is not trivial.
mu.setVal(2.0);
sigma.setVal(1.0);
std::unique_ptr<RooFitResult> fit6{model.fitTo(*modelData, Save(), PrintLevel(-1), SumW2Error(false))};
fit6->Print();
// 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 \f[n\log(n/N)\f], 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 exercise to the reader.
// This counterterms can be enabled by passing the Offset("bin") option to
// RooAbsPdf::fitTo() or RooAbsPdf::createNLL().
std::unique_ptr<RooFitResult> fit7{
model.fitTo(*modelData, Offset("bin"), Save(), PrintLevel(-1), SumW2Error(false))};
fit7->Print();
// You should now see in the last fit result that the fit has converged.
}
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:320
Abstract base class for RooRealVar binning definitions.
virtual void setRange(double xlo, double xhi)=0
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual RooFit::OwningPtr< RooDataHist > generateBinned(const RooArgSet &whatVars, double nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}) const
As RooAbsPdf::generateBinned(const RooArgSet&, const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,...
Definition RooAbsPdf.h:110
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
void setIntegratorConfig()
Remove the specialized numeric integration configuration associated with this object.
virtual void forceNumInt(bool flag=true)
Definition RooAbsReal.h:169
const RooNumIntConfig * getIntegratorConfig() const
Return the numeric integration configuration used for this object.
Efficient implementation of a sum of PDFs of the form.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Exponential PDF.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
static RooMsgService & instance()
Return reference to singleton instance.
StreamConfig & getStream(Int_t id)
Holds the configuration parameters of the various numeric integrators used by RooRealIntegral.
RooPower implements a power law PDF of the form.
Definition RooPower.h:20
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
RooConstVar & RooConst(double val)
RooCmdArg NormSet(Args_t &&... argsOrArgSet)
RooCmdArg IntegrateBins(double precision)
Integrate the PDF over bins.
RooCmdArg Offset(std::string const &mode)
RooCmdArg Save(bool flag=true)
RooCmdArg SumW2Error(bool flag)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Range(const char *rangeName, bool adjustNorm=true)
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Definition fit1.py:1
void removeTopic(RooFit::MsgTopic oldTopic)
[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'bin' created with bounds [0.1,0.6]
RooFitResult: minimized FCN value: 4754.37, estimated distance to minimum: 3.09852e-09
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
c -1.6862e+00 +/- 1.70e-02
[#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: 5.599e-07
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
c -1.8000e+00 +/- 1.87e-02
RooFitResult: minimized FCN value: 15816.4, estimated distance to minimum: 4.97037e-07
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
a -2.6106e-01 +/- 1.06e-02
[#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: 4.99474e-07
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
a -3.1481e-01 +/- 1.15e-02
RooFitResult: minimized FCN value: 15739.6, estimated distance to minimum: 3.92419e-05
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
a -3.0010e-01 +/- 1.07e-02
[#0] PROGRESS:Generation -- RooAbsPdf::generateBinned(model) Performing costly accept/reject sampling. If this takes too long, use extended mode to speed up the process.
RooFitResult: minimized FCN value: -1.47174e+08, estimated distance to minimum: 0.162058
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=-1 HESSE=3
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
RooFitResult: minimized FCN value: 3416.14, estimated distance to minimum: 0.000238265
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.9939e+00 +/- 3.63e-02
nbkg 1.0001e+07 +/- 3.24e+03
nsig 9.2425e+03 +/- 6.92e+02
sigma 4.5747e-01 +/- 2.59e-02
Date
January 2023

Definition in file rf614_binned_fit_problems.C.