# Tree boosting for learning EFT parameters

Suman Chatterjee  Nikolaus Frohner  Lukas Lechner  Robert Schöfbeck  Dennis
Schwarz  Institute of High Energy Physics (HEPHY), Austrian Academy of
Sciences (ÖAW), Nikolsdorfer Gasse 18, 1050 Vienna, Austria  TU Wien,
Karlsplatz 13, 1040 Vienna, Austria

###### Abstract

We present a new tree boosting algorithm designed for the measurement of
parameters in the context of effective field theory (EFT). To construct the
algorithm, we interpret the optimized loss function of a traditional decision
tree as the maximal Fisher information in Poisson counting experiments. We
promote the interpretation to general EFT predictions and develop a suitable
boosting method. The resulting “Boosted Information Tree” algorithm
approximates the score, the derivative of the log-likelihood function with
respect to the parameter. It thus provides a sufficient statistic in the
vicinity of a reference point in parameter space where the estimator is
trained. The training exploits per-event information of likelihood ratios for
different theory parameter values available in the simulated EFT data sets.

###### keywords:

LHC; physics beyond the standard model; machine learning; effective field
theory; boosted decision trees; Fisher information

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: BIT (Boosted Information Trees)
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/HephyAnalysisSW/BIT
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: GPLv3
Programming language: Python2 and Python3
Nature of problem: Providing a discriminator for parameter estimation in the
context of the standard model effective field theory.
Solution method: A tree-based algorithm exploits “augmented” information of
the simulated training data set to regress in the score function and thereby
provides a sufficient test statistic of an EFT parameter.

##  1 Introduction

An important front at the Large Hadron Collider (LHC) experiments is the
search for indirect effects in standard model (SM) measurements induced by
phenomena beyond the SM (BSM). If the BSM energy scale is much higher than the
reach of the LHC, indirect non-resonant effects at lower energies may still be
discernible via dispersed and subtle differences in the spectra of kinematic
observables. These are accessible by performing precision measurements of
differential cross sections. The SM effective field theory (SM-EFT) [1, 2, 3,
4, 5, 6] parameterizes such deviations and has become the leading generative
model in the first attempts of a global interpretation of collider results [7,
8, 9, 10, 11, 12]. The SM-EFT has 2499 degrees of freedom, among which 59 are
flavor diagonal and respect baryon and lepton number conservation [6, 13].

Particle physics event generation, both, for SM-EFT and in general, consists
of consecutive steps involving simulation tools for the hard interaction, the
parton shower, the matrix element (ME) matching, the hadronization of strongly
interacting partons, and the detector response to the stable particles. The
authors of Refs. [14, 15, 16, 17, 18] employ this structure to show how mean
squared error (MSE) loss functions can regress on the SM-EFT likelihood, the
likelihood ratio between different SM-EFT hypotheses, or the gradient in SM-
EFT parameter space of the likelihood function, i.e., the score vector. The
key element is additional information from Monte-Carlo simulations that
characterizes the dependence of the likelihood function on the SM-EFT
parameters. This augmented data is used in a number of sophisticated machine
learning methods based on neural networks.

Tree boosting [19, 20, 21, 22], on the other hand, has received much less
attention in the context of probing EFT operators. Boosting algorithms provide
a strong learner by iteratively training an ensemble of weak learners to the
pseudo-residuals of the previous iteration step. Traditional Boosted Decision
Trees (BDTs) for classification tasks are trained with discretely categorized
data using, e.g., samples of simulated signal and background events and
predict the target label. The BDTs exhibit robust performance in a wide range
of classification problems in high-energy physics (HEP). In the context of
EFT, tree-based discriminators have not been a focus of dedicated development.
A reason may lie in the absence, so far, of a simple strategy to exploit the
augmented data in the training.

In this work, we change that and develop a tree-based algorithm trained on a
simulated event sample generated at a single EFT parameter point. We work with
a single scalar parameter, possibly after a linear basis transformation, such
that the hypotheses correspond to the two scalar values $\theta_{\textrm{0}}$
and $\theta$. The resulting algorithm is called “Boosted Information Tree”
(BIT), consisting of a weak learner and a boosting algorithm. Its loss
function uses the augmented data to regresses to the score function. The
augmented data are per-event weight functions with a simple analytic
dependence on the EFT theory parameters. By approximating the score function,
the BIT discriminator provides a locally sufficient statistic of the
parameter.

The structure of the article is as follows. In Sec. 2, we give an overview of
the relevant previous work. A detailed description of the structure of the
augmented data and the training data set is provided in Sec. 3. In Sec. 4, we
interpret the fit of a traditional decision tree, the weak learner of a BDT,
in terms of an optimization of the Fisher information in Poisson counting
experiments. This interpretation serves as a starting point for constructing
the algorithm in Sec. 5. The approximation of the score function is validated
in analytically tractable toy models, discussed in Sec. 6. Conclusions and
outlook are provided in Sec. 7. Public code implementing the algorithm is
available on GitHub [23], and a brief introduction to the interface is
provided in A.

##  2 Related work

Below we briefly mention the algorithms developed in the recent past and are
closely related to the method we report in this article.

###  2.1 SALLY/SALLINO

In a series of publications, the authors of Refs. [15, 16, 17, 18] introduce
algorithms that use augmented data from the event simulation to regress to the
true likelihood or the true score even if the quantities are intractable. It
is shown how intractable factors in the likelihood, related to the parton
shower, the ME matching, the detector simulation, etc., cancel and the
augmented data is subsequently leveraged to define a number of sophisticated
inference techniques based on neural networks. This methodology is also a
foundation of the present work. At variance with the neural network
approximation, we focus on the training of tree-based discriminators. The BIT
operates in the same setting as the SALLY and SALLINO algorithms.

###  2.2  INFERNO

In Ref. [24], a non-linear summary statistic is constructed by minimizing
inference-motivated losses via stochastic gradient descent. Similar to the
present work, the algorithm uses the Fisher information of the parameters of
interest. In contrast to our work, a neural network is trained, and the
INFERNO training loss is chosen to minimize the expected variance in the
context of parameter regression.

###  2.3  Binary classifier metrics and weight derivative regression

In Refs. [25, 26], various measures for classification and regression
performance that involve the Fisher information in Poisson counting
experiments are introduced. The interpretation of a decision tree with the
Gini impurity as a Fisher information-optimal cross section measurement is the
starting point of the development of the BIT. We do not consider the various
ratios of the Fisher information at the different stages of an estimation
procedure. In the present work, we merely follow up on the proposal to
generalize node split criteria towards loss functions involving the Fisher
information and present a concrete boosted tree algorithm.

##  3 Structure of the training data

At the level of the hard scattering, event simulation of squared MEs provides
a parton-level configuration ${\bm{z}}\sim p({\bm{z}}|{\bm{\theta}})$, where
the vector ${\bm{\theta}}$ denotes the EFT parameters and
$p({\bm{z}}|{\bm{\theta}})$ the parton-level likelihood. The generic structure
of a matrix element in terms of Feynman amplitudes with a single EFT operator
insertion leads to predicted differential cross sections of the general
form111EFT operators that mix with the SM kinetic terms can lead, in
principle, to non-polynomial dependence on $\bm{\theta}$. This behavior does
not spoil the following discussion, and there are SM-EFT models that
conveniently support the truncation of kinetic mixing effects and propagator
corrections to a configurable polynomial order, e.g., Ref. [27].

$\displaystyle\textrm{d}\sigma({\bm{\theta}})$ $\displaystyle\propto$
$\displaystyle|\mathcal{M}_{\textrm{SM}}({\bm{z}})+\theta_{a}\mathcal{M}^{a}_{\textrm{BSM}}({\bm{z}})|^{2}\textrm{d}{\bm{z}}.$
(1)

The leading BSM contribution is
$2\theta_{a}\textrm{Re}\left(\mathcal{M}^{\ast}({\bm{z}})_{\textrm{SM}}\mathcal{M}^{a}_{\textrm{BSM}}({\bm{z}})\right)$
and describes the interference of the BSM and SM amplitudes. It is the only
term in the expansion in $\bm{\theta}$ that does not also receive
contributions from EFT operators with higher mass dimensions. We write the
likelihood $p(\bm{z}|\bm{\theta})$, the likelihood ratio
$r(\bm{z}|\bm{\theta},\bm{\theta}_{0})$, and the score vector
$\bm{t}(\bm{z}|\bm{\theta})$ as

$\displaystyle p(\bm{z}|\bm{\theta})$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma(\bm{\theta})}\frac{\textrm{d}\sigma({\bm{\theta}})}{\textrm{d}\bm{z}},$
(2) $\displaystyle r(\bm{z}|\bm{\theta},\bm{\theta}_{0})$ $\displaystyle=$
$\displaystyle\frac{\textrm{d}\sigma({\bm{\theta}})/\textrm{d}\bm{z}}{\textrm{d}\sigma({\bm{\theta}_{\textrm{0}}})/\textrm{d}\bm{z}}\cdot\frac{\sigma(\bm{\theta}_{\textrm{0}})}{\sigma(\bm{\theta})},\;\textrm{and}$
(3) $\displaystyle\bm{t}(\bm{z}|\bm{\theta})$ $\displaystyle=$
$\displaystyle\bm{\nabla}_{\theta}\log p(\bm{z}|\bm{\theta}),$ (4)

where $\sigma(\bm{\theta})$ and
$\textrm{d}\sigma({\bm{\theta}})/\textrm{d}\bm{z}$ are the total and the
differential parton-level cross sections, respectively.

The simple quadratic structure in Eq. 1 has important practical implications:
If a sample of $N$ events $\\{\bm{z}_{i}\\}_{i=1}^{N}$ is generated from the
probability density function (PDF) in Eq. 2 at a reference parameter point
$\bm{\theta}_{\textrm{ref}}$, the value of the matrix element can be
reevaluated at linearly independent basis parameter points. Each simulated
event is then augmented with a weight function $w_{i}(\bm{\theta})$ that is
polynomial in $\bm{\theta}$ with coefficients determined from Eq. 1, evaluated
at a sufficiently large number of basis parameter points [16]. We can choose
an overall normalization (scaling) of the weight functions to the total number
of predicted events for an arbitrarily chosen luminosity $\mathcal{L}$, i.e.,
$\sum_{i=1}^{N}w_{i}(\bm{\theta})=\mathcal{L}\sigma(\bm{\theta})$. In this
“cross section normalization”, the differential parton-level cross section,
integrated over a small phase space volume $\Delta\bm{z}$, is approximated by

$\int_{\Delta\bm{z}}\frac{\textrm{d}\sigma({\bm{\theta}})}{\textrm{d}\bm{z}}\textrm{d}\bm{z}\approx\frac{\textrm{d}\sigma({\bm{\theta}})}{\textrm{d}\bm{z}}\Delta\bm{z}\approx\frac{1}{\mathcal{L}}\sum_{z_{i}\in\Delta\bm{z}}w_{i}(\bm{\theta})$
(5)

provided that $N$ is sufficiently large and
$\textrm{d}\sigma({\bm{\theta}})/\textrm{d}\bm{z}$ does not vary strongly in
the region $\Delta\bm{z}$. The weight functions $w_{i}(\bm{\theta})$ encode
the $\bm{\theta}$-dependence of all predictions obtained from simulation.

We can also remove the $\bm{\theta}$-dependence of the sample’s normalization
by scaling $\omega_{i}=w_{i}N/(\mathcal{L}\sigma(\bm{\theta}))$, i.e., with a
factor common to all events. In this “PDF normalization”, we obtain an
approximation of the integral over $\Delta\bm{z}$ of the PDF as

$\int_{\Delta\bm{z}}p(\bm{z}|\bm{\theta})\textrm{d}\bm{z}\approx\frac{1}{N}\sum_{z_{i}\in\Delta\bm{z}}\omega_{i}(\bm{\theta}).$
(6)

Via Eq. 1–4, the weight functions are related to the likelihood ratio and the
score function evaluated at the parton-level configuration $\bm{z}_{i}$ by

$\displaystyle r(\bm{z}_{i}|\bm{\theta},\bm{\theta}_{\textrm{ref}})$
$\displaystyle=$
$\displaystyle\frac{\omega_{i}(\bm{\theta})}{\omega_{i}(\bm{\theta}_{\textrm{ref}})}=\frac{w_{i}(\bm{\theta})}{w_{i}(\bm{\theta}_{\textrm{ref}})}\frac{\sigma(\bm{\theta}_{\textrm{ref}})}{\sigma(\bm{\theta})},$
(7) $\displaystyle t(\bm{z}_{i}|\bm{\theta})$ $\displaystyle=$
$\displaystyle\bm{\nabla}_{\theta}\log\omega_{i}(\bm{\theta})=\bm{\nabla}_{\theta}\log\frac{w_{i}(\bm{\theta})}{\sigma(\bm{\theta})}.$
(8)

The gain in computational efficiency from beating the curse of dimensionality
via the simple $\bm{\theta}$-dependence in Eqs. 2–8 is enormous. For example,
the effects from single insertions of the 15 operators in top quark physics
[28] require only 136 float numbers per event to compute
$\sigma(\bm{\theta})$, $\textrm{d}\sigma({\bm{\theta}})/\textrm{d}\bm{z}$, the
likelihood, likelihood ratios, and the score vector in the full 15-dimensional
space of theory parameters. Assuming that a traditional sampling of the theory
parameter space uses a grid of $M$ points in each dimension, a factor of
$M^{15}$ is saved. The BIT algorithm will exploit this efficient structure.

Next, we turn our attention to the detector level, i.e. to observables
including the simulation of the hard interaction, the parton shower, ME
matching, the hadronization of strongly interacting partons, and the detector
response to the stable particles. The event simulation thus provides an event
with parton-level configuration $\bm{z}$ with an observable detector-level
feature vector ${\bm{x}}$. The configuration spaces of both, observed event
features $\bm{x}$ and parton-level configuration $\bm{z}$, are high-
dimensional with large numbers of continuous and discrete variables. Denoting
a detector-level selection by $j$, we obtain expected Poisson yields from the
simulated sample as $\lambda_{j}(\bm{\theta})=\sum_{\bm{x}_{i}\in
j}w_{i}(\bm{\theta})$. The augmented data set in the cross section
normalization (Eq. 5) is then simply

$\\{w_{i}(\bm{\theta}),\mathbf{x}_{i},\mathbf{z}_{i}\\}_{i=1}^{N}.$ (9)

When normalizing to the PDF (Eq. 6), we replace $w_{i}$ by $\omega_{i}$ to
obtain the equivalent data set

$\\{\omega_{i}(\bm{\theta}),\mathbf{x}_{i},\mathbf{z}_{i}\\}_{i=1}^{N}.$ (10)

A conceptual point pertains to the fact that $\bm{x}$ and $\bm{z}$ are jointly
available for each simulated event. Thus, the likelihood
$r({\bm{x}},{\bm{z}}|{\bm{\theta}_{\textrm{1}}},{\bm{\theta}}_{\textrm{2}})$
and the score vector
${\bm{t}}({\bm{x}},{\bm{z}}|{\bm{\theta}})=\nabla_{\bm{\theta}}\log
p({\bm{x}},{\bm{z}}|{\bm{\theta}})$ available from the simulation are jointly
conditional on the vectors $\bm{x}$ and $\bm{z}$ for arbitrary
${\bm{\theta}}$, ${\bm{\theta}_{\textrm{1}}}$, and
${\bm{\theta}_{\textrm{2}}}$. According to the Neyman-Pearson lemma, however,
it is the true likelihood ratio $r(\bm{x}|\bm{\theta}_{1},\bm{\theta}_{2})$
that is the optimal test statistic to discriminate between two hypotheses
$\bm{\theta}_{1}$ and $\bm{\theta}_{2}$. In the case that the parton shower,
the ME matching prescription, and the electromagnetic and strong interactions
in the detector simulation do not depend on $\bm{\theta}$, the true likelihood
can be factorized as
$p(\bm{x}|\bm{\theta})\equiv\int\,p(\bm{x},\bm{z})\textrm{d}\bm{z}=\int\,p(\bm{x}|\bm{z})p(\bm{z}|\bm{\theta})\textrm{d}\bm{z}$,
but the probability $p(\bm{x}|\bm{z})$ to observe a feature vector $\bm{x}$
given the (latent) parton-level configuration $\bm{z}$ is intractable.

So can we still exploit the simple parton-level polynomial structure given
only joint quantities while $p(\bm{x}|\bm{z})$ is intractable? The authors of
Refs. [14, 15, 16, 17, 18] show how MSE loss functions regress on the true
likelihood ratio $r({\bm{x}}|{\bm{\theta}}_{1},{\bm{\theta}}_{2})$ and the
true score ${\bm{t}}({\bm{x}}|{\bm{\theta}})=\bm{\nabla}_{\bm{\theta}}\log
p({\bm{x}}|{\bm{\theta}})$, given only the joint quantities that are also
conditional on the parton-level configuration ${\bm{z}}$. The intractable
factors in the likelihood, related to the parton shower, the ME matching, the
detector simulation, etc., are shown to cancel [16] provided that the EFT
modifications do not enter $p(\bm{x}|\bm{z})$ either directly or indirectly
via loop corrections. We will briefly return to this point in Sec. 5, where we
show that the BIT minimizes a certain MSE loss function. Meanwhile, we drop
$\bm{z}$ in the notation when it is not needed.

##  4 Decision tree fit as Fisher information optimization

Before we construct the new algorithm, it is instructive to interpret the weak
learner of a traditional BDT algorithm in the context of Poisson counting
experiments. It is defined by a maximum number $D$ of consecutive requirements
on the vector of input features $\bm{x}$, collectively denoted by
$\alpha_{j}$, that group the input data in no more than $2^{D}$ terminal nodes
$j\in\mathcal{J}$. All terminal nodes are disjoint and the union of all
$j\in\mathcal{J}$ equals the total feature space. The training minimizes a
loss function with respect to $\alpha_{j}$ which is evaluated on the
categorial input training data $\\{y_{i},\mathbf{x}_{i}\\}_{i=1}^{N}$ . For
binary classification, the labels can be chosen $y_{i}\in\\{-1,1\\}$ and the
prediction is given by the majority $F_{j}=\textrm{sgn}\left(\sum_{i\in
j}y_{i}\right)$ of the training labels in each terminal node $j$. The weak
learner’s binary prediction can then be written as

$F(\mathbf{x})=\sum_{j\in\mathcal{J}}\mathds{1}_{\alpha_{j}}(\mathbf{x})F_{j},$
(11)

where $\mathds{1}_{\alpha_{j}}(\mathbf{x})=1$ if the feature vector
$\mathbf{x}$ satisfies the requirements $\alpha_{j}$ of the terminal node $j$
and is zero otherwise. In HEP, we call training data with $y_{i}=1$
($y_{i}=-1$) signal (background) events and we denote their total number by s
(b); the purity is defined as $\rho=s/(s+b)$ and the total yield is
$\lambda=s+b$. During the training of the decision tree, a loss function

$\Delta_{H}=\lambda_{1}H(\rho_{1})+\lambda_{2}H(\rho_{2})-\lambda H(\rho)$
(12)

is minimized to recursively split nodes. Here, the primary node is
characterized by ($\rho,\lambda$) and the secondary nodes by
($\rho_{1},\lambda_{1}$) and ($\rho_{2},\lambda_{2}$). A particular choice for
$H$, the Gini impurity $H_{G}=\rho(1-\rho)$, leads to

$\Delta_{H_{G}}=-\left(\frac{1}{\lambda_{1}}+\frac{1}{\lambda_{2}}\right)^{-1}(\rho_{1}-\rho_{2})^{2}$
(13)

and implies that optimal node splits result in maximally different secondary
purities with a penalty for small yields.

Equation 13 has an interesting interpretation in HEP, first explained in Ref.
[25]. The Fisher information matrix, defined as the variance of the score
vector,

$I_{ab}(\bm{\theta})=\mathds{E}_{\bm{\theta}}\left[\frac{\partial}{\partial\theta_{a}}\log
p(\bm{x}|\bm{\theta})\frac{\partial}{\partial\theta_{b}}\log
p(\bm{x}|\bm{\theta})\right],$ (14)

has many useful applications in parameter estimation and regression problems.
The Cramér-Rao bound [29, 30], e.g., states that the variance of an unbiased
scalar estimator $\hat{\theta}$ of $\theta$ satisfies

$\textrm{var}(\hat{\theta})\geq\frac{1}{I(\theta)},$ (15)

which we can use as an approximate estimator of its expected variance.

In binned HEP data analyses, we partition the kinematic phase space at the
reconstruction level in disjoint regions and perform Poisson counting
experiments on the resulting yields. For a single Poisson yield
$\lambda(\theta)$, the Fisher information matrix becomes

$I_{ab}^{(\lambda)}(\bm{\theta})=\frac{1}{\lambda}\frac{\partial\lambda}{\partial\theta_{a}}\frac{\partial\lambda}{\partial\theta_{b}}.$
(16)

Let us specify this expression to the simple case of a single-bin measurement
of the signal cross section $\Delta\sigma_{s}$. We can write
$\lambda(\Delta\sigma_{s})=b+s(\Delta\sigma_{s})=b+\Delta\sigma_{s}\mathcal{L}=\frac{1}{\rho}\Delta\sigma_{s}\mathcal{L}$
and thus,

$I^{(\lambda)}_{\Delta\sigma_{s}}=\frac{1}{\lambda}\left(\frac{\partial\lambda}{\partial\Delta\sigma_{s}}\right)^{2}=\frac{\rho\mathcal{L}}{\Delta\sigma_{s}}=\frac{1}{(\Delta\sigma_{s})^{2}}\lambda\rho^{2}.$
(17)

We can drop the irrelevant constant factor $(\Delta\sigma_{s})^{-2}$ and
compare Eq. 17 to the $\lambda H$ terms in Eq. 12. We find a criterion
$H_{F}=-\rho^{2}$ will maximize the Fisher information $I_{\Delta\sigma_{s}}$.
Because $H_{F}$ and $H_{G}$ differ only by terms linear in $\rho$ and such
contributions drop out in Eq. 12, we conclude that the training of a decision
tree with $\Delta_{H_{G}}$ is equivalent to maximizing the Fisher information
in a cross section measurement of the signal.

##  5 Boosted Information Trees

The interpretation stated at the end of Sec. 4 suggests to define a node split
criterion for general $\lambda(\theta)$. We define a training loss in terms of
the Fisher information of a single Poisson yield by generalizing to

$H_{F}(\lambda)=-\left.\left(\frac{1}{\lambda}\frac{\partial\lambda}{\partial\theta}\right)^{2}\right|_{\theta=\theta_{0}},$
(18)

leading to a loss function

$\displaystyle\Delta_{H_{F}}$ $\displaystyle=$
$\displaystyle\lambda_{1}H_{F}(\lambda_{1})+\lambda_{2}H_{F}(\lambda_{2})-\lambda
H_{F}(\lambda)$ (19) $\displaystyle=$
$\displaystyle-I^{(\lambda_{1})}(\theta_{0})-I^{(\lambda_{2})}(\theta_{0})+I^{(\lambda)}(\theta_{0})$
$\displaystyle=$
$\displaystyle-\left.\left(\frac{1}{\lambda_{1}}+\frac{1}{\lambda_{2}}\right)^{-1}\left(\frac{\partial}{\partial\theta}\log\lambda_{1}-\frac{\partial}{\partial\theta}\log\lambda_{2}\right)^{2}\right|_{\theta=\theta_{0}}.$

This node split criterion thus maximizes the gain in the Fisher information of
a scalar parameter $\theta$ of Poisson yields or, equivalently, the difference
in the fractional rate of change
$\partial_{\theta}\log\lambda(\theta)=\lambda^{\prime}(\theta)/\lambda(\theta)$
of the secondary yields with same penalty term as in Eq. 13. Note that a
constant $\theta$-dependent factor, reflecting a $\theta$ dependence of the
total cross section common to all training events, drops out from Eq. 19.
Therefore, either version of the augmented data in Sec. 3 can be used. The
correspondence of Eq. 13 and Eq. 19 can be understood by a loose analogy
between $b\leftrightarrow\lambda(\theta)$ and
$s+b\leftrightarrow\lambda(\theta+\Delta\theta)$ such that
$\rho\leftrightarrow\lambda^{-1}\Delta\lambda/\Delta\theta$ for small
$\Delta\theta$. From Eq. 19 we can thus determine the feature space
partitioning $\mathcal{J}$ of a weak learner. It does, however, not determine
the predicted value, which will become clear when we next consider boosting.

In classification, boosting minimizes a loss function $L(y,F(\mathbf{x}))$ by
iteratively fitting weak learners to the pseudo-residuals of the preceding
iteration step. The AdaBoost.M1 algorithm [31], e.g., has been shown to
minimize the exponential loss $\mathds{E}[\exp(-yF(x))]$ where $y$ denotes the
true classification value [32]. The product $-yF(x)$ is negative for a correct
classification and corresponds to $-\frac{\partial w}{\partial\theta}F(x)$ for
an augmented training event $(w(\theta),\mathbf{x})$. We therefore make the
corresponding linear ansatz

$F(\mathbf{x})=-\sum_{j\in\mathcal{J}}\mathds{1}_{\alpha_{j}}(\mathbf{x})F_{j}\left.\frac{\partial
w}{\partial\theta}\right|_{\theta=\theta_{0}},$ (20)

where the constants $F_{j}$ are the weak learner’s prediction for a new event
that falls into the terminal node $j\in\mathcal{J}$. It is straightforward to
verify that the choice

$F_{j}=\frac{\partial}{\partial\theta}\log\sum_{i\in
j}w_{i}(\theta)\bigg{|}_{\theta=\theta_{0}}=\frac{1}{\lambda_{j}}\left.\frac{\partial\lambda_{j}}{\partial\theta}\right|_{\theta=\theta_{0}}$
(21)

results in

$\mathds{E}(F)=-\sum_{j\in\mathcal{J}}\frac{1}{\lambda_{j}}\left.\left(\frac{\partial\lambda_{j}}{\partial\theta}\right)^{2}\right|_{\theta=\theta_{0}}=-\sum_{j\in\mathcal{J}}I^{(\lambda_{j})}(\theta_{0})$
(22)

and we thus obtain a discriminator that attempts to partition the training
data to maximize the cumulative Poisson Fisher Information of the partitioned
training yields.

Because terms $\mathcal{O}\left((\theta-\theta_{0})^{2}\right)$ in the
expansion of $w(\theta)$ drop out, we simplify the input data further to the
weight and its derivative,

$\\{\mathbf{x}_{i},w_{i},w^{\prime}_{i}\\}=\left\\{\mathbf{x}_{i},w_{i}(\theta_{\textrm{0}}),\left.\frac{\partial
w_{i}}{\partial\theta}\right|_{\theta_{\textrm{0}}}\right\\},$ (23)

evaluated at the parameter point where the discriminator is trained. In
Algorithm 1, we present the pseudo-code for the fitting procedure of our
decision-tree weak learner. It is an adaption of the well-known CART algorithm
by Breiman et al. [33], where the feature space in each recursion step is
further divided by greedily selecting the dimension and cut value combination
that maximizes a local gain. The main differences are that we maximize the
Fisher information as node split criterion as described above and set the log-
derivative of the node yield as prediction corresponding to a terminal node.
Overfitting is mitigated by enforcing a maximum tree depth $D$ and a minimum
number of events $N_{\textrm{min}}$ in a terminal node. Since all events are
partitioned among the nodes of each layer in the tree, the runtime of the
algorithm scales with $\mathcal{O}(dDN\log N)$, where $d$ is the dimension of
the features and $N$ the number of events. It can be efficiently implemented
by sorting the events’ features for each dimension separately and using
cumulative sums over the weights and weight derivatives to find the locally
optimal cut. The terminal nodes $j\in\mathcal{J}$ are encoded each as tuple
consisting of the requirements $\alpha_{j}$ (selecting a subset of the data
$\mathcal{D}_{\alpha}\subset\mathcal{D}$) and prediction $F_{j}$. They
determine the final weak learner
$F(\mathbf{x})=\sum_{j\in\mathcal{J}}\mathds{1}_{\alpha_{j}}(\mathbf{x})\cdot
F_{j}$.

Data: Data set
$\mathcal{D}=\\{\mathbf{x}_{i},w_{i},w^{\prime}_{i}\\}_{i=1}^{N}$.

Input: Tree depth $D$, minimum terminal node size $N_{\mathrm{min}}$

Output: Weak learner $F$

$\alpha\leftarrow(),Q\leftarrow(\alpha),\mathcal{J}\leftarrow\\{\\},P\leftarrow\\{1,\dots,d\\}$

while  _$Q\neq\emptyset$_ do

$\alpha\leftarrow\mathsf{pop}(Q)$ $\bm{\pi}\leftarrow\mathrm{arg\,sort}\
\mathbf{x}\in\mathcal{D}_{\alpha}$ {sorted indices for each dimension $p\in
P$} if  _$|\alpha|\leq D\land|\mathcal{D}_{\alpha}| >2N_{\textrm{min}}$_ then

$K\leftarrow\\{N_{\textrm{min}},\dots,|\mathcal{D}_{\alpha}|-N_{\textrm{min}}\\}$
{allowed cut points} $p^{*},k^{*}\leftarrow\operatorname*{arg\,max}_{p\in
P,k\in
K}\left[\frac{\left(\sum_{i=1}^{k}w^{\prime}_{\bm{\pi}_{p,k}}\right)^{2}}{\sum_{i=1}^{k}w_{\bm{\pi}_{p,k}}}+\frac{\left(\sum_{i=k+1}^{|\mathcal{D}_{\alpha}|}w_{\bm{\pi}_{p,k}}\right)^{2}}{\sum_{i=k+1}^{|\mathcal{D}_{\alpha}|}w_{\bm{\pi}_{p,k}}}\right]$
$c\leftarrow x_{\bm{\pi}_{p^{*},k^{*}},p^{*}}$ from $\mathcal{D}_{\alpha}$ if
_$c$ is a valid cut_ then

$\alpha_{\mathrm{L}}\leftarrow\alpha\cup(p^{*},\leq,c)$
$\alpha_{\mathrm{R}}\leftarrow\alpha\cup(p^{*},>,c)$ $Q\leftarrow
Q\cup(\alpha_{\mathrm{L}},\alpha_{\mathrm{R}})$

else

$\mathcal{J}\leftarrow\mathcal{J}\cup\left(\alpha,\frac{\sum_{w^{\prime}\in\mathcal{D}_{\alpha}}w^{\prime}}{\sum_{w\in\mathcal{D}_{\alpha}}w}\right)$

end if

else

$\mathcal{J}\leftarrow\mathcal{J}\cup\left(\alpha,\frac{\sum_{w^{\prime}\in\mathcal{D}_{\alpha}}w^{\prime}}{\sum_{w\in\mathcal{D}_{\alpha}}w}\right)$

end if

end while

return
$F(\mathbf{x})=\sum_{j\in\mathcal{J}}\mathds{1}_{\alpha_{j}}(\mathbf{x})\cdot
F_{j}$

Algorithm 1 The $\mathsf{fit}$ procedure for the weak learner

Choosing a learning rate $\eta$ and a number of boosting iterations $B$, the
final algorithm is described in Algorithm 2. In each iteration $b$, only a
fraction $\eta f^{(b)}$ is retained from the weak learner for the cumulant
prediction $F^{(b)}$. The replacement $w_{i}^{\prime}\rightarrow
w^{\prime}_{i}-\eta w_{i}F^{(b-1)}$ implements the corresponding reduction in
the local score estimates by reducing the weight derivatives for the next
iteration.

Data: Data set
$\mathcal{D}=\\{\mathbf{x}_{i},w_{i},w^{\prime}_{i}\\}_{i=1}^{N}$.

Input: Number of boosting iterations $B$, learning rate $\eta$

Output: Boosted learner $F^{(B)}$

$F^{(0)}\leftarrow 0$ {initialize boosted learner}

for  _$b\leftarrow 1,\ldots,B$_ do

$f^{(b)}\leftarrow\mathsf{fit}\left(\left\\{\mathbf{x}_{i},w_{i},w^{\prime}_{i}-\eta
w_{i}F^{(b-1)}(\mathbf{x}_{i})\right\\}_{i=1}^{{N}}\right)$ $F^{(b)}\leftarrow
F^{(b-1)}+\eta f^{(b)}$

end for

Algorithm 2 Boosted Information Tree

Equation 21 can indeed be interpreted as an approximation of the score
function in the current node of the decision tree, iteratively refined by the
boosting algorithm. For any region $j$, we can compute the local mean of the
score function of the PDF as

$\frac{\mathds{E}_{\theta_{0}}\left(\mathds{1}_{j}(\bm{x})\cdot\partial_{\theta}\log
p(\bm{x}|\theta)\Big{|}_{\theta=\theta_{\textrm{0}}}\right)}{\mathds{E}_{\theta_{0}}\left(\mathds{1}_{j}(\bm{x})\right)}=\partial_{\theta}\log\lambda\Big{|}_{\theta=\theta_{\textrm{0}}}-\partial_{\theta}\log\sigma\Big{|}_{\theta=\theta_{\textrm{0}}},$
(24)

which allows us to understand the discriminator output as

$F^{(b)}(\bm{x})\simeq
t(\bm{x}|\theta_{0})+\partial_{\theta}\log\sigma(\theta)\Big{|}_{\theta=\theta_{\textrm{0}}}.$
(25)

The last term is a calculable bias term accounting for the $\theta$ dependence
of the inclusive cross section. It does not change the performance of the
discriminator in the analysis setting; in fact, the output can be normalized
to the unit interval and this option is available in the implementation. If
the augmented data set is at the PDF level (i.e. $\omega_{i}$ instead of
$w_{i}$ without the factor $\sigma(\theta)$ as discussed in Sec. 3), the last
term in Eq. 25 is absent.

Training simulation is needed only for a single reference parameter point
which can, but need not be at $\theta_{0}$. The discriminator learns only from
the leading term in the expansion around the reference point in parameter
space; contributions $\mathcal{O}\left((\theta-\theta_{0})^{2}\right)$ do not
enter Algorithm 1. In this linear approximation, the likelihood can be written
as

$p(\bm{x}|\theta)=\frac{1}{Z(\theta)}p(t(\bm{x}|\theta_{0})|\theta_{0})\exp(t(\bm{x}|\theta_{0})\cdot(\theta-\theta_{0})),$
(26)

which is an exponential family. Therefore, the score $t(\bm{x}|\theta)$ is a
sufficient statistic for $p(x|\theta)$, learned by the algorithm [16].

It is interesting to derive the node split criterion from a different starting
point. Let us disregard the Fisher information and instead regress in the
score. For a given node split at a value $\alpha$, we denote the regressed
values in the nodes, defined by $x<\alpha$ and $x\geq\alpha$, as $F_{L}$ and
$F_{R}$, respectively. With an ansatz for the weak regressor
$\hat{t}(x|\alpha,F_{L},F_{R})=F_{L}\Theta(x<\alpha)+F_{R}\Theta(x\geq\alpha)$,
we can minimize the MSE loss

$\textrm{MSE}(\alpha,F_{L},F_{R})=\mathds{E}_{\theta_{0}}\left[\left(\frac{1}{\omega(\theta)}\frac{\partial\omega}{\partial\theta}\Big{|}_{\theta=\theta_{\textrm{0}}}-\hat{t}(x|\alpha,F_{L},F_{R})\right)^{2}\right].$
(27)

The equations $\frac{\partial\textrm{MSE}}{\partial
F_{L}}=\frac{\partial\textrm{MSE}}{\partial F_{R}}=0$ can be solved for
$F_{L}$ and $F_{R}$ and, dropping terms independent of the threshold $\alpha$,
we arrive at

$\displaystyle\textrm{MSE}(\alpha)$ $\displaystyle=$
$\displaystyle-\frac{1}{\lambda_{L}(\alpha)}\left(\frac{\partial\lambda_{L}(\alpha)}{\partial\theta}\right)^{2}\Bigg{|}_{\theta=\theta_{\textrm{0}}}-\frac{1}{\lambda_{R}(\alpha)}\left(\frac{\partial\lambda_{R}(\alpha)}{\partial\theta}\right)^{2}\Bigg{|}_{\theta=\theta_{\textrm{0}}}$
(28) $\displaystyle=$
$\displaystyle-I^{(\lambda_{L}(\alpha))}(\theta_{0})-I^{(\lambda_{R}(\alpha))}(\theta_{0}),$

where $\lambda_{L(,R)}=\sum_{i\in L(,R)}\omega_{i}(\theta)$. The loss in Eq.
27 is of the general form of Eq. 27 in Ref. [16] which implies we regress on
the true score $t(\bm{x}|\theta)$ despite the augmented data being conditional
on $\bm{z}$. Equation 28 shows that regressing in the MSE loss is equivalent
to maximizing the cumulative Fisher information as a function of the threshold
value $\alpha$, as is done in Eq. 22 for the BIT. While the algorithm does
approximate the score, there is no numerical minimization needed for the
regressed value. Computationally, the algorithm behaves as a BDT.

We note that Eq. 27 suggests to extend the algorithm to estimate higher
derivative terms of the likelihood, starting with the quadratic term,
neglected here. Finally, in the case of more than one parameter, we can also
expand the true likelihood ratio, i.e., the optimal test-statistic according
to the Neyman-Pearson lemma, as

$r(\bm{x}|\bm{\theta},\bm{\theta}_{\textrm{ref}})=1+(\bm{\theta}-\bm{\theta}_{\textrm{ref}})\cdot\bm{t}(\bm{x}|\bm{\theta}_{\textrm{ref}})+\mathcal{O}(\bm{\theta}-\bm{\theta}_{\textrm{ref}})^{2}.$
(29)

It implies that the combination
$(\theta_{a}-\theta_{a,\textrm{ref}})t_{a}(\bm{x}|\bm{\theta}_{\textrm{ref}})$
is optimal in the linear approximation for a simple hypothesis test of
$\bm{\theta}$ against $\bm{\theta}_{\textrm{ref}}$. The components
$t_{a}(\bm{x}|\bm{\theta})$ of the score vector can be separately estimated by
BITs, trained with the corresponding partial derivatives of the weight
functions, i.e.,
$\partial_{a}w_{i}(\bm{\theta})|_{\bm{\theta}_{\textrm{ref}}}$ instead of
$w^{\prime}$ in Eq. 23. We leave these developments to future work.

##  6 Fitting analytic toy models

We now test the BIT on a variety of analytic toy models listed in Table 1. We
choose an exponentially falling PDF, a power-law model, a Gaussian model where
the parameter is the mean, a Gaussian model where the parameter is the width,
and a mixture model. In Fig. 1 (left), we show the PDFs that we use to
generate a data set consisting of $10^{5}$ events in each case. In Fig. 1
(center), we compare the theoretical calculation with the predicted score at
different numbers of boosting iterations $b$. The prediction nicely
approximates the theoretically calculated score functions $t(\bm{x}|\theta)$
for increasing $b$. The distributions of the predicted score for the training
data sets and statistically independent test data sets of the same size,
evaluated for $\theta_{0}$ and a nearby parameter point
$\theta_{0}+\Delta\theta$, are shown in Fig. 1 (right). No significant
overtraining is observed in the predicted distributions.

In Fig. 2, we show the same distributions for a mixture of two exponential
PDFs with different parameters $\alpha_{1,2}$,

$p(x|\theta)\propto\left(\exp(-\alpha_{1}(x-x_{0}))+\theta\exp(-\alpha_{2}(x-x_{0}))\right)^{2}.$
(30)

The dependence on the mixing parameter $\theta$ is broadly reminiscent of the
structure in Eq. 1 while the choice of the exponential is arbitrary. The same
conclusions as for the other toy models hold, raising the prospect of good
performance in realistic EFT scenarios.

Table 1: PDF $p(x|\theta)$, score function $t(x|\theta)$, support of the PDF,
constant parameters, and $\theta_{0}$ for the one-dimensional toy models.

Exponential
---
$\displaystyle p(x|\theta)$ $\displaystyle=$ $\displaystyle\theta
e^{-\theta(x-x_{0})},\;\;x>x_{0},\;\;x_{0}=25,\;\;\theta_{0}=0.01$
$\displaystyle t(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\theta}-(x-x_{0})$
Power law
$\displaystyle p(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{\theta-1}{x_{0}}\left(\frac{x}{x_{0}}\right)^{-\theta},\;\;x>x_{0},\;\;x_{0}=100,\;\;\theta_{0}=3$
$\displaystyle t(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\theta-1}-\log\frac{x}{x_{0}}$
Gaussian mean
$\displaystyle p(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2}\left(\frac{x-\theta}{\sigma}\right)^{2}\right),\;\;\sigma=1,\;\;\theta_{0}=0$
$\displaystyle t(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{x-\theta}{\sigma^{2}}$
Gaussian width
$\displaystyle p(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}\theta}\exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\theta}\right)^{2}\right),\;\;\mu=0,\;\;\theta_{0}=1$
$\displaystyle t(x|\theta)$ $\displaystyle=$
$\displaystyle\frac{1}{\theta}\left(\left(\frac{x-\mu}{\theta}\right)^{2}-1\right)$
Mixture
$\displaystyle p(x|\theta)$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2\alpha_{1}}+\frac{2\theta}{\alpha_{1}+\alpha_{2}}+\frac{\theta^{2}}{2\alpha_{2}}\right)^{-1}\left(e^{-\alpha_{1}(x-x_{0})}+\theta
e^{-\alpha_{2}(x-x_{0})}\right)^{2}$ $\displaystyle t(x|\theta)$
$\displaystyle=$
$\displaystyle\frac{2}{\theta+e^{-(x-x_{0})(\alpha_{1}-\alpha_{2})}}-2\frac{\theta\alpha_{1}(\alpha_{1}+2\alpha_{2})+2\alpha_{1}\alpha_{2}}{\alpha_{2}(\alpha_{1}+\alpha_{2})+4\alpha_{1}\alpha_{2}\theta+\alpha_{1}(\alpha_{1}+\alpha_{2})\theta^{2}}$
$\displaystyle\alpha_{1}$ $\displaystyle=$ $\displaystyle
2\alpha_{2}=0.02,\;\;x_{0}=20,\;\;\theta_{0}=0$

Figure 1: The PDF for toy models in Table 1 for $\theta_{0}$ and, for
comparison, for an alternative parameter $\theta_{0}+\Delta\theta$ (left).
Analytically computed score function $t(x,\theta_{0})$ and the output of the
BIT algorithm for different numbers of boosting iterations $b$ (center).
Distributions of the predicted score function for the training data set and a
statistically independent test data set for $\theta_{0}$ and
$\theta_{0}+\Delta\theta$ (right). From top to bottom, the exponential, the
power-law, the Gaussian-mean, and the Gaussian-width model are shown.

Figure 2: The PDF for the mixture model in Table 1 for $\theta_{0}$ and, for
comparison, for an alternative parameter $\theta_{0}+\Delta\theta$ (left).
Analytically computed score function $t(x,\theta_{0})$ and the output of the
BIT algorithm for different numbers of boosting iterations $b$ (center).
Distributions of the predicted score function for the training data set and a
statistically independent test data set for $\theta_{0}$ and
$\theta_{0}+\Delta\theta$ (right).

Figure 3: Overtraining measure $L(\mathcal{D},b)$ as a function of the
boosting iteration (left). There is no significant difference between the
training and the statistically independent test data set. Contour lines of the
analytically computed score of the two-dimensional exponential model (dashed
lines) and the colored contour lines of the predicted score (right).

Overtraining can be monitored by evaluating the loss defined as

$L(\mathcal{D},b)=-\sum_{i\in\mathcal{D}}\left(\frac{\partial\omega_{i}(\theta)}{\partial\theta}F^{(b)}(\mathbf{x}_{i})\right)$
(31)

for the training and a statistically independent test data set. In Fig. 3
(left) we show its evolution as a function of the boosting iteration. After
about 30 iterations, the loss no longer reduces. A statistically independent
data set shows the exact same behavior. The reason for the negligible
overtraining lies in the homogeneity of the training data set in the toy case
which is sampled from a simple analytic PDF. The algorithm cannot place
thresholds based on statistically insignificant fluctuations in the training
data set as in BDT classification: There are no signal and background events
with independent stochastics. We expect the overtraining situation to
deteriorate slightly when the variance of the $w_{i}$ is large in the training
data set. In an extreme situation where $w_{i}^{\prime}/w_{i}$ takes only two
values, we expect the overtraining situation to resemble the classification
case. In all studies conducted so far, a simple requirement $N\geq
N_{\textrm{min}}$ on the events in each node of any of the weak learners is
found to be sufficient to suppress significant overtraining. Further measures
to mitigate overtraining should be informed by studies with realistic SM-EFT
models that we leave to future work.

Because of the similarity of the BIT algorithm to BDT classification, we
expect that other well-known differences of tree-boosting and neural-network
based strategies prevail for the BIT. Besides the reduced risk of overtraining
already discussed, these include, for example, a generally faster execution
time of the trained estimator, no intrinsic limitation in the approximation of
large gradients in feature space, and a protection against extrapolation to
regions in feature space without training data.

Finally, we briefly test the performance in models with higher feature
dimensions. Adding up to 25 random variables per feature does not affect the
output of the one-dimensional analytic models in any way. A more challenging
test of the rectangularly selecting boosted learners is to disperse the
discriminating information in a non-linear way. We consider a two-dimensional
exponential model and train a BIT on $10^{6}$ events with $B=100$. A
comparison of the theoretically calculated contour lines of the score function
and the contour lines of the regressed score is shown in Fig. 3 (right). The
main features are very well described by the discriminator. The test has been
repeated with up to 5-dimensional variants of exponential and power-law
models. After moderate adjustments of training statistic and $B$, the fit
quality to the score function is similar as in the one-dimensional cases.

We leave possible refinements of the algorithm to further mitigate the
overtraining, a possible implementation in a versatile generic tree-boosting
framework such as XGBoost [34], as well as a potential extension to quadratic
terms in the likelihood expansion to future work.

##  7 Conclusions and outlook

We presented the “Boosted Information Tree” algorithm, a boosted tree-based
discriminator that optimizes the Fisher information in the node split of the
weak learner. It regresses on the true score of the model at the computational
cost of a decision tree. It provides a sufficient statistic in the linear
approximation around the reference point in parameter space where the
discriminator is trained and can be used to improve analysis strategies for
measuring EFT parameters. Overtraining is, in general, less critical than in
the case of Boosted Decision Trees. The algorithm is validated on one- and
two-dimensional analytically tractable toy models. An implementation of the
algorithm is available on GitHub [23], and an introduction to the interface is
provided in A.

## Acknowledgements

R. S. and L. L. want to thank Andrea Valassi, author of Refs. [25, 26], for
useful discussions in the early stages of this project.

#### Funding:

The work of D. S. and L.L. was supported by the Austrian Science Fund (FWF)
projects P33771 and P31578, respectively.

##  Appendix A Python implementation and interface

An implementation of the BIT algorithm in both Python 2 and Python 3 is
available on GitHub [23], including a demonstration using Jupyter Notebook.
The python package NumPy [35] is the only external dependency.

#### Preparation

Import the Python modules, for Python 2 import numpy as np from
BoostedInformationTree import BoostedInformationTree and for Python 3: import
numpy as np from BoostedInformationTreeP3 import BoostedInformationTree

Three two-dimensional NumPy arrays are needed as input: the features, the
event weights, and the corresponding weight derivatives as defined in Eq. 23.
Training and test data sets are available for the toy models. The power-law
model is loaded with: data_dir = ’data’ def load_data(name): return
np.loadtxt(’

features = load_data(’training_features_power_law_model’) features =
features.reshape(features.shape[0], -1) weights =
load_data(’training_weights_power_law_model’) diffs =
load_data(’training_diff_weights_power_law_model’)

#### Hyperparameters

The training configuration (hyperparameters) comprises the weak learner’s
maximum depth $D$ (max_depth) and its minimum event count in the terminal
nodes (min_size). The boosting constructs a number $B$ (n_tree) of weak
learners to approximate the score. The learning rate $\eta$ is specified by
the parameter learning_rate. Setting the boolean flag calibrated to True
calibrates the output to the unit interval. If it is set to False, the
approximation of the score function in Eq. 25 is returned. For the toy models,
a reasonable starting point is: learning_rate = 0.02 n_trees = 100
learning_rate = 0.2 max_depth = 2 min_size = 50 calibrated = False

#### Training & Prediction

To start the training, the boost method of a suitably constructed
BoostedInformationTree object is called: bit = BoostedInformationTree(
training_features = training_features, training_weights = training_weights,
training_diff_weights = training_diff_weights, learning_rate = learning_rate,
n_trees = n_trees, max_depth = max_depth, min_size = min_size, calibrated =
calibrated)

bit.boost()

The trained discriminator can be used to predict the score for a new feature
vector or an array of feature vectors:

test_feat = load_data(’test_features_power_law_model’) test_feat =
test_feat.reshape(test_feat.shape[0], -1)

bit.predict(test_feat[0]) # returns one prediction
bit.vectorized_predict(test_feat) # returns all predictions

## References

## References

  * [1] C. P. Burgess, Introduction to effective field theory, Ann. Rev. Nucl. Part. Sci. 57 (2007) 329.  arXiv:hep-th/0701053, doi:10.1146/annurev.nucl.56.080805.140508.
  * [2] G. F. Giudice, C. Grojean, A. Pomarol, R. Rattazzi, The strongly-interacting light Higgs, JHEP 06 (2007) 045.  arXiv:hep-ph/0703164, doi:10.1088/1126-6708/2007/06/045.
  * [3] B. Grinstein, M. B. Wise, Operator analysis for precision electroweak physics, Phys. Lett. B 265 (1991) 326.  doi:10.1016/0370-2693(91)90061-T.
  * [4] I. Brivio, M. Trott, The standard model as an effective field theory, Phys. Rept. 793 (2019) 1.  arXiv:1706.08945, doi:10.1016/j.physrep.2018.11.002.
  * [5] D. de Florian, et al., Handbook of LHC Higgs cross sections: 4. Deciphering the nature of the Higgs sector, CERN Yellow Rep. Monogr. 2 (2016) 1.  arXiv:1610.07922, doi:10.23731/CYRM-2017-002.
  * [6] B. Grzadkowski, M. Iskrzynski, M. Misiak, J. Rosiek, Dimension-six terms in the standard model Lagrangian, JHEP 10 (2010) 085.  arXiv:1008.4884, doi:10.1007/JHEP10(2010)085.
  * [7] A. Falkowski, D. Straub, Flavourful SMEFT likelihood for Higgs and electroweak data, JHEP 04 (2020) 066.  arXiv:1911.07866, doi:10.1007/JHEP04(2020)066.
  * [8] J. De Blas, et al., HEPfit: a code for the combination of indirect and direct constraints on high energy physics models, Eur. Phys. J. C 80 (5) (2020) 456.  arXiv:1910.14012, doi:10.1140/epjc/s10052-020-7904-z.
  * [9] J. Ellis, M. Madigan, K. Mimasu, V. Sanz, T. You, Top, Higgs, diboson and electroweak fit to the standard model effective field theory, JHEP 04 (2021) 279\.  arXiv:2012.02779, doi:10.1007/JHEP04(2021)279.
  * [10] S. Dawson, S. Homiller, S. D. Lane, Putting standard model EFT fits to work, Phys. Rev. D 102 (5) (2020) 055012.  arXiv:2007.01296, doi:10.1103/PhysRevD.102.055012.
  * [11] J. J. Ethier, G. Magni, F. Maltoni, L. Mantani, E. R. Nocera, J. Rojo, E. Slade, E. Vryonidou, C. Zhang, Combined SMEFT interpretation of Higgs, diboson, and top quark data from the LHC, JHEP 11 (2021) 089.  arXiv:2105.00006, doi:10.1007/JHEP11(2021)089.
  * [12] J. J. Ethier, R. Gomez-Ambrosio, G. Magni, J. Rojo, SMEFT analysis of vector boson scattering and diboson data from the LHC Run II, Eur. Phys. J. C 81 (6) (2021) 560.  arXiv:2101.03180, doi:10.1140/epjc/s10052-021-09347-7.
  * [13] W. Buchmuller, D. Wyler, Effective Lagrangian analysis of new interactions and flavor conservation, Nucl. Phys. B 268 (1986) 621.  doi:10.1016/0550-3213(86)90262-2.
  * [14] K. Cranmer, J. Pavez, G. Louppe, Approximating likelihood ratios with calibrated discriminative classifiers, arXiv:1506.02169.
  * [15] J. Brehmer, K. Cranmer, G. Louppe, J. Pavez, Constraining effective field theories with machine learning, Phys. Rev. Lett. 121 (2018) 111801.  arXiv:1805.00013, doi:10.1103/PhysRevLett.121.111801.
  * [16] J. Brehmer, K. Cranmer, G. Louppe, J. Pavez, A guide to constraining effective field theories with machine learning, Phys. Rev. D 98 (2018) 052004.  arXiv:1805.00020, doi:10.1103/PhysRevD.98.052004.
  * [17] J. Brehmer, G. Louppe, J. Pavez, K. Cranmer, Mining gold from implicit models to improve likelihood-free inference, Proc. Nat. Acad. Sci. 117 (2020) 5242.  arXiv:1805.12244, doi:10.1073/pnas.1915980117.
  * [18] J. Brehmer, F. Kling, I. Espejo, K. Cranmer, MadMiner: Machine learning-based inference for particle physics, Comput. Softw. Big Sci. 4 (2020) 3.  arXiv:1907.10621, doi:10.1007/s41781-020-0035-2.
  * [19] J. H. Friedman, Greedy function approximation: A gradient boosting machine., Annals Statist. 29 (2001) 1189.  doi:10.1214/aos/1013203451.
  * [20] J. H. Friedman, Stochastic gradient boosting, Comput. Stat. Data Anal. 38 (2002) 367.  doi:10.1016/S0167-9473(01)00065-2.
  * [21] L. Mason, J. Baxter, P. Bartlett, M. Frean, Boosting algorithms as gradient descent, NIPS’99, MIT Press, 1999.
  * [22] L. Mason, J. Baxter, P. Bartlett, M. Frean, Boosting algorithms as gradient descent in function space (1999).
  * [23] Boosted Information Trees: Implementation in Python2 and Python3, https://github.com/BIT4EFT/BIT (2021).
  * [24] P. De Castro, T. Dorigo, INFERNO: Inference-Aware Neural Optimisation, Comput. Phys. Commun. 244 (2019) 170.  arXiv:1806.04743, doi:10.1016/j.cpc.2019.06.007.
  * [25] A. Valassi, Binary classifier metrics for optimizing HEP event selection, EPJ Web Conf. 214 (2019) 06004.  doi:10.1051/epjconf/201921406004.
  * [26] A. Valassi, Optimising HEP parameter fits via Monte Carlo weight derivative regression, EPJ Web Conf. 245 (2020) 06038.  arXiv:2003.12853, doi:10.1051/epjconf/202024506038.
  * [27] I. Brivio, SMEFTsim 3.0 — A practical guide, JHEP 04 (2021) 073.  arXiv:2012.11343, doi:10.1007/JHEP04(2021)073.
  * [28] D. Barducci, C. Degrande, G. Durieux, F. Maltoni, E. Vryonidou, Interpreting top-quark LHC measurements in the standard-model effective field theory, arXiv:1802.07237.
  * [29] C. R. Rao, Information and the accuracy attainable in the estimation of statistical parameters, Bull. Calcutta Math. Soc. 37 (1945) 81.
  * [30] H. Cramér, Mathematical methods of statistics, Princeton University Press, 1946.
  * [31] Y. Freund, R. E. Schapire, A decision-theoretic generalization of on-line learning and an application to boosting, Journal of Computer and System Sciences 55 (1997) 119.  doi:10.1006/jcss.1997.1504.
  * [32] J. Friedman, T. Hastie, R. Tibshirani, Additive logistic regression: A statistical view of boosting (with discussion and a rejoinder by the authors), Annals Statist. 28 (2000) 337.  doi:10.1214/aos/1016218223.
  * [33] L. Breiman, J. H. Friedman, R. A. Olshen, C. J. Stone, Classification and regression trees, Wadsworth Publishing Company, 1984.
  * [34] T. Chen, C. Guestrin, XGBoost: A scalable tree boosting system, KDD ’16, Association for Computing Machinery, 2016.  arXiv:1603.02754, doi:10.1145/2939672.2939785.
  * [35] C. R. Harris, et al., Array programming with NumPy, Nature 585 (2020) 357.  doi:10.1038/s41586-020-2649-2.