Introduction to zfit¶

In this notebook, we will have a walk through the main components of zfit and their features. Especially the extensive model building part will be discussed separately. zfit consists of 5 mostly independent parts. Other libraries can rely on this parts to do plotting or statistical inference, such as hepstats does. Therefore we will discuss two libraries in this tutorial: zfit to build models, data and a loss, minimize it and get a fit result and hepstats, to use the loss we built here and do inference.

Data¶

This component in general plays a minor role in zfit: it is mostly to provide a unified interface for data.

Preprocessing is therefore not part of zfit and should be done beforehand. Python offers many great possibilities to do so (e.g. Pandas).

zfit Data can load data from various sources, most notably from Numpy, Pandas DataFrame, TensorFlow Tensor and ROOT (using uproot). It is also possible, for convenience, to convert it directly to_pandas. The constructors are named from_numpy, from_root etc.

import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import zfit
from zfit import z

A Data needs not only the data itself but also the observables: the human readable string identifiers of the axes (corresponding to “columns” of a Pandas DataFrame). It is convenient to define the Space not only with the observable but also with a limit: this can directly be re-used as the normalization range in the PDF.

First, let’s define our observables

obs = zfit.Space('obs1', (-5, 10))

This Space has limits. Next to the effect of handling the observables, we can also play with the limits: multiple Spaces can be added to provide disconnected ranges. More importantly, Space offers functionality:

limit1d: return the lower and upper limit in the 1 dimensional case (raises an error otherwise)
rect_limits: return the n dimensional limits
area(): calculate the area (e.g. distance between upper and lower)
inside(): return a boolean Tensor corresponding to whether the value is inside the Space
filter(): filter the input values to only return the one inside

size_normal = 10000
data_normal_np = np.random.normal(size=size_normal, scale=2)

data_normal = zfit.Data.from_numpy(obs=obs, array=data_normal_np)

The main functionality is

nevents: attribute that returns the number of events in the object
data_range: a Space that defines the limits of the data; if outside, the data will be cut
n_obs: defines the number of dimensions in the dataset
with_obs: returns a subset of the dataset with only the given obs
weights: event based weights

Furthermore, value returns a Tensor with shape (nevents, n_obs).

To retrieve values, in general z.unstack_x(data) should be used; this returns a single Tensor with shape (nevents) or a list of tensors if n_obs is larger then 1.

print(f"We have {data_normal.nevents} events in our dataset with the minimum of {np.min(data_normal.unstack_x())}")  # remember! The obs cut out some of the data

We have 9931 events in our dataset with the minimum of -4.994790077305963

data_normal.n_obs

Model¶

Building models is by far the largest part of zfit. We will therefore cover an essential part, the possibility to build custom models, in an extra chapter. Let’s start out with the idea that you define your parameters and your observable space; the latter is the expected input data.

There are two types of models in zfit:

functions, which are rather simple and “underdeveloped”; their usage is often not required.
PDF that are function which are normalized (over a specified range); this is the main model and is what we gonna use throughout the tutorials.

A PDF is defined by

(1)¶\[\begin{align} \mathrm{PDF}_{f(x)}(x; \theta) = \frac{f(x; \theta)}{\int_{a}^{b} f(x; \theta)} \end{align}\]

where a and b define the normalization range (norm_range), over which (by inserting into the above definition) the integral of the PDF is unity.

zfit has a modular approach to things and this is also true for models. While the normalization itself (e.g. what are parameters, what is normalized data) will already be pre-defined in the model, models are composed of functions that are transparently called inside. For example, a Gaussian would usually be implemented by writing a Python function def gauss(x, mu, sigma), which does not care about the normalization and then be wrapped in a PDF, where the normalization and what is a parameter is defined.

In principle, we can go far by using simply functions (e.g. TensorFlowAnalysis/AmpliTF by Anton Poluektov uses this approach quite successfully for Amplitude Analysis), but this design has limitations for a more general fitting library such as zfit (or even TensorWaves, being built on top of AmpliTF). The main thing is to keep track of the different ordering of the data and parameters, especially the dependencies.

Let’s create a simple Gaussian PDF. We already defined the Space for the data before, now we only need the parameters. This are a different object than a Space.

Parameter¶

A Parameter (there are different kinds actually, more on that later) takes the following arguments as input: Parameter(human readable name, initial value[, lower limit, upper limit]) where the limits are recommended but not mandatory. Furthermore, step_size can be given (which is useful to be around the given uncertainty, e.g. for large yields or small values it can help a lot to set this). Also, a floating argument is supported, indicating whether the parameter is allowed to float in the fit or not (just omitting the limits does not make a parameter constant).

Parameters have a unique name. This is served as the identifier for e.g. fit results. However, a parameter cannot be retrieved by its string identifier (its name) but the object itself should be used. In places where a parameter maps to something, the object itself is needed, not its name.

mu = zfit.Parameter('mu', 1, -3, 3, step_size=0.2)
sigma_num = zfit.Parameter('sigma42', 1, 0.1, 10, floating=False)

These attributes can be changed:

print(f"sigma is float: {sigma_num.floating}")
sigma_num.floating = True
print(f"sigma is float: {sigma_num.floating}")

sigma is float: False

sigma is float: True

PITFALL NOTEBOOKS: since the parameters have a unique name, a second parameter with the same name cannot be created; the behavior is undefined and therefore it raises an error. While this does not pose a problem in a normal Python script, it does in a Jupyter-like notebook, since it is an often practice to “rerun” a cell as an attempt to “reset” things. Bear in mind that this does not make sense, from a logic point of view. The parameter already exists. Best practice: write a small wrapper, do not rerun the parameter creation cell or simply rerun the notebook (restart kernel & run all). For further details, have a look at the discussion and arguments here

Now we have everything to create a Gaussian PDF:

gauss = zfit.pdf.Gauss(obs=obs, mu=mu, sigma=sigma_num)

Since this holds all the parameters and the observables are well defined, we can retrieve them

gauss.n_obs  # dimensions

gauss.obs

('obs1',)

gauss.space

<zfit Space obs=('obs1',), axes=(0,), limits=(array([[-5.]]), array([[10.]]))>

gauss.norm_range

<zfit Space obs=('obs1',), axes=(0,), limits=(array([[-5.]]), array([[10.]]))>

As we’ve seen, the obs we defined is the space of Gauss: this acts as the default limits whenever needed (e.g. for sampling). gauss also has a norm_range, which equals by default as well to the obs given, however, we can explicitly change that with set_norm_range.

We can also access the parameters of the PDF in two ways, depending on our intention: either by name (the parameterization name, e.g. mu and sigma, as defined in the Gauss), which is useful if we are interested in the parameter that describes the shape

gauss.params

OrderedDict([('mu', <zfit.Parameter 'mu' floating=True value=1>),
             ('sigma', <zfit.Parameter 'sigma42' floating=True value=1>)])

or to retrieve all the parameters that the PDF depends on. As this now may sounds trivial, we will see later that models can depend on other models (e.g. sums) and parameters on other parameters. There is one function that automatically retrieves all dependencies, get_params. It takes three arguments to filter:

floating: whether to filter only floating parameters, only non-floating or don’t discriminate
is_yield: if it is a yield, or not a yield, or both
extract_independent: whether to recursively collect all parameters. This, and the explanation for why independent, can be found later on in the Simultaneous tutorial.

Usually, the default is exactly what we want if we look for all free parameters that this PDF depends on.

gauss.get_params()

OrderedSet([<zfit.Parameter 'mu' floating=True value=1>, <zfit.Parameter 'sigma42' floating=True value=1>])

The difference will also be clear if we e.g. use the same parameter twice:

gauss_only_mu = zfit.pdf.Gauss(obs=obs, mu=mu, sigma=mu)
print(f"params={gauss_only_mu.params}")
print(f"get_params={gauss_only_mu.get_params()}")

params=OrderedDict([('mu', <zfit.Parameter 'mu' floating=True value=1>), ('sigma', <zfit.Parameter 'mu' floating=True value=1>)])

get_params=OrderedSet([<zfit.Parameter 'mu' floating=True value=1>])

Functionality¶

PDFs provide a few useful methods. The main features of a zfit PDF are:

pdf: the normalized value of the PDF. It takes an argument norm_range that can be set to False, in which case we retrieve the unnormalized value
integrate: given a certain range, the PDF is integrated. As pdf, it takes a norm_range argument that integrates over the unnormalized pdf if set to False
sample: samples from the pdf and returns a Data object

integral = gauss.integrate(limits=(-1, 3))  # corresponds to 2 sigma integral
integral

<tf.Tensor: shape=(1,), dtype=float64, numpy=array([0.95449974])>

Tensors¶

As we see, many zfit functions return Tensors. This is however no magical thing! If we’re outside of models, than we can always safely convert them to a numpy array by calling zfit.run(...) on it (or any structure containing potentially multiple Tensors). However, this may not even be required often! They can be added just like numpy arrays and interact well with Python and Numpy:

np.sqrt(integral)

array([0.97698502])

They also have shapes, dtypes, can be slices etc. So do not convert them except you need it. More on this can be seen in the talk later on about zfit and TensorFlow 2.0.

sample = gauss.sample(n=1000)  # default space taken as limits
sample

<zfit.core.data.SampleData at 0x7f84b4f2e8b0>

sample.unstack_x()[:10]

<tf.Tensor: shape=(10,), dtype=float64, numpy=
array([ 1.2152566 ,  2.81183774,  0.3578351 ,  0.14692911, -0.46053768,
        1.0086027 ,  0.24786822,  0.79486641, -0.42387978,  0.31100973])>

sample.n_obs

sample.obs

('obs1',)

We see that sample returns also a zfit Data object with the same space as it was sampled in. This can directly be used e.g.

probs = gauss.pdf(sample)
probs[:10]

<tf.Tensor: shape=(10,), dtype=float64, numpy=
array([0.38980595, 0.07728027, 0.32461143, 0.27725891, 0.13730869,
       0.39892752, 0.30065567, 0.39063626, 0.14476329, 0.31465064])>

NOTE: In case you want to do this repeatedly (e.g. for toy studies), there is a way more efficient way (see later on)

Plotting¶

so far, we have a dataset and a PDF. Before we go for fitting, we can make a plot. This functionality is not directly provided in zfit (but can be added to zfit-physics). It is however simple enough to do it:

def plot_model(model, data, scale=1, plot_data=True):  # we will use scale later on

    nbins = 50

    lower, upper = data.data_range.limit1d
    x = tf.linspace(lower, upper, num=1000)  # np.linspace also works
    y = model.pdf(x) * size_normal / nbins * data.data_range.area()
    y *= scale
    plt.plot(x, y)
    data_plot = zfit.run(z.unstack_x(data))  # we could also use the `to_pandas` method
    if plot_data:
        plt.hist(data_plot, bins=nbins)

plot_model(gauss, data_normal)

We can of course do better (and will see that later on, continuously improve the plots), but this is quite simple and gives us the full power of matplotlib.

Different models¶

zfit offers a selection of predefined models (and extends with models from zfit-physics that contain physics specific models such as ARGUS shaped models).

print(zfit.pdf.__all__)

['BasePDF', 'BaseFunctor', 'Exponential', 'CrystalBall', 'DoubleCB', 'Gauss', 'Uniform', 'TruncatedGauss', 'WrapDistribution', 'Cauchy', 'Poisson', 'Chebyshev', 'Legendre', 'Chebyshev2', 'Hermite', 'Laguerre', 'RecursivePolynomial', 'ProductPDF', 'SumPDF', 'GaussianKDE1DimV1', 'FFTConvPDFV1', 'ConditionalPDFV1', 'ZPDF', 'SimplePDF', 'SimpleFunctorPDF']

To create a more realistic model, we can build some components for a mass fit with a

signal component: CrystalBall
combinatorial background: Exponential
partial reconstructed background on the left: Kernel Density Estimation

mass_obs = zfit.Space('mass', (0, 1000))

# Signal component

mu_sig = zfit.Parameter('mu_sig', 400, 100, 600)
sigma_sig = zfit.Parameter('sigma_sig', 50, 1, 100)
alpha_sig = zfit.Parameter('alpha_sig', 200, 100, 400, floating=False)  # won't be used in the fit
n_sig = zfit.Parameter('n sig', 4, 0.1, 30, floating=False)
signal = zfit.pdf.CrystalBall(obs=mass_obs, mu=mu_sig, sigma=sigma_sig, alpha=alpha_sig, n=n_sig)

# combinatorial background

lam = zfit.Parameter('lambda', -0.01, -0.05, -0.001)
comb_bkg = zfit.pdf.Exponential(lam, obs=mass_obs)

part_reco_data = np.random.normal(loc=200, scale=150, size=700)
part_reco_data = zfit.Data.from_numpy(obs=mass_obs, array=part_reco_data)  # we don't need to do this but now we're sure it's inside the limits

part_reco = zfit.pdf.GaussianKDE1DimV1(obs=mass_obs, data=part_reco_data, bandwidth='adaptive')

Composing models¶

We can also compose multiple models together. Here we’ll stick to one dimensional models, the extension to multiple dimensions is explained in the “custom models tutorial”.

Here we will use a SumPDF. This takes pdfs and fractions. If we provide n pdfs and:

n - 1 fracs: the nth fraction will be 1 - sum(fracs)
n fracs: no normalization attempt is done by SumPDf. If the fracs are not implicitly normalized, this can lead to bad fitting behavior if there is a degree of freedom too much

sig_frac = zfit.Parameter('sig_frac', 0.3, 0, 1)
comb_bkg_frac = zfit.Parameter('comb_bkg_frac', 0.25, 0, 1)
model = zfit.pdf.SumPDF([signal, comb_bkg, part_reco], [sig_frac, comb_bkg_frac])

In order to have a corresponding data sample, we can just create one. Since we want to fit to this dataset later on, we will create it with slightly different values. Therefore, we can use the ability of a parameter to be set temporarily to a certain value with

print(f"before: {sig_frac}")
with sig_frac.set_value(0.25):
    print(f"new value: {sig_frac}")
print(f"after 'with': {sig_frac}")

before: <zfit.Parameter 'sig_frac' floating=True value=0.3>

new value: <zfit.Parameter 'sig_frac' floating=True value=0.25>

after 'with': <zfit.Parameter 'sig_frac' floating=True value=0.3>

While this is useful, it does not fully scale up. We can use the zfit.param.set_values helper therefore. (Sidenote: instead of a list of values, we can also use a FitResult, the given parameters then take the value from the result)

with zfit.param.set_values([mu_sig, sigma_sig, sig_frac, comb_bkg_frac, lam], [370, 34, 0.18, 0.15, -0.006]):
    data = model.sample(n=10000)

plot_model(model, data);

Plotting the components is not difficult now: we can either just plot the pdfs separately (as we still can access them) or in a generalized manner by accessing the pdfs attribute:

def plot_comp_model(model, data):
    for mod, frac in zip(model.pdfs, model.params.values()):
        plot_model(mod, data, scale=frac, plot_data=False)
    plot_model(model, data)

plot_comp_model(model, data)

Now we can add legends etc. Btw, did you notice that actually, the frac params are zfit Parameters? But we just used them as if they were Python scalars and it works.

print(model.params)

OrderedDict([('frac_0', <zfit.Parameter 'sig_frac' floating=True value=0.3>), ('frac_1', <zfit.Parameter 'comb_bkg_frac' floating=True value=0.25>), ('frac_2', <zfit.ComposedParameter 'Composed_autoparam_2' params=OrderedDict([('param_0', <zfit.Parameter 'sig_frac' floating=True value=0.3>), ('param_1', <zfit.Parameter 'comb_bkg_frac' floating=True value=0.25>)]) value=0.45>)])

Extended PDFs¶

So far, we have only looked at normalized PDFs that do contain information about the shape but not about the absolute scale. We can make a PDF extended by adding a yield to it.

The behavior of the new, extended PDF does NOT change, any methods we called before will act the same. Only exception, some may require an argument less now. All the methods we used so far will return the same values. What changes is that the flag model.is_extended now returns True. Furthermore, we have now a few more methods that we can use which would have raised an error before:

get_yield: return the yield parameter (notice that the yield is not added to the shape parameters params)
ext_{pdf,integrate}: these methods return the same as the versions used before, however, multiplied by the yield
sample is still the same, but does not require the argument n anymore. By default, this will now equal to a poissonian sampled n around the yield.

The SumPDF now does not strictly need fracs anymore: if all input PDFs are extended, the sum will be as well and use the (normalized) yields as fracs

The preferred way to create an extended PDf is to use PDF.create_extended(yield). However, since this relies on copying the PDF (which may does not work for different reasons), there is also a set_yield(yield) method that sets the yield in-place. This won’t lead to ambiguities, as everything is supposed to work the same.

yield_model = zfit.Parameter('yield_model', 10000, 0, 20000, step_size=10)
model_ext = model.create_extended(yield_model)

alternatively, we can create the models as extended and sum them up

sig_yield = zfit.Parameter('sig_yield', 2000, 0, 10000, step_size=1)
sig_ext = signal.create_extended(sig_yield)

comb_bkg_yield = zfit.Parameter('comb_bkg_yield', 6000, 0, 10000, step_size=1)
comb_bkg_ext = comb_bkg.create_extended(comb_bkg_yield)

part_reco_yield = zfit.Parameter('part_reco_yield', 2000, 0, 10000, step_size=1)
part_reco.set_yield(part_reco_yield)  # unfortunately, `create_extended` does not work here. But no problem, it won't change anyting.
part_reco_ext = part_reco

model_ext_sum = zfit.pdf.SumPDF([sig_ext, comb_bkg_ext, part_reco_ext])

zfit

Introduction to zfit¶

Data¶

Model¶

Parameter¶

Functionality¶

Tensors¶

Plotting¶

Different models¶

Composing models¶

Extended PDFs¶

Loss¶

Extended loss¶

Minimization¶

Pre-fit parts of the PDF¶

Fit result¶

Estimating uncertainties¶

What is ‘new_result’?¶

A simple profile¶

End of zfit¶

Beginning of hepstats¶