Logo ROOT  

Reference Guide
Loading...
Searching...
No Matches
rf611_weightedfits.C File Reference
Tutorials » RooFit packages tutorials » RooFit

Detailed Description

View in nbviewer Open in SWAN
Likelihood and minimization: Parameter uncertainties for weighted unbinned ML fits

Parameter uncertainties for weighted unbinned ML fits

Based on example from https://arxiv.org/abs/1911.01303

This example compares different approaches to determining parameter uncertainties in weighted unbinned maximum likelihood fits. Performing a weighted unbinned maximum likelihood fits can be useful to account for acceptance effects and to statistically subtract background events using the sPlot formalism. It is however well known that the inverse Hessian matrix does not yield parameter uncertainties with correct coverage in the presence of event weights. Three approaches to the determination of parameter uncertainties are compared in this example:

  1. Using the inverse weighted Hessian matrix [SumW2Error(false)]
  2. Using the expression [SumW2Error(true)]

    \[ V_{ij} = H_{ik}^{-1} C_{kl} H_{lj}^{-1} \]

    where H is the weighted Hessian matrix and C is the Hessian matrix with squared weights
  3. The asymptotically correct approach (for details please see https://arxiv.org/abs/1911.01303) [Asymptotic(true)]

    \[ V_{ij} = H_{ik}^{-1} D_{kl} H_{lj}^{-1} \]

    where H is the weighted Hessian matrix and D is given by

    \[ D_{kl} = \sum_{e=1}^{N} w_e^2 \frac{\partial \log(P)}{\partial \lambda_k}\frac{\partial \log(P)}{\partial \lambda_l} \]

    with the event weight \(w_e\).

The example performs the fit of a second order polynomial in the angle cos(theta) [-1,1] to a weighted data set. The polynomial is given by

\[ P = \frac{ 1 + c_0 \cdot \cos(\theta) + c_1 \cdot \cos(\theta) \cdot \cos(\theta) }{\mathrm{Norm}} \]

The two coefficients \( c_0 \) and \( c_1 \) and their uncertainties are to be determined in the fit.

The per-event weight is used to correct for an acceptance effect, two different acceptance models can be studied:

  • acceptancemodel==1: eff = \( 0.3 + 0.7 \cdot \cos(\theta) \cdot \cos(\theta) \)
  • acceptancemodel==2: eff = \( 1.0 - 0.7 \cdot \cos(\theta) \cdot \cos(\theta) \) The data is generated to be flat before the acceptance effect.

The performance of the different approaches to determine parameter uncertainties is compared using the pull distributions from a large number of pseudoexperiments. The pull is defined as \( (\lambda_i - \lambda_{gen})/\sigma(\lambda_i) \), where \( \lambda_i \) is the fitted parameter and \( \sigma(\lambda_i) \) its uncertainty for pseudoexperiment number i. If the fit is unbiased and the parameter uncertainties are estimated correctly, the pull distribution should be a Gaussian centered around zero with a width of one.

#include "TH1D.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TLegend.h"
#include "RooRealVar.h"
#include "RooFitResult.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
using namespace RooFit;
void rf611_weightedfits(int acceptancemodel = 2)
{
// I n i t i a l i s a t i o n a n d S e t u p
//------------------------------------------------
// plotting options
gStyle->SetPaintTextFormat(".1f");
gStyle->SetEndErrorSize(6.0);
gStyle->SetTitleSize(0.05, "XY");
gStyle->SetLabelSize(0.05, "XY");
gStyle->SetTitleOffset(0.9, "XY");
gStyle->SetTextSize(0.05);
gStyle->SetPadLeftMargin(0.125);
gStyle->SetPadBottomMargin(0.125);
gStyle->SetPadTopMargin(0.075);
gStyle->SetPadRightMargin(0.075);
gStyle->SetMarkerStyle(20);
gStyle->SetMarkerSize(1.0);
gStyle->SetHistLineWidth(2.0);
gStyle->SetHistLineColor(1);
// initialise TRandom3
TRandom3 *rnd = new TRandom3();
rnd->SetSeed(191101303);
// accepted events and events weighted to account for the acceptance
TH1 *haccepted = new TH1D("haccepted", "Generated events;cos(#theta);#events", 40, -1.0, 1.0);
TH1 *hweighted = new TH1D("hweighted", "Generated events;cos(#theta);#events", 40, -1.0, 1.0);
// histograms holding pull distributions
std::array<TH1 *, 3> hc0pull;
std::array<TH1 *, 3> hc1pull;
std::array<TH1 *, 3> hntotpull;
std::array<std::string, 3> methodLabels{"Inverse weighted Hessian matrix [SumW2Error(false)]",
"Hessian matrix with squared weights [SumW2Error(true)]",
"Asymptotically correct approach [Asymptotic(true)]"};
auto makePullXLabel = [](std::string const &pLabel) {
return "Pull (" + pLabel + "^{fit}-" + pLabel + "^{gen})/#sigma(" + pLabel + ")";
};
for (std::size_t i = 0; i < 3; ++i) {
std::string const &iLabel = std::to_string(i);
// using the inverse Hessian matrix
std::string hc0XLabel = methodLabels[i] + ";" + makePullXLabel("c_{0}") + ";";
std::string hc1XLabel = methodLabels[i] + ";" + makePullXLabel("c_{1}") + ";";
std::string hntotXLabel = methodLabels[i] + ";" + makePullXLabel("N_{tot}") + ";";
hc0pull[i] = new TH1D(("hc0pull" + iLabel).c_str(), hc0XLabel.c_str(), 20, -5.0, 5.0);
// using the correction with the Hessian matrix with squared weights
hc1pull[i] = new TH1D(("hc1pull" + iLabel).c_str(), hc1XLabel.c_str(), 20, -5.0, 5.0);
// asymptotically correct approach
hntotpull[i] = new TH1D(("hntotpull" + iLabel).c_str(), hntotXLabel.c_str(), 20, -5.0, 5.0);
}
// number of pseudoexperiments (toys) and number of events per pseudoexperiment
constexpr std::size_t ntoys = 500;
constexpr std::size_t nstats = 500;
// parameters used in the generation
constexpr double c0gen = 0.0;
constexpr double c1gen = 0.0;
// Silence fitting and minimisation messages
auto &msgSv = RooMsgService::instance();
msgSv.getStream(1).removeTopic(RooFit::Minimization);
msgSv.getStream(1).removeTopic(RooFit::Fitting);
std::cout << "Running " << ntoys * 3 << " toy fits ..." << std::endl;
// M a i n l o o p : r u n p s e u d o e x p e r i m e n t s
//----------------------------------------------------------------
for (std::size_t i = 0; i < ntoys; i++) {
// S e t u p p a r a m e t e r s a n d P D F
//-----------------------------------------------
// angle theta and the weight to account for the acceptance effect
RooRealVar costheta("costheta", "costheta", -1.0, 1.0);
RooRealVar weight("weight", "weight", 0.0, 1000.0);
// initialise parameters to fit
RooRealVar c0("c0", "0th-order coefficient", c0gen, -1.0, 1.0);
RooRealVar c1("c1", "1st-order coefficient", c1gen, -1.0, 1.0);
c0.setError(0.01);
c1.setError(0.01);
// create simple second-order polynomial as probability density function
RooPolynomial pol("pol", "pol", costheta, {c0, c1}, 1);
double ngen = nstats;
if (acceptancemodel == 1)
ngen *= 2.0 / (23.0 / 15.0);
else
ngen *= 2.0 / (16.0 / 15.0);
RooRealVar ntot("ntot", "ntot", ngen, 0.0, 2.0 * ngen);
RooExtendPdf extended("extended", "extended pdf", pol, ntot);
int npoisson = rnd->Poisson(nstats);
// G e n e r a t e d a t a s e t f o r p s e u d o e x p e r i m e n t i
//-------------------------------------------------------------------------------
RooDataSet data("data", "data", {costheta, weight}, WeightVar("weight"));
// generate nstats events
for (std::size_t j = 0; j < npoisson; j++) {
bool finished = false;
// use simple accept/reject for generation
while (!finished) {
costheta = 2.0 * rnd->Rndm() - 1.0;
// efficiency for the specific value of cos(theta)
double eff = 1.0;
if (acceptancemodel == 1)
eff = 1.0 - 0.7 * costheta.getVal() * costheta.getVal();
else
eff = 0.3 + 0.7 * costheta.getVal() * costheta.getVal();
// use 1/eff as weight to account for acceptance
weight = 1.0 / eff;
// accept/reject
if (10.0 * rnd->Rndm() < eff * pol.getVal())
finished = true;
}
haccepted->Fill(costheta.getVal());
hweighted->Fill(costheta.getVal(), weight.getVal());
data.add({costheta, weight}, weight.getVal());
}
auto fillPulls = [&](std::size_t i) {
hc0pull[i]->Fill((c0.getVal() - c0gen) / c0.getError());
hc1pull[i]->Fill((c1.getVal() - c1gen) / c1.getError());
hntotpull[i]->Fill((ntot.getVal() - ngen) / ntot.getError());
};
// F i t t o y u s i n g t h e t h r e e d i f f e r e n t a p p r o a c h e s t o u n c e r t a i
// n t y d e t e r m i n a t i o n
//-------------------------------------------------------------------------------------------------------------------------------------------------
// this uses the inverse weighted Hessian matrix
extended.fitTo(data, SumW2Error(false), PrintLevel(-1));
fillPulls(0);
// this uses the correction with the Hesse matrix with squared weights
extended.fitTo(data, SumW2Error(true), PrintLevel(-1));
fillPulls(1);
// this uses the asymptotically correct approach
extended.fitTo(data, AsymptoticError(true), PrintLevel(-1));
fillPulls(2);
}
std::cout << "... done." << std::endl;
// P l o t o u t p u t d i s t r i b u t i o n s
//--------------------------------------------------
// plot accepted (weighted) events
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);
TCanvas *cevents = new TCanvas("cevents", "cevents", 800, 600);
cevents->cd(1);
hweighted->SetMinimum(0.0);
hweighted->SetLineColor(2);
hweighted->Draw("hist");
haccepted->Draw("same hist");
TLegend *leg = new TLegend(0.6, 0.8, 0.9, 0.9);
leg->AddEntry(haccepted, "Accepted");
leg->AddEntry(hweighted, "Weighted");
leg->Draw();
cevents->Update();
// plot pull distributions
TCanvas *cpull = new TCanvas("cpull", "cpull", 1200, 800);
cpull->Divide(3, 3);
std::vector<TH1 *> pullHistos{hc0pull[0], hc0pull[1], hc0pull[2], hc1pull[0], hc1pull[1],
hc1pull[2], hntotpull[0], hntotpull[1], hntotpull[2]};
gStyle->SetOptStat(1100);
gStyle->SetOptFit(11);
for (std::size_t i = 0; i < pullHistos.size(); ++i) {
cpull->cd(i + 1);
pullHistos[i]->Fit("gaus");
pullHistos[i]->Draw("ep");
}
cpull->Update();
}
gStyle
externTStyle * gStyle
Definition TStyle.h:442
RooDataSet
Container class to hold unbinned data.
Definition RooDataSet.h:32
RooExtendPdf
RooExtendPdf is a wrapper around an existing PDF that adds a parameteric extended likelihood term to ...
Definition RooExtendPdf.h:22
RooMsgService::instance
static RooMsgService & instance()
Return reference to singleton instance.
Definition RooMsgService.cxx:344
RooPolynomial
RooPolynomial implements a polynomial p.d.f of the form.
Definition RooPolynomial.h:25
RooRealVar
Variable that can be changed from the outside.
Definition RooRealVar.h:37
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:44
TCanvas
The Canvas class.
Definition TCanvas.h:23
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:716
TCanvas::Update
void Update() override
Update canvas pad buffers.
Definition TCanvas.cxx:2486
TH1D
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:926
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3393
TH1::Draw
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3097
TH1::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:653
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1295
TRandom3
Random number generator class based on M.
Definition TRandom3.h:27
TRandom3::Rndm
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom3.cxx:107
TRandom3::SetSeed
void SetSeed(ULong_t seed=0) override
Set the random generator sequence.
Definition TRandom3.cxx:216
TRandom::Poisson
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:403
RooFit::WeightVar
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
Definition RooGlobalFunc.cxx:398
RooFit::AsymptoticError
RooCmdArg AsymptoticError(bool flag)
Definition RooGlobalFunc.cxx:739
RooFit::SumW2Error
RooCmdArg SumW2Error(bool flag)
Definition RooGlobalFunc.cxx:735
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition RooGlobalFunc.cxx:659
c1
return c1
Definition legend1.C:41
leg
leg
Definition legend1.C:34
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:72
rf503_wspaceread.data
data
Definition rf503_wspaceread.py:37
Running 1500 toy fits ...
... done.
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 12.2745
NDf = 12
Edm = 3.84566e-06
NCalls = 53
Constant = 82.0512 +/- 4.52865
Mean = -0.042844 +/- 0.0572716
Sigma = 1.19191 +/- 0.0392558 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 8.83881
NDf = 9
Edm = 6.19877e-08
NCalls = 53
Constant = 104.798 +/- 5.86613
Mean = -0.0132039 +/- 0.0432773
Sigma = 0.938314 +/- 0.0320986 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 7.28099
NDf = 10
Edm = 6.30801e-07
NCalls = 53
Constant = 103.063 +/- 5.55774
Mean = -0.0233563 +/- 0.0435444
Sigma = 0.954358 +/- 0.0285394 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 23.1247
NDf = 14
Edm = 2.27317e-06
NCalls = 53
Constant = 69.5383 +/- 3.95119
Mean = -0.160071 +/- 0.0652133
Sigma = 1.37036 +/- 0.0473756 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 7.18474
NDf = 15
Edm = 9.01145e-06
NCalls = 61
Constant = 51.3838 +/- 3.49309
Mean = -0.985792 +/- 0.0760226
Sigma = 1.36366 +/- 0.0605719 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 11.9456
NDf = 11
Edm = 4.88766e-08
NCalls = 61
Constant = 94.7614 +/- 5.45747
Mean = -0.0978895 +/- 0.0483575
Sigma = 1.03309 +/- 0.0382346 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 20.9563
NDf = 16
Edm = 8.77656e-06
NCalls = 53
Constant = 69.8331 +/- 3.99247
Mean = -0.0767062 +/- 0.0640105
Sigma = 1.36919 +/- 0.0474577 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 8.84994
NDf = 10
Edm = 1.85865e-06
NCalls = 53
Constant = 101.665 +/- 5.77203
Mean = -0.0558063 +/- 0.0444665
Sigma = 0.96468 +/- 0.0336688 (limited)
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 9.73174
NDf = 10
Edm = 1.95484e-06
NCalls = 53
Constant = 99.7164 +/- 5.69932
Mean = -0.0653755 +/- 0.0459247
Sigma = 0.982023 +/- 0.0349408 (limited)
Date
November 2019
Author
Christoph Langenbruch

Definition in file rf611_weightedfits.C.