INFERNO: Inference-Aware Neural Optimisation

Complex computer simulations are commonly required for accurate data modelling in many scientific disciplines, making statistical inference challenging due to the intractability of the likelihood evaluation for the observed data. Furthermore, sometimes one is interested on inference drawn over a subset of the generative model parameters while taking into account model uncertainty or misspecification on the remaining nuisance parameters. In this work, we show how non-linear summary statistics can be constructed by minimising inference-motivated losses via stochastic gradient descent such they provided the smallest uncertainty for the parameters of interest. As a use case, the problem of confidence interval estimation for the mixture coefficient in a multi-dimensional two-component mixture model (i.e. signal vs background) is considered, where the proposed technique clearly outperforms summary statistics based on probabilistic classification, which are a commonly used alternative but do not account for the presence of nuisance parameters.


Introduction
Simulator-based inference is currently at the core of many scientific fields, such as population genetics, epidemiology, and experimental particle physics. In many cases the implicit generative procedure defined in the simulation is stochastic and/or lacks a tractable probability density p(x|θ), where θ ∈ Θ is the vector of model parameters. Given some experimental observations D = {x 0 , ..., x n }, a problem of special relevance for these disciplines is statistical inference on a subset of model parameters ω ∈ Ω ⊆ Θ. This can be approached via likelihood-free inference algorithms such as Approximate Bayesian Computation (ABC) [1], simplified synthetic likelihoods [2] or density estimation-by-comparison approaches [3].
Because the relation between the parameters of the model and the data is only available via forward simulation, most likelihood-free inference algorithms tend to be computationally expensive due to the need of repeated simulations to cover the parameter space. When data are high-dimensional, likelihood-free inference can rapidly become inefficient, so low-dimensional summary statistics s(D) are used instead of the raw data for tractability. The choice of summary statistics for such cases becomes critical, given that naive choices might cause loss of relevant information and a corresponding degradation of the power of resulting statistical inference.
As a motivating example we consider data analyses at the Large Hadron Collider (LHC), such as those carried out to establish the discovery of the Higgs boson [4,5]. In that framework, the ultimate aim is to extract information about Nature from the large amounts of high-dimensional data on the subatomic particles produced by energetic collision of protons, and acquired by highly complex detectors built around the collision point. Accurate data modelling is only available via stochastic simulation of a complicated chain of physical processes, from the underlying fundamental interaction to the subsequent particle interactions with the detector elements and their readout. As a result, the density p(x|θ) cannot be analytically computed. The inference problem in particle physics is commonly posed as hypothesis testing based on the acquired data. An alternate hypothesis H 1 (e.g. a new theory that predicts the existence of a new fundamental particle) is tested against a null hypothesis H 0 (e.g. an existing theory, which explains previous observed phenomena). The aim is to check whether the null hypothesis can be rejected in favour of the alternate hypothesis at a certain confidence level surpassing 1 − α, where α, known as the Type I error rate, is commonly set to α = 3 × 10 −7 for discovery claims. Because α is fixed, the sensitivity of an analysis is determined by the power 1 − β of the test, where β is the probability of rejecting a false null hypothesis, also known as Type II error rate.
Due to the high dimensionality of the observed data, a low-dimensional summary statistic has to be constructed in order to perform inference. A well-known result of classical statistics, the Neyman-Pearson lemma [6], establishes that the likelihood-ratio Λ(x) = p(x|H 0 )/p(x|H 1 ) is the most powerful test when two simple hypotheses are considered. As p(x|H 0 ) and p(x|H 1 ) are not available, simulated samples are used in practice to obtain an approximation of the likelihood ratio by casting the problem as supervised learning classification.
In many cases, the nature of the generative model (a mixture of different processes) allows the treatment of the problem as signal (S) vs background (B) classification [7], when the task becomes one of effectively estimating an approximation of p S (x)/p B (x) which will vary monotonically with the likelihood ratio. While the use of classifiers to learn a summary statistic can be effective and increase the discovery sensitivity, the simulations used to generate the samples which are needed to train the classifier often depend on additional uncertain parameters (commonly referred as nuisance parameters). These nuisance parameters are not of immediate interest but have to be accounted for in order to make quantitative statements about the model parameters based on the available data. Classification-based summary statistics cannot easily account for those effects, so their inference power is degraded when nuisance parameters are finally taken into account.
In this work, we present a new machine learning method to construct non-linear sample summary statistics that directly optimises the expected amount of information about the subset of parameters of interest using simulated samples, by explicitly and directly taking into account the effect of nuisance parameters. In addition, the learned summary statistics can be used to build synthetic sample-based likelihoods and perform robust and efficient classical or Bayesian inference from the observed data, so they can be readily applied in place of current classification-based or domain-motivated summary statistics in current scientific data analysis workflows.

Problem Statement
Let us consider a set of n i.i.d. observations D = {x 0 , ..., x n } where x ∈ X ⊆ R d , and a generative model which implicitly defines a probability density p(x|θ) used to model the data. The generative model is a function of the vector of parameters θ ∈ Θ ⊆ R p , which includes both relevant and nuisance parameters. We want to learn a function s : D ⊆ R d×n → S ⊆ R b that computes a summary statistic of the dataset and reduces its dimensionality so likelihood-free inference methods can be applied effectively. From here onwards, b will be used to denote the dimensionality of the summary statistic s(D).
While there might be infinite ways to construct a summary statistic s(D), we are only interested in those that are informative about the subset of interest ω ∈ Ω ⊆ Θ of the model parameters. The concept of statistical sufficiency is especially useful to evaluate whether summary statistics are informative. In the absence of nuisance parameters, classical sufficiency can be characterised by means of the factorisation criterion: where h and g are non-negative functions. If p(D|ω) can be factorised as indicated, the summary statistic s(D) will yield the same inference about the parameters ω as the full set of observations D. When nuisance parameters have to be accounted in the inference procedure, alternate notions of sufficiency are commonly used such as partial or marginal sufficiency [8,9]. Nonetheless, for the problems of relevance in this work, the probability density is not available in closed form so the general task of finding a sufficient summary statistic cannot be tackled directly. Hence, alternative methods to build summary statistics have to be followed.
For simplicity, let us consider a problem where we are only interested on statistical inference on a single one-dimensional model parameter ω = {ω 0 } given some observed data. Be given a summary statistic s and a statistical procedure to obtain an unbiased interval estimate of the parameter of interest which accounts for the effect of nuisance parameters. The resulting interval can be characterised by its width ∆ω 0 =ω + 0 −ω − 0 , defined by some criterion so as to contain on average, upon repeated samping, a given fraction of the probability density, e.g. a central 68.3% interval. The expected size of the interval depends on the summary statistic s chosen: in general, summary statistics that are more informative about the parameters of interest will provide narrower confidence or credible intervals on their value. Under this figure of merit, the problem of choosing an optimal summary statistic can be formally expressed as finding a summary statistic s * that minimises the interval width: s * = argmin s ∆ω 0 .
(2) The above construction can be extended to several parameters of interest by considering the interval volume or any other function of the resulting confidence or credible regions.

Method
In this section, a machine learning technique to learn non-linear sample summary statistics is described in detail. The method seeks to minimise the expected variance of the parameters of interest obtained via a non-parametric simulation-based synthetic likelihood. A graphical description of the technique is depicted on Fig. 1. The parameters of a neural network are optimised by stochastic gradient descent within an automatic differentiation framework, where the considered loss function accounts for the details of the statistical model as well as the expected effect of nuisance parameters.  The family of summary statistics s(D) considered in this work is composed by a neural network model applied to each dataset observation f (x; φ) : X ⊆ R d → Y ⊆ R b , whose parameters φ will be learned during training by means of stochastic gradient descent, as will be discussed later. Therefore, using set-builder notation the family of summary statistics considered can be denoted as: will reduce the dimensionality from the input observations space X to a lowerdimensional space Y. The next step is to map observation outputs to a dataset summary statistic, which will in turn be calibrated and optimised via a non-parametric likelihood L(D; θ, φ) created using a set of simulated observations G s = {x 0 , ..., x g }, generated at a certain instantiation of the simulator parameters θ s .
In experimental high energy physics experiments, which are the scientific context that initially motivated this work, histograms of observation counts are the most commonly used non-parametric density estimator because the resulting likelihoods can be expressed as the product of Poisson factors, one for each of the considered bins. A naive sample summary statistic can be built from the output of the neural network by simply assigning each observation x to a bin corresponding to the cardinality of the maximum element of f (x; φ), so each element of the sample summary will correspond to the following sum: which can in turn be used to build the following likelihood, where the expectation for each bin is taken from the simulated sample G s : where the n/g factor accounts for the different number of observations in the simulated samples. In cases where the number of observations is itself a random variable providing information about the parameters of interest, or where the simulated observation are weighted, the choice of normalisation of L may be slightly more involved and problem specific, but nevertheless amenable.
In the above construction, the chosen family of summary statistics is non-differentiable due to the argmax operator, so gradient-based updates for the parameters cannot be computed. To work around this problem, a differentiable approximationŝ(D; φ) is considered. This function is defined by means of a sof tmax operator:ŝ where the temperature hyper-parameter τ will regulate the softness of the operator. In the limit of τ → 0 + , the probability of the largest component will tend to 1 while others to 0, and thereforê s(D; φ) → s(D; φ). Similarly, let us denote byL(D; θ, φ) the differentiable approximation of the non-parametric likelihood obtained by substituting s(D; φ) withŝ(D; φ). Instead of using the observed data D, the value ofL may be computed when the observation for each bin is equal to its corresponding expectation based on the simulated sample G s , which is commonly denoted as the Asimov likelihood [10]L A : for which it can be easily proven that argmax θ∈θ (L A (θ; φ)) = θ s , so the maximum likelihood estimator (MLE) for the Asimov likelihood is the parameter vector θ s used to generate the simulated dataset G s . In Bayesian terms, if the prior over the parameters is flat in the chosen metric, then θ s is also the maximum a posteriori (MAP) estimator. By taking the negative logarithm and expanding in θ around θ s , we can obtain the Fisher information matrix [11] for the Asimov likelihood: which can be computed via automatic differentiation if the simulation is differentiable and included in the computation graph or if the effect of varying θ over the simulated dataset G s can be effectively approximated. While this requirement does constrain the applicability of the proposed technique to a subset of likelihood-free inference problems, it is quite common for e.g. physical sciences that the effect of the parameters of interest and the main nuisance parameters over a sample can be approximated by the changes of mixture coefficients of mixture models, translations of a subset of features, or conditional density ratio re-weighting.
Ifθ is an unbiased estimator of the values of θ, the covariance matrix fulfils the Cramér-Rao lower bound [12,13]: cov θ (θ) ≥ I(θ) −1 (9) and the inverse of the Fisher information can be used as an approximate estimator of the expected variance, given that the bound would become an equality in the asymptotic limit for MLE. If some of the parameters θ are constrained by independent measurements characterised by their likelihoods {L 0 C (θ), ..., L c C (θ)}, those constraints can also be easily included in the covariance estimation, simply by considering the augmented likelihoodL A instead ofL A in Eq. 8: In Bayesian terminology, this approach is referred to as the Laplace approximation [14] where the logarithm of the joint density (including the priors) is expanded around the MAP to a multidimensional normal approximation of the posterior density: which has already been approached by automatic differentiation in probabilistic programming frameworks [15]. While a histogram has been used to construct a Poisson count sample likelihood, non-parametric density estimation techniques can be used in its place to construct a product of observation likelihoods based on the neural network output f (x; φ) instead. For example, an extension of this technique to use kernel density estimation (KDE) should be straightforward, given its intrinsic differentiability.
The loss function used for stochastic optimisation of the neural network parameters φ can be any function of the inverse of the Fisher information matrix at θ s , depending on the ultimate inference aim. The diagonal elements I −1 ii (θ s ) correspond to the expected variance of each of the φ i under the normal approximation mentioned before, so if the aim is efficient inference about one of the parameters ω 0 = θ k a candidate loss function is: which corresponds to the expected width of the confidence interval for ω 0 accounting also for the effect of the other nuisance parameters in θ. This approach can also be extended when the goal is inference over several parameters of interest ω ⊆ θ (e.g. when considering a weighted sum of the relevant variances). A simple version of the approach just described to learn a neural-network based summary statistic employing an inference-aware loss is summarised in Algorithm 1.
Sample a representative mini-batch G s from g(θ s ).

7:
Update network parameters φ → SGD(∇ φ U ). 8: end for 4 Related Work Classification or regression models have been implicitly used to construct summary statistics for inference in several scientific disciplines. For example, in experimental particle physics, the mixture model structure of the problem makes it amenable to supervised classification based on simulated datasets [16,17]. While a classification objective can be used to learn powerful feature representations and increase the sensitivity of an analysis, it does not take into account the details of the inference procedure or the effect of nuisance parameters like the solution proposed in this work.
The first known effort to include the effect of nuisance parameters in classification and explain the relation between classification and the likelihood ratio was by Neal [18]. In the mentioned work, Neal proposes training of classifier including a function of nuisance parameter as additional input together with a per-observation regression model of the expectation value for inference. Cranmer et al. [3] improved on this concept by using a parametrised classifier to approximate the likelihood ratio which is then calibrated to perform statistical inference. At variance with the mentioned works, we do not consider a classification objective at all and the neural network is directly optimised based on an inference-aware loss. Additionally, once the summary statistic has been learnt the likelihood can be trivially constructed and used for classical or Bayesian inference without a dedicated calibration step. Furthermore, the approach presented in this work can also be extended, as done by Baldi et al. [19] by a subset of the inference parameters to obtain a parametrised family of summary statistics with a single model.
Recently, Brehmer et al. [20,21,22] further extended the approach of parametrised classifiers to better exploit the latent-space space structure of generative models from complex scientific simulators. Additionally they propose a family of approaches that include a direct regression of the likelihood ratio and/or likelihood score in the training losses. While extremely promising, the most performing solutions are designed for a subset of the inference problems at the LHC and they require considerable changes in the way the inference is carried out. The aim of this work is different, as we try to learn sample summary statistics that may act as a plug-in replacement of classifier-based dimensionality reduction and can be applied to general likelihood-free problems where the effect of the parameters can be modelled or approximated.
Within the field of Approximate Bayesian Computation (ABC), there have been some attempts to use neural network as a dimensionality reduction step to generate summary statistics. For example, Jiang et al. [23] successfully employ a summary statistic by directly regressing the parameters of interest and therefore approximating the posterior mean given the data, which then can be used directly as a summary statistic.
A different path is taken by Louppe et al. [24], where the authors present a adversarial training procedure to enforce a pivotal property on a predictive model. The main concern of this approach is that a classifier which is pivotal with respect to nuisance parameters might not be optimal, neither for classification nor for statistical inference. Instead of aiming for being pivotal, the summary statistics learnt by our algorithm attempt to find a transformation that directly reduces the expected effect of nuisance parameters over the parameters of interest.

Experiments
In this section, we first study the effectiveness of the inference-aware optimisation in a synthetic mixture problem where the likelihood is known. We then compare our results with those obtained by standard classification-based summary statistics. All the code needed to reproduce the results presented the results presented here is available in an online repository [25], extensively using TENSORFLOW [26] and TENSORFLOW PROBABILITY [15,27] software libraries.

3D Synthetic Mixture
In order to exemplify the usage of the proposed approach, evaluate its viability and test its performance by comparing to the use of a classification model proxy, a three-dimensional mixture example with two components is considered. One component will be referred as background f b (x|λ) and the other as signal f s (x); their probability density functions are taken to correspond respectively to: so that (x 0 , x 1 ) are distributed according to a multivariate normal distribution while x 2 follows an independent exponential distribution both for background and signal, as shown in Fig. 2a. The signal distribution is fully specified while the background distribution depends on r, a parameter which shifts the mean of the background density, and a parameter λ which specifies the exponential rate in the third dimension. These parameters will be the treated as nuisance parameters when benchmarking different methods. Hence, the probability density function of observations has the following mixture structure: where µ is the parameter corresponding to the mixture weight for the signal and consequently (1 − µ) is the mixture weight for the background. The low-dimensional projections from samples from the mixture distribution for a small µ = 50/1050 is shown in Fig. 2b. In this toy problem, we consider a case where the underlying model predicts that the total number of observations are Poisson distributed with a mean s + b, where s and b are the expected number of signal and background observations. Thus the following parametrisation will be more convenient for building sample-based likelihoods: This parametrisation is common for physics analyses at the LHC, because theoretical calculations provide information about the expected number of observations. If the probability density is known, but the expectation for the number of observed events depends on the model parameters, the likelihood can be extended [28] with a Poisson count term as: which will be used to provide an optimal inference baseline when benchmarking the different approaches. Another quantity of relevance is the conditional density ratio, which would correspond to the optimal classifier (in the Bayes risk sense) separating signal and background events in a balanced dataset (equal priors): (18) noting that this quantity depends on the parameters that define the background distribution r and λ, but not on s or b that are a function of the mixture coefficients. It can be proven (see appendix A ) that s * (x) is a sufficient summary statistic with respect to an arbitrary two-component mixture model if the only unknown parameter is the signal mixture fraction µ (or alternatively s in the chosen parametrisation). In practise, the probability density functions of signal and background are not known analytically, and only forward samples are available through simulation, so alternative approaches are required.
While the synthetic nature of this example allows to rapidly generate training data on demand, a training dataset of 200,000 simulated observations has been considered, in order to study how the proposed method performs when training data is limited. Half of the simulated observations correspond to the signal component and half to the background component. The latter has been generated using r = 0.0 and λ = 3.0. A validation holdout from the training dataset of 200,000 observations is only used for computing relevant metrics during training and to control over-fitting. The final figures of merit that allow to compare different approaches are computed using a larger dataset of 1,000,000 observations. For simplicity, mini-batches for each training step are balanced so the same number of events from each component is taken both when using standard classification or inference-aware losses.
An option is to pose the problem as one of classification based on a simulated dataset. A supervised machine learning model such a neural network can be trained to discriminate signal and background observations, considering a fixed parameters r and λ. The output of such a model typically consist in class probabilities c s and c b given an observation x, which will tend asymptotically to the optimal classifier from Eq. 18 given enough data, a flexible enough model and a powerful learning rule. The conditional class probabilities (or alternatively the likelihood ratio f s (x)/f b (x)) are powerful learned features that can be used as summary statistic; however their construction ignores the effect of the nuisance parameters r and λ on the background distribution. Furthermore, some kind of nonparametric density estimation (e.g. a histogram) has to be considered in order to build a calibrated statistical model using the classification-based learned features, which will in turn smooth and reduce the information available for inference.
To exemplify the use of this family of classification-based summary statistics, a histogram of a deep neural network classifier output trained on simulated data and its variation computed for different values of r and λ are shown in Fig. 3a. The details of the training procedure will be provided later in this document. The classifier output can be directly compared with s(x|r = 0.0, λ = 3.0) evaluated using the analytical distribution function of signal and background according to Eq. 18, which is shown in Fig. 3b and corresponds to the optimal classifier. The trained classifier approximates very well the optimal classifier. The summary statistic distribution for background depends considerably on the value of the nuisance parameters both for the trained and the optimal classifier, which will in turn cause an important degradation on the subsequent statistical inference. The statistical model described above has up to four unknown parameters: the expected number of signal observations s, the background mean shift r, the background exponential rate in the third dimension λ, and the expected number of background observations. The effect of the expected number of signal and background observations s and b can be easily included in the computation graph by weighting the signal and background observations. This is equivalent to scaling the resulting vector of Poisson counts (or its differentiable approximation) if a non-parametric counting model as the one described in Sec. 3 is used. Instead the effect of r and λ, both nuisance parameters that will define the background distribution, is more easily modelled as a transformation of the input data x. In particular, r is a nuisance parameter that causes a shift on the background along the first dimension and its effect can accounted for in the computation graph by simply adding (r, 0.0, 0.0) to each observation in the mini-batch generated from the background distribution. Similarly, the effect of λ can be modelled by multiplying x 2 by the ratio between the λ 0 used for generation and the one being modelled. These transformations are specific for this example, but alternative transformations depending on parameters could also be accounted for as long as they are differentiable or substituted by a differentiable approximation.
For this problem, we are interested in carrying out statistical inference on the parameter of interest s. In fact, the performance of inference-aware optimisation as described in Sec. 3 will be compared with classification-based summary statistics for a series of inference benchmarks based on the synthetic problem described above that vary in the number of nuisance parameters considered and their constraints: • Benchmark 0: no nuisance parameters are considered, both signal and background distributions are taken as fully specified (r = 0.0, λ = 3.0 and b = 1000.).
• Benchmark 1: r is considered as an unconstrained nuisance parameter, while λ = 3.0 and b = 1000 are fixed.
• Benchmark 2: r and λ are considered as unconstrained nuisance parameters, while b = 1000 is fixed.
When using classification-based summary statistics, the construction of a summary statistic does depend on the presence of nuisance parameters, so the same model is trained independently of the benchmark considered. In real-world inference scenarios, nuisance parameters have often to be accounted for and typically are constrained by prior information or auxiliary measurements. For the approach presented in this work, inference-aware neural optimisation, the effect of the nuisance parameters and their constraints can be taken into account during training. Hence, 5 different training procedures for INFERNO will be considered, one for each of the benchmarks, denoted by the same number.
The same basic network architecture is used both for cross-entropy and inference-aware training: two hidden layers of 100 nodes followed by ReLU activations. The number of nodes on the output layer is two when classification proxies are used, matching the number of mixture classes in the problem considered. Instead, for inference-aware classification the number of output nodes can be arbitrary and will be denoted with b, corresponding to the dimensionality of the sample summary statistics. The final layer is followed by a softmax activation function and a temperature τ = 0.1 for inference-aware learning to ensure that the differentiable approximations are closer to the true expectations. Standard mini-batch stochastic gradient descent (SGD) is used for training and the optimal learning rate is fixed and decided by means of a simple scan; the best choice found is specified together with the results.
In Fig. 4a, the dynamics of inference-aware optimisation are shown by the validation loss, which corresponds to the approximate expected variance of parameter s, as a function of the training step for 10 random-initialised instances of the INFERNO model corresponding to Benchmark 2. All inference-aware models were trained during 200 epochs with SGD using mini-batches of 2000 observations and a learning rate γ = 10 −6 . All the model initialisations converge to summary statistics that provide low variance for the estimator of s when the nuisance parameters are accounted for.
To compare with alternative approaches and verify the validity of the results, the profiled likelihoods obtained for each model are shown in Fig. 4b. The expected uncertainty if the trained models are used for subsequent inference on the value of s can be estimated from the profile width when ∆L = 0.5. Hence, the average width for the profile likelihood using inference-aware training, 16.97 ± 0.11, can  be compared with the corresponding one obtained by uniformly binning the output of classificationbased models in 10 bins, 24.01 ± 0.36. The models based on cross-entropy loss were trained during 200 epochs using a mini-batch size of 64 and a fixed learning rate of γ = 0.001.
A more complete study of the improvement provided by the different INFERNO training procedures is provided in Table 1, where the median and 1-sigma percentiles on the expected uncertainty on s are provided for 100 random-initialised instances of each model. In addition, results for 100 randominitialised cross-entropy trained models and the optimal classifier and likelihood-based inference are also included for comparison. The confidence intervals obtained using INFERNO-based summary statistics are considerably narrower than those using classification and tend to be much closer to those expected when using the true model likelihood for inference. Much smaller fluctuations between initialisations are also observed for the INFERNO-based cases. The improvement over classification increases when more nuisance parameters are considered. The results also seem to suggest the inclusion of additional information about the inference problem in the INFERNO technique leads to comparable or better results than its omission.
Given that a certain value of the parameters θ s has been used to learn the summary statistics as described in Algorithm 1 while their true value is unknown, the expected uncertainty on s has also been computed for cases when the true value of the parameters θ true differs. The variation of the expected uncertainty on s when either r or λ is varied for classification and inference-aware summary  statistics is shown in Fig. 5 for Benchmark 2. The inference-aware summary statistics learnt for θ s work well when θ true = θ s in the range of variation explored.
This synthetic example demonstrates that the direct optimisation of inference-aware losses as those described in the Sec. 3 is effective. The summary statistics learnt accounting for the effect of nuisance parameters compare very favourably to those obtained by using a classification proxy to approximate the likelihood ratio. Of course, more experiments are needed to benchmark the usefulness of this technique for real-world inference problems as those found in High Energy Physics analyses at the LHC.

Conclusions
Classification-based summary statistics for mixture models often suffer from the need of specifying a fixed model of the data, thus neglecting the effect of nuisance parameters in the learning process. The effect of nuisance parameters is only considered downstream of the learning phase, resulting in sub-optimal inference on the parameters of interest.
In this work we have described a new approach for building non-linear summary statistics for likelihood-free inference that directly minimises the expected variance of the parameters of interest, which is considerably more effective than the use of classification surrogates when nuisance parameters are present.
The results obtained for the synthetic experiment considered clearly demonstrate that machine learning techniques, in particular neural networks, can be adapted for learning summary statistics that match the particularities of the inference problem at hand, greatly increasing the information available for subsequent inference. The application of INFERNO to non-synthetic examples where nuisance parameters are relevant, such as the systematic-extended Higgs dataset [29], are left for future studies.
Furthermore, the technique presented can be applied to arbitrary likelihood-free problems as long as the effect of parameters over the simulated data can be implemented as a differentiable transformations. As a possible extension, alternative non-parametric density estimation techniques such as kernel density could very well be used in place Poisson count models. like to acknowledge Gilles Louppe and Joeri Hermans for some useful discussions directly related to this work. This work is part of a more general effort to develop new statistical and machine learning techniques to use in High Energy Physics analyses within within the AMVA4NewPhysics project, which is supported by the European Union's Horizon 2020 research and innovation programme under Grant Agreement number 675440. CloudVeneto is also acknowledged for the use of computing and storage facilities provided.