LHCb Starter Kit 2023: Statistics I (extended version)

Hans Dembinski, TU Dortmund

Scikit-HEP project
Reading simple ROOT files
Histograms
Data visualization
Fitting

About me

Working with the LHCb experiment since 2016
Before astroparticle physicist working with Pierre Auger Observatory and IceCube
Statistics convener in LHCb (2019-2021, 2024-now)
Interests
- Statistics
- High-performance computing
- API design: “Make Interfaces Easy to Use Correctly and Hard to Use Incorrectly” Scott Meyers

Scikit-HEP project

Provides high-performance data analysis tools in Python
Partially funded by NSF grants
Many LHCb developers: Eduardo Rodrigues, Chris Burr, myself, …
Core values
- Specialized tools that integrate well with existing Python ecosystem
- Easy to use for beginners
- Flexible for power users
- Good documentation
- Easy installation

The power of combinatorics

Build everything from basic components that fit together

Image credit: Toni Zaat on Unsplash

Important packages

Python scientific stack
- Numpy: fast computation with arrays
- Scipy: integration, statistical distributions, special functions, …
- Matplotlib: scientific plotting
- Numba: JIT Compiler for Python (very fast)
- Pandas: Convenient processing of tabular data
Scikit-HEP Project
- uproot: fast reading and writing of ROOT Trees
- particle: get particle properties from PDG IDs
- boost-histogram: multi-dimensional generalised histograms
- iminuit: fitting and error computation package
- resample: Easy bootstrapping to compute uncertainty and bias
- pyhf: for histogram fitting/limit setting and preserving likelihoods
Other packages
- jacobi: Easy numerical error propagation
- numba-stats: Fast statistical distributions for fitting
- tabulate: Convert data tables to various formats (LaTeX, HTML, Markdown, …)
- scikit-learn: Basic machine learning tools
- xgboost: Gradient-boosting machine library that won many Kaggle competitions

Install everything with pip install LIBRARY

Working with Jupyter notebooks

I recommend to use vscode
- Very good Jupyter notebook support
- Remote explorer: edit and run code on remote machines

Prepare computing environment

In terminal, create and enter virtual environment (you need Python-3.8 or later)
```
python3 -m venv .venv
source .venv/bin/activate
```

Install packages

python -m pip install --upgrade numba matplotlib uproot boost-histogram iminuit scipy particle numba_stats ipywidgets

uproot provides fast and convenient access to ROOT trees

uproot 5.2.2

name                 | typename                 | interpretation                
---------------------+--------------------------+-------------------------------
trk_len              | int32_t                  | AsDtype('>i4')
mc_trk_len           | int32_t                  | AsDtype('>i4')
trk_imc              | int32_t[]                | AsJagged(AsDtype('>i4'))
trk_px               | float[]                  | AsJagged(AsDtype('>f4'))
trk_py               | float[]                  | AsJagged(AsDtype('>f4'))
trk_pz               | float[]                  | AsJagged(AsDtype('>f4'))
mc_trk_px            | float[]                  | AsJagged(AsDtype('>f4'))
mc_trk_py            | float[]                  | AsJagged(AsDtype('>f4'))
mc_trk_pz            | float[]                  | AsJagged(AsDtype('>f4'))
mc_trk_pid           | int32_t[]                | AsJagged(AsDtype('>i4'))

Tree contains fake simulated LHCb events
Truth in simulation
- mc_trk_len: number of true particles in event
- mc_trk_px: x-component of true momentum of particle
- mc_trk_py: y-component of true momentum of particle
- mc_trk_pz: z-component of true momentum of particle
- mc_trk_pid: particle ID (PDG ID)
Reconstructed simulation
- trk_len: number reconstructed tracks in event
- trk_px: x-component of momentum of reconstructed track
- trk_py: y-component of momentum of reconstructed track
- trk_pz: z-component of momentum of reconstructed track
Special
- trk_imc: index of matched true particle or -1

[6,
 7,
 2,
 6,
 7,
 7,
 6,
 7,
 4,
 6,
 ...,
 2,
 4,
 13,
 3,
 6,
 3,
 4,
 6,
 3]
------------------
type: 1000 * int32

Exercise: Explore some other arrays from event

Some branches are 2D (e.g. trk_px)
First index iterates over events
Second index iterates over tracks per event

0 6 [-0.979, 0.232, -0.464, 0.629, 0.0287, 0.156]
1 7 [-0.59, 0.102, -0.282, -0.585, 0.0525, -0.249, -0.0836]
2 2 [0.124, 0.309]
3 6 [-0.287, 0.0082, -0.671, -0.933, -0.011, 0.0323]
4 7 [-0.086, -0.101, 0.0523, 0.204, -0.0287, 0.115, 0.0954]

Mathematical operators and slicing works on these arrays similar to Numpy

[6,
 7,
 2,
 6,
 7,
 7,
 6,
 7,
 4,
 6]
----------------
type: 10 * int32

[False,
 False,
 False,
 False,
 False,
 False,
 False,
 False,
 False,
 False,
 ...,
 False,
 False,
 True,
 False,
 False,
 False,
 False,
 False,
 False]
-----------------
type: 1000 * bool

Many numpy functions work on awkward arrays

1.24.4

Exercise

Use numpy functions to count how many events have zero tracks
Hint: use np.sum and array mask

'5 events with zero tracks'

For some operations and awkward arrays, we need the awkward library

2.6.1

Exercise * Try to plot trk_px * Can you figure out from the error message why it does not work? * Fix the issue with the function ak.flatten

Reading data from a large tree

Tree often do not fit into memory (tree size > 100 GB and more possible)
Iterate over chunks

Accelerate Python code with Numba

Just-In-Time compiler for scientific Python code
Very powerful and very fast
Fast for loops
Not all Python code supported, works best with numerical code
Only use when needed

numba 0.57.1

let’s subtract reconstructed from true momentum to check momentum resolution
need to select tracks matched to MC particles
this is relatively slow in Python, even if we use array computations

963 ms ± 35.4 ms per loop (mean ± std. dev. of 3 runs, 1 loop each)

4.17 ms ± 1.17 ms per loop (mean ± std. dev. of 3 runs, 1 loop each)

Created array delta_px could still become very large in memory
Better: Create histogram object first and then fill iteratively with chunks

Python frontend to Boost.Histogram library in C++ from the Boost project
Very fast and flexible
Multi-dimensional histograms and other binned statistics (profiles!)
Supports weighted data
And much much more, see docs

boost-histogram 1.4.0

Histogram(Regular(20, -2, 2), storage=Double()) # Sum: 5032.0

                 ┌───────────────────────────────────────────────────────────┐
[-inf,   -2) 0   │                                                           │
[  -2, -1.8) 1   │▏                                                          │
[-1.8, -1.6) 3   │▎                                                          │
[-1.6, -1.4) 17  │█▎                                                         │
[-1.4, -1.2) 25  │█▊                                                         │
[-1.2,   -1) 78  │█████▌                                                     │
[  -1, -0.8) 181 │████████████▉                                              │
[-0.8, -0.6) 308 │█████████████████████▊                                     │
[-0.6, -0.4) 453 │████████████████████████████████▏                          │
[-0.4, -0.2) 647 │█████████████████████████████████████████████▉             │
[-0.2,    0) 819 │██████████████████████████████████████████████████████████ │
[   0,  0.2) 782 │███████████████████████████████████████████████████████▍   │
[ 0.2,  0.4) 629 │████████████████████████████████████████████▌              │
[ 0.4,  0.6) 480 │██████████████████████████████████                         │
[ 0.6,  0.8) 323 │██████████████████████▉                                    │
[ 0.8,    1) 159 │███████████▎                                               │
[   1,  1.2) 77  │█████▌                                                     │
[ 1.2,  1.4) 37  │██▋                                                        │
[ 1.4,  1.6) 7   │▌                                                          │
[ 1.6,  1.8) 3   │▎                                                          │
[ 1.8,    2) 3   │▎                                                          │
[   2,  inf) 0   │                                                           │
                 └───────────────────────────────────────────────────────────┘

Exercise * Use iterative loop to fill h_mc_px and h_delta_px, use calc from before * Plot all three histograms overlaid into one figure with plt.stairs

Histograms can be reduced
Shrink axis range with slices and using loc
Make binning coarser with rebin

Beyond 1D histograms

Supports advanced axis types (e.g. category axis)
Supports higher dimensional histograms
Supports generalized histograms with binned statistics (aka “profile” in ROOT)

Histogram(
  Regular(50, -2, 2),
  IntCategory([-321, 211, 2212, -2212, -211, 321], growth=True),
  storage=Mean()) # Sum: Mean(count=3240, value=0.000433322, variance=0.263444)

We use particle library to look-up properties of particles based on their particle ID
Get name, mass, life-time, quark content, …

particle 0.23.1

Exercise * Only show bins with sufficient entries, use p_2d.counts()

We find that the reconstruction resolution is the same for all particle species
Which means we can remove this axis by summing over its entries (called projection in ROOT)
boost-histogram keeps track of events in overflow and underflow bins by default to make summing exact

Fits

Typical analysis work flow (often automated with Snakemake)
1. Pre-select data
2. Make histograms or profiles from pre-selected data
3. Fit histograms or profiles to extract physical parameters
Many specialized fitting tools for individual purposes, e.g.: scipy.optimize.curve_fit
Generic libraries
- RooFit
- zfit
- pyhf
- iminuit
- …
Fitting (“estimation” in statistics)
- Adjust parametric statistical model to dataset and find model parameters which match dataset best
- AI community calls this “learning”
- Estimate: result of fit
Dataset: samples \(\{ \vec x_i \}\)
- Sample can be one-dimensional or multidimensional
- Values can be continuous (e.g. candidate mass) or discrete (e.g. bin counts)
Model: Probability density or probability mass function \(f(\vec x; \vec p)\) which depends on unknown parameters \(\vec p\)
Conjecture: maximum-likelihood estimate (MLE) \(\hat{\vec p} = \text{argmax}_{\vec p} L(\vec p)\) is optimal, with \(L(\vec p) = \prod_i f(\vec x_i; \vec p)\)
- Proven: Asymptotically optimal, unbiased with minimum variance
  - No proof for finite data sets, but unchallenged
- Universal: MLE can be computed for any problem
- Equivalent: minimize negative log-likelihood (NLL), \(\hat{\vec p} = \text{argmin}_{\vec p}(-\ln L(\vec p))\)
Limiting case of maximum-likelihood fit: least-squares fit aka chi-square fit

Python frontend to Minuit2 C++ library maintained by ROOT team at CERN
Maximum-likelihood and least-squares fits with error estimation
Easy to use for beginners, flexible for power users
Comes with builtin cost functions for common fits
Large number of tutorials

iminuit 2.25.2

iminuit vs. other packages

pyhf
- Better choice when you need to compute limits and preserve the likelihood
zfit and RooFit
- Help you build statistical models: automatic normalization & convolutions
  - This has to be done “by hand” in iminuit
- Good for standard tasks, but restrict your freedom
iminuit
- Has good documentation
- Well-designed API
- Is very flexible and open
- Is very fast with numba
Typical fit in particle physics: mass distribution of decay candidates with signal peak over smooth background
Task: get the number of signal candidates = signal yield
Can be solved with a fit

To apply maximum-likelihood method, we need statistical model
Full sample is additive mixture of signal sample and background sample
Assumption 1: signal is normal-distributed; pdf is \(\mathcal{N}(\mu, \sigma)\) with parameters \(\mu\) and \(\sigma\)
Assumption 2: background is truncated exponential distribution

Migrad
FCN = -1790	Nfcn = 209
EDM = 762 (Goal: 0.0001)
INVALID Minimum	ABOVE EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	0.500	0.012
1	mu	0.500	0.003
2	sigma	0.155	0.005
3	slope	1.00	0.07

	z	mu	sigma	slope
z	0.000149	-1e-6 (-0.037)	0.033e-3 (0.508)	-0.25e-3 (-0.304)
mu	-1e-6 (-0.037)	9.19e-06	-0e-6 (-0.025)	-68e-6 (-0.335)
sigma	0.033e-3 (0.508)	-0e-6 (-0.025)	2.91e-05	-0.053e-3 (-0.147)
slope	-0.25e-3 (-0.304)	-68e-6 (-0.335)	-0.053e-3 (-0.147)	0.00442

Fit did not work
Plotting 1D or 2D slices of likelihood may help

Increase print level

W FCN result is NaN for [ 0.406325 0.495731 -0.0437361 1.13819 ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W VariableMetricBuilder No improvement in line search
W VariableMetricBuilder Iterations finish without convergence; Edm 761.614 Requested 0.0001
W VariableMetricBuilder No convergence; Edm 761.614 is above tolerance 0.001
W FCN result is NaN for [ 0.406325 0.495731 -0.0437361 1.13819 ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W VariableMetricBuilder No improvement in line search
W VariableMetricBuilder Iterations finish without convergence; Edm 761.614 Requested 0.0001
W VariableMetricBuilder No convergence; Edm 761.614 is above tolerance 0.001
W FCN result is NaN for [ 0.406325 0.495731 -0.0437361 1.13819 ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W VariableMetricBuilder No improvement in line search
W VariableMetricBuilder Iterations finish without convergence; Edm 761.614 Requested 0.0001
W VariableMetricBuilder No convergence; Edm 761.614 is above tolerance 0.001
W FCN result is NaN for [ 0.406325 0.495731 -0.0437361 1.13819 ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W VariableMetricBuilder No improvement in line search
W VariableMetricBuilder Iterations finish without convergence; Edm 761.614 Requested 0.0001
W VariableMetricBuilder No convergence; Edm 761.614 is above tolerance 0.001
W FCN result is NaN for [ 0.406325 0.495731 -0.0437361 1.13819 ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W FCN result is NaN for [ nan nan nan nan ]
W VariableMetricBuilder No improvement in line search
W VariableMetricBuilder Iterations finish without convergence; Edm 761.614 Requested 0.0001
W VariableMetricBuilder No convergence; Edm 761.614 is above tolerance 0.001

fit failed because parameters need limits

Exercise

Repeat fit, but with parameter limits.
Set parameter limits by calling m.limits["parameter"] = (lower, upper) before m.migrad()
You can set identical limits for multiple parameters, for example, m.limits["a", "b"] = (0, 1)

Migrad
FCN = -2095	Nfcn = 481
EDM = 1.78e-06 (Goal: 0.0001)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	0.498	0.009			0	1
1	mu	0.4968	0.0020			0	1
2	sigma	0.0994	0.0020			0
3	slope	1.03	0.06			0

	z	mu	sigma	slope
z	8.94e-05	-1e-6 (-0.063)	10e-6 (0.547)	-0.12e-3 (-0.213)
mu	-1e-6 (-0.063)	4.21e-06	-0e-6 (-0.061)	-27e-6 (-0.225)
sigma	10e-6 (0.547)	-0e-6 (-0.061)	3.94e-06	-18e-6 (-0.156)
slope	-0.12e-3 (-0.213)	-27e-6 (-0.225)	-18e-6 (-0.156)	0.00345

compute asymmetric errors with Minos algorithm (profile likelihood method)

Migrad
FCN = -2095	Nfcn = 692
EDM = 1.78e-06 (Goal: 0.0001)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	0.498	0.009	-0.009	0.009	0	1
1	mu	0.4968	0.0021	-0.0020	0.0020	0	1
2	sigma	0.0994	0.0020	-0.0019	0.0020	0
3	slope	1.03	0.06	-0.06	0.06	0

	z		mu		sigma		slope
Error	-0.009	0.009	-0.002	0.002	-0.0019	0.0020	-0.06	0.06
Valid	True	True	True	True	True	True	True	True
At Limit	False	False	False	False	False	False	False	False
Max FCN	False	False	False	False	False	False	False	False
New Min	False	False	False	False	False	False	False	False

	z	mu	sigma	slope
z	8.94e-05	-1e-6 (-0.063)	10e-6 (0.547)	-0.12e-3 (-0.213)
mu	-1e-6 (-0.063)	4.21e-06	-0e-6 (-0.061)	-27e-6 (-0.225)
sigma	10e-6 (0.547)	-0e-6 (-0.061)	3.94e-06	-18e-6 (-0.156)
slope	-0.12e-3 (-0.213)	-27e-6 (-0.225)	-18e-6 (-0.156)	0.00345

you can compute 2D confidence regions with mncontour and draw them with draw_mncontour

Exercise

Draw 1,2,3 sigma contours, use keyword cl

for nice overview: compute matrix of likelihood profiles with draw_mnmatrix

Writing unbinned negative-log-likelihood is repetitive, only model ever changes
Use built-in cost function from iminuit to get nice extras

Migrad
FCN = -4189	Nfcn = 169
EDM = 8.37e-05 (Goal: 0.0002)	time = 0.4 sec
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	0.498	0.009			0	1
1	mu	0.4968	0.0020			0	1
2	sigma	0.0994	0.0020			0
3	slope	1.03	0.06			0

	z	mu	sigma	slope
z	8.62e-05	-1e-6 (-0.037)	10e-6 (0.529)	-0.11e-3 (-0.211)
mu	-1e-6 (-0.037)	4.17e-06	-0e-6 (-0.039)	-28e-6 (-0.230)
sigma	10e-6 (0.529)	-0e-6 (-0.039)	3.83e-06	-17e-6 (-0.152)
slope	-0.11e-3 (-0.211)	-28e-6 (-0.230)	-17e-6 (-0.152)	0.00345

Exercise

Use interactive mode via m.interactive() to check how fit reacts to parameter changes
Useful for debugging and to find starting values by hand

Faster fits?

Option 1: accelerate NLL computation with Numba
Option 2: use binned data

Migrad
FCN = -4.216e+05	Nfcn = 169
EDM = 0.000135 (Goal: 0.0002)	time = 5.9 sec
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	499.7e-3	0.9e-3			0	1
1	mu	499.98e-3	0.20e-3			0	1
2	sigma	99.59e-3	0.20e-3			0
3	slope	0.997	0.006			0

	z	mu	sigma	slope
z	8.87e-07	-0.01e-6 (-0.049)	0.10e-6 (0.548)	-1.2e-6 (-0.229)
mu	-0.01e-6 (-0.049)	4.18e-08	-0 (-0.060)	-0.26e-6 (-0.226)
sigma	0.10e-6 (0.548)	-0 (-0.060)	4.09e-08	-0.18e-6 (-0.164)
slope	-1.2e-6 (-0.229)	-0.26e-6 (-0.226)	-0.18e-6 (-0.164)	3.06e-05

General insights * Time spend in minimizer is negligible * Bottleneck is evaluating NLL function * Bottleneck inside NLL function is usually evaluation of model * NLL computation is trivially parallelizable: log(pdf) computed independently for each data point * Numba allows us to exploit auto-vectorization (SIMD instructions) and parallel computing on multiple cors

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Migrad
FCN = -2.108e+05	Nfcn = 242
EDM = 2.48e-06 (Goal: 0.0001)	time = 1.4 sec
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	499.7e-3	0.9e-3			0	1
1	mu	499.98e-3	0.20e-3			0	1
2	sigma	99.59e-3	0.20e-3			0
3	slope	0.997	0.006			0

	z	mu	sigma	slope
z	8.87e-07	-0.01e-6 (-0.049)	0.10e-6 (0.548)	-1.2e-6 (-0.229)
mu	-0.01e-6 (-0.049)	4.18e-08	-0 (-0.060)	-0.26e-6 (-0.226)
sigma	0.10e-6 (0.548)	-0 (-0.060)	4.09e-08	-0.18e-6 (-0.164)
slope	-1.2e-6 (-0.229)	-0.26e-6 (-0.226)	-0.18e-6 (-0.164)	3.06e-05

Often simpler: fit histograms
No bias from fitting histograms if done correctly, only loss in precision
Need model cdf instead of model pdf: \(F(x; \vec p) = \int_{-\infty}^x f(x'; \vec p) \text{d}x'\)
Recommended to use iminuit.cost.BinnedNLL, comes with many useful features

Migrad
FCN = 84.34 (χ²/ndof = 0.9)	Nfcn = 162
EDM = 1.36e-07 (Goal: 0.0002)	time = 0.4 sec
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	z	499.7e-3	0.9e-3
1	mu	499.99e-3	0.20e-3
2	sigma	99.58e-3	0.20e-3
3	slope	0.997	0.006

	z	mu	sigma	slope
z	8.98e-07	-0.01e-6 (-0.049)	0.10e-6 (0.535)	-1.2e-6 (-0.236)
mu	-0.01e-6 (-0.049)	4.22e-08	-0 (-0.052)	-0.26e-6 (-0.228)
sigma	0.10e-6 (0.535)	-0 (-0.052)	4.15e-08	-0.18e-6 (-0.161)
slope	-1.2e-6 (-0.236)	-0.26e-6 (-0.228)	-0.18e-6 (-0.161)	3.09e-05

See study of binned fits vs. unbinned fit for detailed comparison
With sufficient many bins (>= 20 in this case), increase in parameter uncertainties remains negligible
Binned fits are (100-200)x faster than unbinned fits without using numba

How well is model matching data?

Judgement requires binning of data
Look at pull distribution
- Compute (data - model) / error for each bin
Compare \(\chi^2\) value (goodness-of-fit test statistic) against the degrees of freedom
- Simple check: reduced chi-square \(\chi^2/n_\text{dof}\) should be about 1
- Better check: p-value \(= \int_{\chi^2_\text{observed}}^{\infty} f(\chi^2; n_\text{dof}) \, \text{d}\chi^2\) (survival function)
- If p-value very small (< 1%), wrong model or underestimated errors

'chi2/ndof = 0.88 p-value = 0.80'

Exercise

Plot histogram of pulls, should look like standard normal distribution
Compute mean and standard deviation (use np.mean and np.std)

Other fits

To fit yields (counts) instead of fractions, use ExtendedUnbinnedNLL or ExtendedBinnedNLL
To fit (x, y) pairs, use LeastSquares
To fit simulated templates, use Template
See documentation how models should be constructed then

Template fits

Sometimes PDF/CDF not known in parametric form, but can be Monte-Carlo simulated
Idea: fit a histogram of simulation output (template) to data histogram
Additional source of uncertainty from finite size of simulation
Correct error propagation with Barlow-Beeston method, Barlow and Beeston, Comput.Phys.Commun. 77 (1993) 219-228
iminuit.cost.Template
- Extension to weighted data and simuation, H. Dembinski, A. Abdelmotteleb, Eur.Phys.J.C 82 (2022) 11, 1043
- Approximation to Barlow-Beeston method for unweighted data
- Faster computation
- Allows you to mix parametric models and templates

Migrad
FCN = 8.352 (χ²/ndof = 0.6)	Nfcn = 132
EDM = 2.75e-07 (Goal: 0.0002)	time = 0.5 sec
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	770	70			0
1	x1	220	40			0

	x0	x1
x0	4.3e+03	-0.8e3 (-0.354)
x1	-0.8e3 (-0.354)	1.32e+03

Exercise

With make_data2 generate templates with 1 000 000 simulated points
Fit this data with the Template cost function
Compare parameter uncertainties

Migrad
FCN = 4.329 (χ²/ndof = 0.3)	Nfcn = 110
EDM = 2.66e-06 (Goal: 0.0002)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	820	32			0
1	x1	201	20			0

	x0	x1
x0	1.03e+03	-0.2e3 (-0.326)
x1	-0.2e3 (-0.326)	415

With Template you can also fit a mix of template and parametric model
Often template only needed for background(s)

Migrad
FCN = 2.526 (χ²/ndof = 0.2)	Nfcn = 106
EDM = 1.52e-06 (Goal: 0.0002)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	819	34			0
1	x1_s	202	23
2	x1_mu	0.492	0.006			0	1
3	x1_sigma	0.049	0.007			0

	x0	x1_s	x1_mu	x1_sigma
x0	1.16e+03	-0.3e3 (-0.428)	8.13e-3 (0.038)	-79.31e-3 (-0.330)
x1_s	-0.3e3 (-0.428)	541	-8.17e-3 (-0.056)	79.27e-3 (0.483)
x1_mu	8.13e-3 (0.038)	-8.17e-3 (-0.056)	3.99e-05	-0 (-0.069)
x1_sigma	-79.31e-3 (-0.330)	79.27e-3 (0.483)	-0 (-0.069)	4.98e-05

Overfitting

Do not use models with degenerate parameters
Avoid overfitting

Migrad
FCN = 23.81 (χ²/ndof = 4.8)	Nfcn = 51
EDM = 4.25e-21 (Goal: 0.0002)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	0.950	0.011
1	x1	0.920	0.027
2	x2	0.21	0.04

	x0	x1	x2
x0	0.000127	-0.21e-3 (-0.680)	0.14e-3 (0.321)
x1	-0.21e-3 (-0.680)	0.000737	-0.7e-3 (-0.666)
x2	0.14e-3 (0.321)	-0.7e-3 (-0.666)	0.00147

Exercise

Fit with increasing number of parameters (increase the length of the starting parameters)
- What do you observe? Look out for changes in chi2/ndof and the parameter correlations

Try fitting this model

def model(x, p):
    return p[0] + bernstein.density(x, p[1:], -0.4, 0.4)

What do you observe?

Migrad
FCN = 1.069 (χ²/ndof = 0.5)	Nfcn = 152
EDM = 2.41e-15 (Goal: 0.0002)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	0.956	0.026
1	x1	1.02	0.13
2	x2	0.47	0.27
3	x3	1.30	0.34
4	x4	0.29	0.25
5	x5	0.10	0.08

	x0	x1	x2	x3	x4	x5
x0	0.000693	-3.0e-3 (-0.885)	5.7e-3 (0.800)	-6.2e-3 (-0.687)	3.6e-3 (0.551)	-0.7e-3 (-0.307)
x1	-3.0e-3 (-0.885)	0.0161	-0.033 (-0.961)	0.037 (0.852)	-0.022 (-0.698)	0.004 (0.394)
x2	5.7e-3 (0.800)	-0.033 (-0.961)	0.0721	-0.09 (-0.953)	0.05 (0.822)	-0.010 (-0.480)
x3	-6.2e-3 (-0.687)	0.037 (0.852)	-0.09 (-0.953)	0.116	-0.08 (-0.938)	0.016 (0.583)
x4	3.6e-3 (0.551)	-0.022 (-0.698)	0.05 (0.822)	-0.08 (-0.938)	0.062	-0.014 (-0.721)
x5	-0.7e-3 (-0.307)	0.004 (0.394)	-0.010 (-0.480)	0.016 (0.583)	-0.014 (-0.721)	0.00649

Other Formats

About me

Scikit-HEP project

The power of combinatorics

Important packages

Working with Jupyter notebooks

Prepare computing environment

Reading data from a large tree

Accelerate Python code with Numba

Beyond 1D histograms

Fits

iminuit vs. other packages

Faster fits?

How well is model matching data?

Other fits

Template fits

Overfitting