\subsection{Model}
\label{sec:model}
Building models is the core competence of \zfit{}. As seen in the
example in Sec.
\ref{sec:quickstart}, this can be done in a simple manner by using
already implemented models and possibly combining them, but the models can also
be completely custom built. Within \zfit{}, there are two basic types of models
to cover most cases:
Functions and Probability Density Functions (PDF).
The basic features of a model include
\begin{itemize}
\item Each model is defined inside a \zspace{}. Its dimension are
``observables", simple string identifiers as previously discussed.
\item A model implements a function that returns its value
depending on some data. This is either \pyth{pdf} for a PDF or \pyth{func}
for a Func.
\item Full as well as partial integration over a model is
possible. This includes numerical as well as analytical integration, if
available.
\item Generating a sample following the models shape using
numerical or analytical methods, the latter only if available.
\end{itemize}
The value function as well as the integration and sampling are implemented to
return pure
Tensors. Depending on the task, higher-level methods providing either a more
intuitive, imperative behaviour or a significantly more
performant execution for certain cases such as repetitive pseudo-experiments
are also implemented.
The main differences between a PDF and a Function concern the normalisation
and the output dimensionality. This leads to a few subtle differences.
\begin{description}
\item[PDF] A PDF $f(x)$ is only well-defined with a given normalisation
range.
This defines the normalisation constant so that, with the limits from
$lower$ to $upper$, the integral over the PDF equals to one as
\begin{equation}
\label{eq:norm_range}
\sum_{i}\int_{l_i}^{u_i} f(x) dx = 1
\end{equation}
with $i$ indexing all the limits that make up the domain, $l_i$ and $u_i$
are the lower respectively upper limit.
A PDF object has three special attributes, which are
\begin{itemize}
\item
the probability density function \pyth{pdf}. It returns always a rank
one Tensor,
i.e. a simple vector, with the length number of data points.
\item the probability density function
\textit{without} the normalisation constant, \pyth{unnormalised_pdf}.
In cases where only the
shape is
needed, using this function is less expensive, especially
if the normalisation has to be computed numerically.
\item the limits that define the normalisation constant of Eq.
\ref{eq:norm_range}. They can be set using the
\pyth{set_norm_range}
method.
\end{itemize}
\item[Function] A Function is in a way more simple and general than a PDF.
It takes the same values as a PDF but returns something with dimensionality
$\mathbb{R}^m$. It can be used to transform values or to use it as a
building block for more complex expressions.
It has the method \pyth{func} to evaluate its value and no other special
attributes.
\end{description}
\subsubsection{Parametrization}
\label{sec:parametrization}
Models can be parametrized by \pyth{Parameter} objects which can be used in the
implementation of the shape function like
any other
Tensor-like object when building the model. In
the following we will first have a look at the \pyth{Parameter} itself.
Afterwards we will see how they are exactly used with a model.
There is a distinction between dependent and independent parameters as only the
latter can be changed directly and have limits while the former are any
arbitrary combination of them.
An independent \pyth{Parameter}
\begin{itemize}
\item has a name with purely descriptive purpose;
\item has an initial value;
\item maybe has lower and upper limits;
\item is either floating or not, independent of whether limits were
specified;
\item has a step size indicating the order of magnitude of the parameter
which can be given. Otherwise, it is automatically inferred. A
well chosen step size improves the minimization process and can be
critical, mostly in the absence of limits and with a weak dependence on the
model, so that large steps will be required to change the model.
\end{itemize}
Currently, the shape of parameters is implicitly restricted to a
scalar.\footnote{This
comes from the fact that the PDF will have different sized data as input which
is not controlled by the user, such as when doing numerical integration. Any
parameter therefore has to be able to broadcast seamlessly.} As a consequence,
a
parameter cannot simply be treated like data as it is possible in \roofit. A
function \pyth{Parameter} which depends on the data itself will most likely be
available in the future.
The name of a parameter, as any other name for a single object in \zfit{}, has
purely descriptive character. There is purposely\footnote{Contrary to the
\root framework. If the need ever arises, adding this as an additional feature
is relatively simple.} no direct way provided to access parameters by name:
instead of using names the actual parameter object is passed around. This
has the advantage of avoiding double bookkeeping the parameters, since then a
reference on an object as well as on a string would be needed. The name is
mandatory for parameters, as opposed to other objects in \zfit{}, since
matching the value of a parameter to the right name is a critical task during
and after a minimization, when reading off the right value.
Composed parameters are dependent parameters. In general they can be any
Tensor, that is a result of any kind of operation. They can depend on zero, one
or
more independent parameters, as composed parameters are arbitrary functions;
therefore, operations such as shifting and scaling are included as a trivial
subset. For more details on the
actual implementation and the dependency management see Appendix
\ref{appendix:dependency management}.
Composed parameters can neither be floating nor have limits currently, since
the independent parameters a composed parameter may depends on can have
arbitrary relations and have to be restricted themselves. In order to restrict
a parameter, arbitrary constraints can be given to the loss instead.
Every model can depend on multiple parameters, both dependent and independent.
Each parameter that parametrises the model has a name
specific to the model and is given on instantiation. For example a
\texttt{Gauss} has parameters named \texttt{mu} and \texttt{sigma} as in the
example in Sec. \ref{sec:quickstart}. They are stored in a mapping attribute
named \pyth{params}.
\begin{center}
\begin{minipage}{\textwidth}
\begin{python}
mu_param = gauss.params["mu"]
sigma_param = gauss.params["sigma"]
\end{python}
\end{minipage}
\end{center}
Notice that this is not contradictory to the statement above that single
objects \textit{cannot} be accessed by \textit{their} name. Each object has a
unique identifier, the name, but objects can have names for their constituents
that are \textit{not} unique, like \texttt{mu}. This simply describes a part of
the PDF and any \texttt{Gauss} will have a parameter \texttt{mu}.
\subsubsection{Implementing a custom PDF}
An essential feature of \zfit{} is the ability to simply create custom models.
There is a large freedom in building models from the \pyth{BaseModel} class,
since
it takes care of most boilerplate and has well defined entry points than can be
customized. The full implementation details and possibilities for
customizations are described in Appendix \ref{appendix:basemodel}.
For the most
common use cases though, there exists a simple way of creating a custom model.
The \pyth{ZPDF}, basically a more user
friendly wrapper around \pyth{BasePDF}, can be be used as a base class in these
simple cases.
The following function will be used as the PDF shape
\begin{equation}
\label{eq:simple_linear_model}
f(x, y) = a \cdot x^2 + b \cdot y^4.
\end{equation}
To implement this function,
\pyth{_unnormalized_pdf} has to be
overwritten. For the vast majority of custom models, this is the only method
to be overwritten. Changing other methods, especially \pyth{_pdf}, is an
advanced feature and only needed in special cases. For more details on the
customization and the possibility of hooking into the calls see Appendix
\ref{appendix:basemodel}.
This is a two dimensional PDF with two parameters. To implement
it in \zfit{}, a new class is created
\begin{center}
\begin{minipage}{\textwidth}
\label{py:custom_pdf_linear_model}
\begin{python}
class EvenPolyPDF(zfit.pdf.ZPDF):
"""Implementation of f(x, y) = a*x^2 + b*y^4"""
_PARAMS = ['a', 'b']
_N_OBS = 2 # since two dimensional
def _unnormalized_pdf(x):
xdata, ydata = x.unstack_x()
a = self.params['a']
b = self.params['b']
return a * xdata ** 2 + b * ydata ** 4
\end{python}
\end{minipage}
\end{center}
Note that we only need the \textit{shape} of the function but do not need to
take care of the normalisation, as numerical methods are already implemented in
the base class.
We can see the advantage of the preprocessing that is done by the base class,
especially the reordering of the data \pyth{x}. It is a \zdata{} object and
calling the \pyth{unstack_x} method returns a list\footnote{Or a single Tensor
fo the one dimensional case if not specified differently in the arguments.}
containing a Tensor for each column sorted
according to the models observables. The length of the list has to
coincide with the specified \pyth{_N_OBS}, the number of observables. It is not
mandatory to specify this field in a Model and is sometimes not possible to
know previously, in
which case it can simply be left away.
The naming of the parametrization of the function is defined with
\pyth{_PARAMS}. These exact names have to be used when creating an instance of
the model. The parameters are then stored in the \pyth{params} dictionary and
extracting them
is usually the first step inside \pyth{_unnormalized_pdf}.
Documentation plays an important role here: it defines the name of the
parameters that have to be used and what they represent in the function, but it
is also crucial to communicate the ordering of the data. In this case, a user
can see that the first observable and the corresponding column in the
data will be used as \textit{x} in the function.
This PDF is already complete and works out of the box. We can create an
instance and use its methods. As an example, the observables of the instance
will be called "xobs" and "yobs" with the normalisation range going from $0$ to
$10$ and from $0$ to $5$, respectively.
\begin{center}
\begin{minipage}{\textwidth}
\label{py:instance_custom_pdf_linear_model}
\begin{python}
param_a = zfit.Parameter("a", 1.)
param_b = zfit.Parameter("b", 2.)
x_obs = zfit.Space("xobs", limits=(0, 10))
y_obs = zfit.Space("yobs", limits=(0, 5))
obs = x_obs * y_obs
poly_model = EvenPolyPDF(a=param_a, b=param_b, obs=obs)
prob = poly_model.pdf(...) # some data needed
x_limits = zfit.Space("xobs", limits=(3, 5))
y_limits = zfit.Space("yobs", limits=(1, 3))
integral_limits = x_limits * y_limits
integral = poly_model.integrate(integral_limits)
sample = poly_model.sample(n=1000)
\end{python}
\end{minipage}
\end{center}
Using the \pyth{integrate} method, we obtain the integral $i$ of our normalized
PDF $f_{norm}$
\begin{align}
\label{eq:model integral}
i &= \int_3^5 \mathrm{d}x \int_{1}^{3} \mathrm{d}y ~f_{normed} (x, y) \\
&\stackrel{(\ref{eq:pdf from func})}{=} \frac{\int_3^5 \mathrm{d}x
\int_{1}^{3}
\mathrm{d}y ~f(x,
y)}{\int_0^{10}
\mathrm{d}x \int_{0}^{5} \mathrm{d}y ~f(x, y)}
\end{align}
with $x$ and $y$ corresponding to the observables \pyth{x} and \pyth{y},
respectively, and $f(x, y)$ is the unnormalised PDF as defined in Eq.
\ref{eq:simple_linear_model}. Using Eq. \ref{eq:pdf from func} we end up with a
precise formulation of what is
\textit{actually} executed in \zfit{}.
In other words, the integral was
calculated over the limits $3$ to $5$ and $1$
to $3$
respectively and normalised over the normalisation range \texttt{obs}. The
latter is defined on initialisation and taken as the default
normalisation range. As mentioned before, the limits could
also be multiplied in a different order resulting in \texttt{limits} having the
order (``yobs", ``xobs"). The model takes care internally (see also Appendix
\ref{sec:public_methods}) that the right limits are at the right place.
The \pyth{sample} method is used to generate 1000 points from the model.
With the instance created, also probability densities \pyth{prob} for points
can
be calculated by calling the method with some \zdata{} object called
\pyth{data} as follows
\begin{center}
\begin{minipage}{\textwidth}
\begin{python}
probs = poly_model.pdf(data)
\end{python}
\end{minipage}
\end{center}
It is worth noting that none of
the
operations are executed yet and what is returned by the methods are
\text{Tensors}, as described in Sec. \ref{sec:implementation}. To run the
actual
computations, it is necessary to call \pyth{zfit.run(...)} on any Tensor.
Running
\begin{center}
\begin{minipage}{\textwidth}
\begin{python}
probs_np = zfit.run(probs_np)
integral_np = zfit.run(integral)
\end{python}
\end{minipage}
\end{center}
returns an actual numpy array for \pyth{probs_np} and a Python float for
\pyth{integral_np}.
For further, advanced customization of the PDF, methods can be overridden as
described in Appendix \ref{sec:public_methods}. However, in most use cases
there exist better ways: for example, to add an analytic
integral to the model, overwriting \pyth{_analytic_integrate} should be the
last resort. Instead, integrals
defined over a specific range only or the whole range and over all dimensions
or just partially can be registered with a model. A priority attribute allows
to specify preferences on one
over other methods.
A typical use case for this feature is a special integral
that is known
exactly, for example the value of the integral over the full space of a
Gaussian shaped model\footnote{As already hinted in Eq. \ref{eq:model
integral}, it is important to
note that the integral \textit{to be implemented} is over the
unnormalised PDF as implemented in the \pyth{_unnormalized_pdf} method.} An
integral over both dimensions is registered but partial integrals
could be added as well
\begin{equation}
\label{integral_xy}
\int_{x_0}^{x_1} \int_{y_0}^{y_1} f(x, y) \mathrm{d}x \mathrm{d}y = (1 / 3
\cdot a
\cdot x^3 + 1 / 5 \cdot b \cdot y^5) \Big|_{x_0}^{x_1} \Big|_{y_0}^{y_1}.
\end{equation}
In case a full integral is requested but no analytic
integral
over all dimensions is available, the fallback of \pyth{integrate} looks for
partial integrals. If available, it uses them to numerically integrate over the
remaining dimensions instead of using the unnormalised PDF.
The implementation of the integral with \zfit{} is done by registering it to
the PDF
\begin{center}
\begin{minipage}{\textwidth}
\begin{python}
def integral_full(x, limits, norm_range, params, model):
lower, upper = limits.limit1d
a = params['a']
b = params['b']
lower = ztf.convert_to_tensor(lower)
upper = ztf.convert_to_tensor(upper)
def indef_integral(limit):
return 1 / 3 * a ** 3 + 1 / 5 * b * y ** 5
return indef_integral(upper) - indef_integral(lower)
lower = ((zfit.Space.ANY_LOWER, zfit.Space.ANY_LOWER),)
upper = ((zfit.Space.ANY_UPPER, zfit.Space.ANY_UPPER),)
limits = zfit.Space.from_axes(axes=(0, 1), limits=(lower, upper))
EvenPolyPDF.register_analytic_integral(func=integral_full,
limits=limits)
\end{python}
\end{minipage}
\end{center}
Since the observables that will be assigned to each axis are unknown, the
\zspace{} has to be defined with axes, not with observables. This follows the
principle of using axes inside the model as explained in \ref{sec:spaces and
dimensions}. Instead of \pyth{ANY_LOWER} and \pyth{ANY_UPPER}, it could also be
defined over a specific domain
by using plain Python floats. The PDF will now automatically use this integral
if possible. Otherwise, as for example when creating a partial integral, it
will silently fall back to numerical integration.
This example demonstrates how to implement a shape that does only depend on the
data and parameters. Models, that also depend on other models, such as the
\pyth{SumPDF} from Sec. \ref{sec:quickstart}, require to be built from a
\pyth{Functor}. This implies as a minor change only an additional argument to
the base
class, the depending models, in order to track the dependencies
correctly. This is described in more details in
Appendix \ref{appendix:functor}.
\subsubsection{Sampling}
Sampling from a model can be done in two different ways to cover two distinct
use cases.
On one hand, the method \pyth{sample} returns a pure Tensor, which behaves as
any other
sampling algorithm in TF, for example \pyth{tf.random.uniform}. This can be
useful for specific and
mostly advanced cases, but he overall behaviour can be unintuitive for
inexperienced users, as discussed in Appendix \ref{appendix:tensorflow}.
On the other hand, sampling from a
model is typically used to perform toy studies, which consists of the
generation
of events
according to the model and afterwards a fit of the model with randomly
initialized parameters to the sampled data. With \pyth{create_sampler}, a
\zdata -like object which handles
the
correct storage and the sampling from the model is created. As any other
\zdata{}, it can
be used to build a loss. The actual data is sampled by invoking \pyth{resample}
and stays unchanged until it is invoked again. This sampler keeps a dependency
on the original model and uses it to sample, keeping the original
parameter values it was created with. If desired, any parameter that has
changed in between will be at that new value when \pyth{resample} is called,
effectively changing the sampling shape.
\footnote{Just to be clear: the default behaviour is that the sampler
samples from the exact distribution that it was created with \textit{in terms
of parameter values}.}
To illustrate the use of the toy sampler, let's look at a toy study example.
Note that the actual (large)
\textit{objects} needed
such as the model, minimiser, loss etc. are created \textit{outside} of the
loop and only the necessary \textit{methods} are called \textit{inside}. In
general, this removes boilerplate code from the additional object creation but
is especially important when working with a graph based backend.
We assume to have previously built a model already and a minimizer as for
example shown in Sec. \ref{sec:quickstart}, and we generate 1000 events
according to our model.
\begin{center}
\begin{minipage}{\textwidth}
\begin{python}
sampler = model.create_sampler(n=1000)
nll = zfit.loss.UnbinnedNLL(model, sampler)
for _ in range(n_toys):
sampler.resample() # now the sampling happens
for param in model.get_dependents(only_floating=True):
param.set_value(np.random.normal()) # any initial value
result = minimizer.minimize(nll)
if result.converged: # check if the fit was successful
... # safe results in a list or similar
\end{python}
\end{minipage}
\end{center}
The sampling is implemented with the accept-reject method that works with
arbitrary shapes. If an analytic
inverse integral is registered, this will be used, providing a more efficient
way of sampling. More details about the different techniques and the
use of importance sampling are described in Appendix
\ref{appendix:sampling techniques}.
\subsubsection{Extended PDFs}
In addition to the shape, a PDF can carry more information, namely a yield, in
which case we refer to it as an Extended PDF.
The yield is a scale that describes the magnitude of the PDF, typically
reflecting
the number of events in the data sample. It can be used to multiply the output
of \pyth{pdf}, the normal probability density function, which results in a
number probability. This implies that when multiplying the integral with the
yield, the number of events is retrieved instead of the probability over a
certain range.\footnote{\textit{Currently}, the integral of an
extended PDF is multiplied by the yield \textit{by default} but this behaviour
is meant to change in a future version}.
To count for example how many signal particles are in
the sample, a composite PDF existing of a background and a signal component,
both extended,
can be fitted to the data. The \pyth{SumPDF} used to build this composition
does not need a fraction if the PDFs are already extended, but uses the yields,
normalised to one, as fractions. In this case, integrating the PDF that
represents the signal
shape over the whole range returns the number of signal events.
To create an extended PDF, \pyth{create_extended} can be used. This method
returns a copy of the PDF and adds a yield to it.
