Optimal statistical inference in the presence of systematic uncertainties using neural network optimization based on binned Poisson likelihoods with nuisance parameters

Data analysis in science, e.g., high-energy particle physics, is often subject to an intractable likelihood if the observables and observations span a high-dimensional input space. Typically the problem is solved by reducing the dimensionality using feature engineering and histograms, whereby the latter technique allows to build the likelihood using Poisson statistics. However, in the presence of systematic uncertainties represented by nuisance parameters in the likelihood, the optimal dimensionality reduction with a minimal loss of information about the parameters of interest is not known. This work presents a novel strategy to construct the dimensionality reduction with neural networks for feature engineering and a differential formulation of histograms so that the full workflow can be optimized with the result of the statistical inference, e.g., the variance of a parameter of interest, as objective. We discuss how this approach results in an estimate of the parameters of interest that is close to optimal and the applicability of the technique is demonstrated with a simple example based on pseudo-experiments and a more complex example from high-energy particle physics.


Introduction
Measurements in many areas of research like, e.g., highenergy particle physics, are typically based on the statistical inference of one or more parameters of interest defined by the likelihood L(D, θ) with the observables x ∈ X ⊆ R d building the dataset D = {x 0 , ..., x n } ⊆ R n×d and the parameters θ of the statistical model. The likelihood would have to be evaluated for the dataset D spanning a high-dimensional input space, which is computational expensive and typically unfeasible. The dimension of D can be reduced by the engineering of high-level observables and the usage of summary statistics. Analysts create high-level observables to reduce the dimension d of a single observation to k, ideally without losing information about the parameters θ. An example from high-energy particle physics is the usage of an invariant mass of a decay system instead of the kinematic properties of all its constituents. The dimension n of D can be reduced with the application of a summary statistic, for which histograms are frequently used so that the statistical model can be expressed in form of a likelihood, based on Poisson statistics. The dimension is thus reduced from the number of observations n to the number of bins h in the histogram, whereby the analyst tries to lose as little information as possible by optimizing the location and number of the bins. Applying both methods, the initial dimension of D ⊆ R n×d is reduced to R h×k . This paper discusses an analysis strategy using machine learning techniques, by which the suboptimal performance introduced by the reduction of dimensionality can be avoided resulting in estimates of the parameters of interest µ ∈ θ close to optimal. We put emphasis on the applicability of this approach to analyses commonly performed in high-energy particle physics at the Large Hadron Collider (LHC) [1,2] like the discovery of the Higgs boson in 2012 [3,4].
Section 2 presents the method in detail and section 3 puts the proposed technique in context of related work. Section 4 shows the performance of the method with a simple example using pseudo-experiments of a twocomponent mixture model with signal and background and section 5 applies the same approach to a more complex example from high-energy particle physics.

Method
The method is built on top of an initial dataset D ⊆ R n×d used for the statistical inference of the parameters of interest with n being the number of observations and d the number of observables. To simplify the statistical evaluation, we want to reduce the number of observables by the engineering of high-level observables. Besides manual crafting of such features, a suited approach taken from machine learning is using a neu- with ω being the free parameters. After application of the NN, we get a transformed dataset D NN ⊆ R n×k with k the number of output nodes of the NN architecture.
To reduce the dataset D NN further, the number of observations n is compressed using a histogram. Histograms are widely used as a summary statistic since counts are well described by the Poisson statistic and therefore well suited to build the statistical model of the analysis. For example in high-energy particle physics, many statistical models and well established methods for describing systematic uncertainties are based on binned Poisson likelihoods and could only be replaced with an enormous effort. The resulting dataset is D H ⊆ R h×k using h number of bins for the k-dimensional histogram. The count operation for a single bin in the histogram can be written as C = n i=0 S(x i ) with So that we can propagate the gradient from the result of the statistical inference to the free parameters ω of the NN, the histogram has to be differentiable. Since the derivative of S is ill-defined on the edges of the bin and otherwise zero, the gradient is not suitable for optimization. Therefore, we use a smoothed approximation of the gradient [5] shown in figure 1 for a one-dimensional bin. The approximation uses the similarity of S to a Gaussian function G normalized to max(G) = 1 with the standard deviation being the half-width of the bin. We replace only the gradient of the operation S and not the calculation of the count itself. On top of the reduced dataset D H , we build the statistical model using a binned likelihood L(D H , θ) with θ being the parameters of the statistical model. For a mixture model with the two processes signal s and background b, the binned likelihood describing the statistical component is given by with P being the Poisson distribution, d the observation and µ ∈ θ the parameter of interest scaling the expectation of the signal process s. Moreover, the formulation of the statistical model allows to implement systematic uncertainties by adding nuisance parameters to the set of parameters θ. For the model in equation 2, a single nuisance parameter η controlling a systematic variation ∆ of the expected bin contents results in with N being a normalized Gaussian constraining the nuisance η. If the systematic variation is asymmetric, the additional nuisance term can be written as max (η, 0)∆ up + min (η, 0)∆ down (4) or with any other differential formulation [6]. The performance of an analysis is measured in terms of the variance of the estimate for the parameters of interest, for example in our case the variance of the estimated signal strength µ. We built a differential estimate of the variance using the Fisher information [7] of the likelihood in equation 3 given by Because the maximum likelihood estimator is asymptotically efficient [8,9], the variance of the estimates for θ is asymptotically close to Assuming the first diagonal element to correspond to the parameter of interest µ, without loss of generality, the loss function to optimize the variance of the estimate for µ with respect to the free parameters ω of the NN function f is V 00 .
To be independent of the statistical fluctuations of the observation, the optimization is performed on an Asimov dataset [1]. This artificial dataset replaces the observation d with the nominal expectation for s and b serving as representative for the median expected outcome of the analysis in presence of the signal plus background hypothesis.
Given the assumption that the dimensionality reduction performed by the NN together with the histogram is a sufficient statistic, the optimization can find a function for f that gives the best estimate for the parameter of interest µ, similar to a statistical inference performed on the initial high-dimensional dataset D with an unbinned likelihood.
A graphical overview of the proposed method is given in figure 2.

Related work
The approach in [10] identifies the problem that a histogram has no suitable derivative and therefore replaces the summary statistic by means of a softmax function, which shares similarities with a count. This does not allow to construct trivially a likelihood on top of the summary and is therefore only suited for likelihood-free inference problems. Further, the approach parametrizes the systematic variations in the input space D and not on the reduced space D H after the application of the summary statistic, which omits the possibility to describe systematic variations with reweighting techniques in form of statistical weights. These points constrain the applicability of the approach. We circumvent this problem by leaving the computation of the histogram unchanged but approximating the gradient.
The strategy to allow a NN to find the best compression of the data has been also discussed in [11]. This approach shows that the NN is able to learn a summary statistic that is a close approximation of a sufficient statistic, yielding a powerful statistical inference. Similar to [10], the method is based on likelihoodfree inference restricting the applicability for our target use-cases in high-energy particle physics at the LHC.
A related approach to include systematic uncertainties in the training of the NN is the explicit decorrelation against the systematic variation. For example, the idea has been discussed on the basis of an adversarial architecture [12] and an approach penalizing the variation using approximated bin counts [5]. These strategies are not aware of the analysis objective such as the variance of a parameter of interest and therefore the decorrelation is subject to manual optimization. For a large number of nuisances, this optimization procedure is computational expensive and typically unfeasible.

Application to a simple example based on pseudo-experiments
A simple example based on pseudo-experiments and a known likelihood in the input space R n×d is used to illustrate our approach. The distributions of the signal and background components in the input space are shown in figure 3. We assume a systematic uncertainty on the mean of the background process modelled by the shifts x 2 ± 1, representing the systematic variations in equation 4.
The NN architecture is a fully-connected feed-forward network with 100 nodes in one hidden layer. The initialization follows the Glorot algorithm [13] and the activation function is a rectified linear unit [14]. The output layer has a single node with a sigmoid activation function.
We use eight bins for the histogram of the NN output and compute the variance of the estimate for the parameter of interest µ denoted by V 00 . The operations are implemented using TensorFlow as computational graph library [15,16] and we use the provided automatic differentiation and the Adam algorithm [17] to optimize the free parameters ω with the objective to minimize V 00 . The systematic variations ∆ can be implemented with reweighting techniques using statistical weights or duplicates of the nominal dataset with the simulated variations, whereas we chose the latter solution. Each gradient step is performed on the full dataset with 10 5 simulated events for each process. The training is stopped if the loss has not improved for 100 gradient steps eventually using the model with the smallest loss on the validation dataset for further analysis. We found that the convergence is more stable if the model is first optimized only on the statistical part of the likelihood Analytic gradient optimizing ω w.r.t. σ(µ) Fig. 2 Graphical overview of the proposed method to optimize the reduction of the dataset used for the statistical inference of the parameters of interest from end to end. The number of observables d in the initial dataset with n observations is reduced to a set of k observables by the neural network function f with the free parameters ω. The dataset is compressed further by summarizing the n observations using a k-dimensional histogram with h bins. Eventually the free parameters ω are optimized with the variance of the parameter of interest µ as objective, which is made possible by an approximated gradient for the histogram. shown in equation 2 and therefore apply the pretraining for 30 gradient steps. We apply statistical weights to scale the expectation of signal and background to 50 and 10 3 , respectively. The dataset is split in half for training and validation, and all results are computed from a statistically independent dataset of the same size as the original one.
The best possible expected result in terms of the variance of the estimate for µ is given by a fit of the unbinned statistical model without dimensionality reduction. Alternatively, we can get an asymptotically close result by using a binned likelihood with sufficiently large number of bins in the two-dimensional input space. The  dataset. Further, we find the uncertainty of µ in all fits by profiling the likelihood [18] rather than using the approximation by the covariance matrix in equation 6. We obtain all results in this paper with validated statistical tools, RooFit and RooStats [19,20,21], such as used by most publications analyzing data of the LHC experiments.
The first comparison to this best-possible result is done by training the NN not on the variance of the estimate for µ, V 00 , but on the cross entropy loss with signal and background weighted to the same expectation. This approach has been used in multiple analyses in highenergy particle physics [22,23]. The NN function f is a sufficient statistic -and therefore optimal -if no systematic uncertainties have to be considered for the statistical inference such as the likelihood in equation 2 [10]. The resulting function f is shown in the input space and by the distribution of the output in figure 5. The NN learns to project the two-dimensional space spanned by x 1 and x 2 on the diagonal, which is trivially the optimal dimensionality reduction in this simple example. If we apply the statistical model including the systematic uncertainty on the histograms in figure 5, the parameter of interest is fitted as µ = 1.0 +0.45 −0.44 with an uncertainty worse by 19 % than the best possible result obtained above.
As a consistency check for our new strategy described in section 2, we train the NN on the variance of the estimate for µ given by V 00 in equation 6 but without adding the nuisance parameter η modelling the systematic uncertainty. The resulting NN function f in the input space, the distribution of the outputs and the profile of the likelihood are shown in figure 6. As expected, the plane of the function f in the input space is qualitatively similar, resulting with µ = 1.0 +0.47 −0.46 in a comparable performance than the training on the cross entropy loss. It should be noted that the systematic uncertainty has been included again for the statistical inference.
When adding the nuisance parameter η to the likelihood, the training of the NN results in the function f shown in figure 7. The uncertainty of the parameter of interest is with the fit result µ = 1.0 +0.39 −0.36 considerably decreased and lowers the residual difference to the optimal result from 19 % to 4 %. The function f in the input space in figure 7 shows that the training identified successfully the signal-enriched region with less contribution of the systematic uncertainty resulting in counts in the histogram yielding high signal statistics with a small uncertainty from the variation of the background process. Figure 7 shows also that the NN function is decorrelated against the systematic uncertainty because the profile of the likelihood changes only little if we remove the systematic uncertainty from the statistical model. The proposed method shares this feature with other approaches for decorrelation of the NN function such as discussed in section 3. The difference is that the strength of the decorrelation is not a hyperparameter but controlled by the higher objective V 00 , which enables us to find directly the best trade-off between statistical and systematic uncertainty contributing to the estimate of µ. The correlation of the parameter of interest µ to the parameter η controlling the systematic variation is reduced from 64 % for the training on the cross entropy loss to 13 % for the training on the variance of the parameter of interest V 00 .
5 Application to a more complex analysis task typical for high-energy particle physics In this section, we apply the proposed method to a problem typical for data analysis in high-energy particle physics at the LHC. We use a subset of the dataset published for the Higgs boson machine learning challenge [24] extended by a systematic variation. The goal of the challenge is to achieve the best possible significance for the signal process representing Higgs boson decays to two tau leptons overlaid by the background simulated as a mixture of different physical processes [24]. We pick from the dataset four variables, namely PRI met, DER mass vis, DER pt h and DER deltaeta jet jet and select only events, which have all of these features defined. In addition to the event weights provided with the dataset, we scale the signal expectation with a factor of two. The final dataset has 244.0 and 35140.1 (106505 and 131480) weighted (unweighted) events for the signal and background process, respectively. The systematic uncertainty in the dataset is assumed as a 10 % uncertainty on the missing transverse energy implemented with the transformation PRI met · (1.0 ± 0.1) and propagated to the other variables using reweighting. The distributions of the variables including the systematic variations are shown in figures 8 to 10. The NN is trained only on three of the four variables, excluding the missing transverse energy. The systematic variations propagated to the remaining variables are thus correlated via a hidden variable, representing a more complex scenario than the simple example in section 4. We split the dataset using one third for training and validation of the NN, and two thirds for the results presented in this paper. The NN architecture and the training procedure are the same as implemented for the simple example in section 4 with the difference that we apply a standardization of the input ranges following the rule (x − x)/σ(x) with the mean x and standard variation σ(x) of the input x.
An (asymptotically) optimal result as derived for the previous example is not available since the likelihood in the input space is not known. Instead we use the training on the cross entropy loss as reference with µ = 1.0 +0.69   For the assessment of the distributions of the NN output, it should be noted that in contrast to the training based on the cross entropy loss, for the training based on V 00 no preference is given for signal (background) events to obtain values close to 1 (0). Similar to the result from the simple example in section 4, the profiles of the likelihood for all scenarios show that the training on V 00 removes the dependence on the systematic uncertainty yielding a smaller variance on µ. On the other hand, the training on the cross entropy optimizes best the estimate of µ in the absence of systematic uncertainties, as expected from our previous discussion. With the proposed strategy, the NN function learns to decorrelate against the systematic uncertainty, visible in the correlation of the signal strength µ to the parameter η controlling the systematic variation, which drops from 69 % for the training on the cross entropy to 4 % for the training on the variance of the parameter of interest V 00 , based on the full likelihood information as given in equation 3.
To improve the estimate of µ for the approach with the NN trained on the cross entropy loss, a possible strategy could be to increase the number of histogram bins to exploit better the separation between the signal and background process. Figure 15 shows the development of the performance with the number of bins for the training on the cross entropy loss and the training on the likelihood via V 00 . The training on the cross entropy loss results in an estimate of µ with a mean correlation to the nuisance parameter η of 66 % and a falling uncertainty in µ with an average distance of 0.18 between the result for taking only the statistical uncertainties and statistical and systematic uncertainties into account for the statistical inference of µ. In contrast, the strategy with the NN trained on V 00 shows a reduction of the correlation between µ and η of 0.35 when moving from two to eight bins for the input histogram for the statistical inference. The estimate remains robust against the systematic variation for all tested configurations, yielding a smaller variance for the estimate of µ compared to the training on the cross entropy loss. The average distance between the inference using only the statistical part of the likelihood and the full statistical model is 0.01. Including the systematic uncertainty in the inference, the comparison of the estimate of µ between the training based on V 00 and the training based on the cross entropy shows an improved variance of µ by 0.07 on average, yielding a stable average improvement of 10 %.
It should be noted that in practice the granularity of the binning is limited by the statistics of data and the simulation. Limited statistical precision in the simulation is usually taken into account by introducing dedicated systematic uncertainties in the statistical model that typically degrade the performance of the analysis for a large number of bins.

Summary
We have presented a novel approach to optimize statistical inference in the presence of systematic uncertainties, when using dimensionality reduction of the dataset and likelihoods based on Poisson statistics. Neural networks and the differential approximation for the gradient of a histogram enables us to optimize directly the variance of the estimate of the parameters of interest in consideration of the nuisance parameters representing the systematic uncertainties of the measurement. The proposed method yields an improved performance for data analysis influenced by systematic uncertainties in comparison to conventional strategies using classification-based objectives for the dimensionality reduction. The improvements are discussed using a simple example based on pseudo-experiments with a known likelihood in the input space and we show that the technique is able to perform a statistical inference close to optimal by leveraging the given information about the systematic uncertainties. The applicability of the method for more complex analyses is demonstrated with an example typical for data analyses in high-energy particle physics.