Learning Likelihood Ratios with Neural Network Classifiers
The likelihood ratio is a crucial quantity for statistical inference in
science that enables hypothesis testing, construction of confidence intervals,
reweighting of distributions, and more. Many modern scientific applications,
however, make use of data- or simulation-driven models for which computing the
likelihood ratio can be very difficult or even impossible. By applying the
so-called ``likelihood ratio trick,'' approximations of the likelihood ratio
may be computed using clever parametrizations of neural network-based
classifiers. A number of different neural network setups can be defined to
satisfy this procedure, each with varying performance in approximating the
likelihood ratio when using finite training data. We present a series of
empirical studies detailing the performance of several common loss functionals
and parametrizations of the classifier output in approximating the likelihood
ratio of two univariate and multivariate Gaussian distributions as well as
simulated high-energy particle physics datasets.
Shahzar Rizvi, Mariel Pettee, Benjamin Nachman
# Learning Likelihood Ratios with Neural Network Classifiers

###### Abstract

The likelihood ratio is a crucial quantity for statistical inference in science that enables hypothesis testing, construction of confidence intervals, reweighting of distributions, and more. Many modern scientific applications, however, make use of data- or simulation-driven models for which computing the likelihood ratio can be very difficult or even impossible. By applying the so-called "likelihood ratio trick," approximations of the likelihood ratio may be computed using clever parametrizations of neural network-based classifiers. A number of different neural network setups can be defined to satisfy this procedure, each with varying performance in approximating the likelihood ratio when using finite training data. We present a series of empirical studies detailing the performance of several common loss functionals and parametrizations of the classifier output in approximating the likelihood ratio of two univariate and multivariate Gaussian distributions as well as simulated high-energy particle physics datasets.

## I Introduction

Claiming a scientific discovery requires a hypothesis test, i.e. a statistical threshold for claiming that one's experimental data reject the null hypothesis in favor of an alternative hypothesis. This might involve two probability densities:

* \(H_{0}\) (the null hypothesis)
* \(H_{1}\) (the alternative hypothesis)

By the Neyman-Pearson lemma [1], the strongest ("uniformly most powerful") measure of whether the experimental data \(x\) support \(H_{0}\) vs. \(H_{1}\) is a likelihood ratio test. These tests are particularly widespread in reporting results in High-Energy Physics (HEP), but are also commonly used for statistical analyses across astrophysics, biology, medicine, and other scientific domains concerned with hypothesis testing or confidence intervals. The need for likelihood ratios goes beyond hypothesis testing, too - they can also be used to reweight a distribution to align with a target distribution, such as reweighting simulation samples to match real data [2; 3; 4; 5; 6; 7; 8; 9].

In the simplest form of a likelihood ratio test, where \(H_{0}\) and \(H_{1}\) are fully-defined by parameters \(\theta_{0}\) and \(\theta_{1}\), the background-only hypothesis is either rejected (or not) depending on the value of the ratio of likelihoods \(p(\theta_{0}\mid x)\) under \(H_{0}\) and \(p(\theta_{1}\mid x)\) under \(H_{1}\) in relation to the desired significance level.

In practice, however, the probability densities \(H_{0}\) and \(H_{1}\) may not be explicitly known. Worse, they might be nearly impossible to compute, such as in instances where they are generated by a complex simulation model. In these cases, we can use machine learning to directly approximate the likelihood ratio itself, bypassing the need to approximate the individual probability densities.

A classifier function \(f(x)\) (for instance, from a neural network) designed to distinguish data sampled from \(H_{0}\) (\(f(x)\to 0\)) vs. \(H_{1}\) (\(f(x)\to 1\)) can be used to approximate the likelihood ratio by minimizing a proper loss functional (defined in Section II):

\[\operatorname*{argmin}_{f}L[f]=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{1})}= \mathcal{L}(x). \tag{1}\]

For instance, in the familiar case of training a classifier by minimizing the binary cross-entropy loss (see I), the optimal decision function \(f(x)\) is:

\[f(x)=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{0})+p(x\mid\theta_{1})}. \tag{2}\]

We can then approximate the likelihood ratio with a monotonic transformation of the neural network output \(f(x)\)1:

Footnote 1: This notation assumes balanced training sets for simplicity. With imbalanced classes, one would need to modify the likelihood ratio to include prior factors \(p(\theta_{i})\), though the likelihood ratio trick will still apply [10].

\[\frac{f(x)}{1-f(x)} =\frac{\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{0})+p(x\mid\theta _{1})}}{1-\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{0})+p(x\mid\theta_{1})}} \tag{3}\] \[=\frac{p(x\mid\theta_{0})}{\underline{p(x\!+\!\theta_{0})}+p(x \mid\theta_{1})-\underline{p(x\!+\!\theta_{0})}}\] (4) \[=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{1})}=\mathcal{L}(x). \tag{5}\]
This procedure, sometimes called the "likelihood ratio trick", is well-known in statistics (see e.g. [11; 12; 13]) and has been frequently used in particle physics [14; 15; 16; 2; 6; 10; 26].

A number of different loss functionals beyond binary cross-entropy can be defined to satisfy this setup, but in practice, not all such classifiers will perform equally well when approximating the likelihood ratio. In this paper, we perform a series of empirical studies to understand how different choices of loss functional and parametrization of the resulting classifier affect the performance of likelihood ratio approximation for pairs of distributions.

Several recent works have investigated some improved configurations for the likelihood ratio trick in certain scientific contexts. [27] introduces a new likelihood estimation procedure as an extension of [14] using binary cross-entropy loss with SELU [28] activation. [29] notes that for one- and two-dimensional toy simulations of particle physics datasets, the maximum likelihood classifier (MLC) loss performed better than the binary cross-entropy loss when estimating the likelihood ratio - the first application of MLC loss in particle physics. [10] directly compares linear and exponential parameterizations of maximum likelihood classifier loss with binary cross-entropy loss for one-dimensional Gaussians. [14] uses calibrated classifiers to improve likelihood ratio estimation, and [18; 19] define several different approaches to likelihood ratio estimation, including augmenting the likelihood ratio trick with score regression (Rascal, Sally, etc.). [30] introduces modified versions of the cross-entropy loss that show stronger performance under limited training dataset sizes than the typical cross-entropy loss, while [31] compares the estimation of the likelihood ratio via mean square loss with ELU [32] activation, cross-entropy loss with sigmoid activation, and a proposed exponential loss with no activation function on univariate Gaussian distributions. Still other methods use normalizing flows to determine the likelihood ratio by modeling the individual densities [33; 34] or to obviate the need for the likelihood ratio approximation for reweighting distributions [35].

In light of these existing studies, this work serves as a detailed comparison of a wide range of configurations of loss functionals and output parametrizations across datasets including one-dimensional Gaussians, multi-dimensional Gaussians, and simulated high-energy particle physics datasets. We aim to highlight some best practices and serve as a guide for approximating likelihood ratios with neural network classifiers in the wider scientific community, and particularly within the domains of particle physics and astrophysics.

This paper is organized as follows. In Section II, we summarize the theoretical foundation for learning likelihood ratios with neural network classifiers. In Section III, we present a series of studies focused on optimizing likelihood ratio estimation for one-dimensional Gaussian distributions where the true likelihood ratio is exactly known. In Section IV, we extend these studies to multi-dimensional Gaussian distributions. In Section V, we present some more realistic examples using simulated high-energy physics data where the true likelihood ratio is approximated using a Normalizing Flow model [33]. Finally, we summarize our conclusions and recommendations for further studies in Section VI.

## II Learning likelihood ratios

Let the parameters \(\theta_{0}\) and \(\theta_{1}\) define two distributions, \(p(x\mid\theta_{0})\) and \(p(x\mid\theta_{1})\), as described in Section II.A. of [10]. The goal is to determine or approximate the likelihood ratio

\[\mathcal{L}(x)=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{1})} \tag{6}\]

between the two distributions.

Consider the general loss functional that depends on a learnable function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}\) and rescaling functions \(A:\mathbb{R}\rightarrow\mathbb{R}\) and \(B:\mathbb{R}\rightarrow\mathbb{R}\):

\[\begin{split} L[f]=-\int\mathrm{d}x\bigg{(}& p(x\mid\theta_{0})A(f(x))\\ &+p(x\mid\theta_{1})B(f(x))\bigg{)}.\end{split} \tag{7}\]

We can take the functional derivative of the loss functional to show that the extremum can be transformed to obtain the likelihood ratio:

\[\frac{\delta L}{\delta f} =-\frac{\partial}{\partial f}\Big{(}p(x\mid\theta_{0})A(f(x))+p (x\mid\theta_{1})B(f(x))\Big{)} \tag{8}\] \[=-\bigg{(}p(x\mid\theta_{0})A^{\prime}(f(x))\cdot f^{\prime}(x)\] (9) \[\qquad+p(x\mid\theta_{0})B^{\prime}(f(x))\cdot f^{\prime}(x)\bigg{)}\] \[=0\iff-\frac{B^{\prime}(f(x))}{A^{\prime}(f(x))}=\frac{p(x\mid \theta_{0})}{p(x\mid\theta_{1})}=\mathcal{L}(x). \tag{10}\]

Given that \(-B^{\prime}(f)/A^{\prime}(f)\) is a monotonic rescaling of \(f\) and \(L[f]\) is convex, the learned function \(f\) is an optimal classifier.

In this paper, we first consider the four loss functionals defined by the rescaling functions in Table 1. While this is by no means an exhaustive list of all possible loss functionals, it includes a diverse array of different loss configurations. As detailed in Sec. III.3, we also consider generalized forms of two of these four loss functionals.

A neural network parametrizes the learned function \(f\) as \(\phi(z)\), where \(z\) is the pre-activation output of the network and \(\phi\) is the final activation function. For the binary cross entropy (BCE) and mean squared error (MSE) losses,

\[\mathcal{L}(x)=-\frac{B^{\prime}(f)}{A^{\prime}(f)}=\frac{f}{1-f}, \tag{11}\]
so the likelihood ratio is the odds ratio of the learned function. That is, minimizing the BCE and MSE losses defines a classifier that computes

\[\operatorname*{argmin}_{f}L[f]=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{0})+p(x \mid\theta_{1})}\in(0,1) \tag{12}\]

. To parametrize \(f\) such that the likelihood ratio is non-negative, we require that \(\phi:\mathbb{R}\rightarrow(0,1)\).

However, for the maximum likelihood classifier (MLC) and square root (SQR) losses,

\[\mathcal{L}(x)=-\frac{B^{\prime}(f)}{A^{\prime}(f)}=f, \tag{13}\]

so the likelihood ratio is the learned function, without transformation.

\[\operatorname*{argmin}_{f}L[f]=\frac{p(x\mid\theta_{0})}{p(x\mid\theta_{1})} \tag{14}\]

In this case the loss-minimizing classifier computes the likelihood ratio \(\mathcal{L}(x)\in(0,\infty)\). The requirement on \(\phi\) is that \(\phi:\mathbb{R}\rightarrow(0,\infty)\).

## III Univariate Gaussians

In our first case study, we consider two Gaussian distributions with slightly different means and unit variances: \(X_{0}\sim\text{Normal}(+0.1,1)\) and \(X_{1}\sim\text{Normal}(-0.1,1)\). We also considered univariate Beta and Gamma distributions - these results can be found in Appendix A.

While one could in principle use Boosted Decision Trees (BDTs) instead of neural networks for the classifiers, we found that neural networks outperformed BDTs across a variety of test cases, as shown in Appendix C. All of our classifiers are therefore implemented as neural networks using Keras [36] with a Tensorflow [37] backend and Adam[38] optimizer. Each classifier consists of three hidden layers with 64, 128, and 64 nodes, sequentially. Rectified Linear Unit (ReLU) activation functions are used for the intermediate layers, with the activation for the output layer depending on the loss used to train the neural network and the parametrization being tested. Each of the three hidden layers is followed by a dropout layer with dropout probability of 10%.

Unless otherwise stated, the networks were trained with 1,000,000 samples (750,000 used for training and 250,000 used for validation). 100,000 separate samples were used to evaluate the networks' performances (in particular, to calculate their mean absolute errors). Each network was trained for up to 100 epochs with a batch size of 10%, as in [10]. If the validation loss did not decrease for 10 consecutive epochs, the training was stopped (early stopping with a patience of 10). No detailed hyperparameter optimization was done.

### Naive Implementation

#### iii.1.1 Motivation

The naive parametrization for \(\phi(z)\) in the case of the BCE and MSE losses is \(\phi=\sigma\), the logistic function commonly used as the activation for classification tasks. In the case of the MLC and SQR losses, the most common parametrization would be \(\phi=\text{ReLU}\), the rectified linear unit activation. We chose these parametrizations for our naive implementation.

To better understand how these common parametrizations of the classifiers affect their ability to learn the likelihood ratio, we implemented neural network architecture with each of the four losses, trained them to classify between the two Gaussian distributions. Since the true likelihood ratio is known, we can compare how well each of the four classifiers learns the likelihood ratio function.

#### iii.1.2 Methods

We implemented each classifier using an identical neural network architecture, differing only in the final activation, which acted as either the logistic (for the BCE and MSE classifiers) or ReLU (for the MLC and SQR classifiers) parametrizations for the learned function.

We then trained each of the four classifier architecture on the dataset 100 times each, using the classifier's corresponding loss functional. Each classifier was evaluated on the interval \((-6,6)\) and transformed into the likelihood ratio over that same interval using the appropriate transformation from equations 11 and 13. We averaged the resulting 100 predictions for the likelihood ratio.

To numerically compare the performances of different classifiers in learning the likelihood ratio, we computed their empirical mean absolute errors over 100,000 samples. For \(\hat{\mathcal{L}}\) the estimated likelihood ratio, the mean absolute error is defined as

\[\text{MAE}[\hat{\mathcal{L}}]=\mathbb{E}\left[\big{|}\mathcal{L}(X)-\hat{ \mathcal{L}}(X)\big{|}\right]. \tag{15}\]

We computed this for each classifier as an empirical average over the 100 different likelihood ratio predictors to

\begin{table}
\begin{tabular}{l c c} Loss Name & \(A(f)\) & \(B(f)\) \\ \hline \hline Binary Cross-Entropy & \(\ln(f)\) & \(\ln(1-f)\) \\ Mean Squared Error & \(-(1-f)^{2}\) & \(-f^{2}\) \\ Maximum Likelihood Classifier & \(\ln(f)\) & \(1-f\) \\ Square Root & \(-\frac{1}{\sqrt{f}}\) & \(-\sqrt{f}\) \\ \end{tabular}
\end{table}
Table 1: The rescaling functions \(A\) and \(B\) used to assemble the four different loss functionals considered.
get a numerical measure of how well each predictor approximated the likelihood ratio.

Next, we examined how varying the amount of data upon which the classifiers were trained affected their performance. In particular, for each loss, we trained 100 classifiers for each \(N\in\{10^{2},10^{3},10^{4},10^{5},10^{6},10^{7}\}\). For each value of \(N\), \(0.75N\) observations were used for training and \(0.25N\) observations were used for validation. The value of \(N=10^{6}\) corresponds to our default sample size. As before, 100,000 samples were used to estimate the MAE for each value of \(N\).

#### iii.1.3 Results

Figure 1 displays the likelihood ratio fits averaged over 100 models for each of the four classifiers, compared against the true likelihood ratio. The largest deviations here are in regions far outside the bulk of the training data, where the models will largely be extrapolating. We are primarily concerned with evaluating the likelihood ratio approximation where the data has good coverage: approximately \(x\in[-3,3]\).

In Fig. 2, we show how the expected error for classifiers trained with each choice of loss functional decreases as the sample size increases.

#### iii.1.4 Discussion

The four losses result in similarly performing fits near \(x=0\); however, the MLC and SQR losses rapidly diverge from the true likelihood ratio in regions for which there is little data coverage. By comparison, the BCE and MSE perform much better, staying within 3% of the true likelihood ratio even in regions far outside the bulk of the data (\(|x|>4\)).

The performance of these classifiers varies with the size of the training dataset \(N\). For relatively small training sample sizes (\(N<1000\)), the scale of the mean absolute error is dominated by the inductive bias present in each activation function: BCE and MSE losses (both using \(\sigma(z)\) activation) are nearly identical in size, while MLC and SQR losses (both using \(\mathrm{ReLU}(z)\) activation) are similarly clustered. As \(N\) increases, the MLC and SQR classifier performances approach those of the BCE and MSE classifiers. However, even for values of \(N\) larger than \(10^{5}\), the SQR classifier's MAE remains at least 0.015 above the average performance of the BCE/MSE classifiers.

### Parametrizing \(f\)

The parametrization of the learned function can be adjusted. In the naive implementation, the BCE and MSE neural networks use a logistic activation function, while the MLC and SQR neural networks use a ReLU activation function.

Let \(z(x)\) be the function that the neural network represents. Then \(f=\phi(z)\) is our classifier, where \(\phi\) is some parametrization of the learned function. In the cases described before, we have either \(f=\sigma(z)\) (for BCE and MSE) or \(f=\mathrm{ReLU}(z)\) (for MLC and SQR).

However, for BCE and MSE, any function \(\phi:\mathbb{R}\rightarrow(0,1)\) will suffice. Two readily available such functions

Figure 1: Average likelihood ratio fits for the four different losses. The MAEs are 0.0083, 0.0081, 0.0150, and 0.0254, for the BCE, MSE, MLC, and SQR likelihood ratio models, respectively.

Figure 2: Mean absolute errors computed for the four different losses trained with increasingly larger sample sizes \(N\).
are hyperbolic tangent and arctangent, adjusted to the appropriate range:

\[f(z) =\frac{1}{2}\left(\tanh z+1\right), \tag{16}\] \[f(z) =\frac{1}{\pi}\left(\arctan z+\frac{\pi}{2}\right). \tag{17}\]

In Fig. 3, we show the likelihood ratio fits, averaged over 100 models, for the logistic, hyperbolic tangent, and arctangent parametrizations of the BCE and MSE classifiers. In both cases, the default logistic parametrization performs the best, followed closely by the hyperbolic tangent parametrization, and followed distantly by the arctangent parametrization. This result is not surprising, as the logistic function is known to be well-suited for classification.

For the MLC and SQR losses, we instead require any function \(\phi:\mathbb{R}\rightarrow(0,\infty)\). While the ReLU function is the default, there are other functions with such ranges, including:

\[f(z) =z^{2}, \tag{18}\] \[f(z) =\exp z. \tag{19}\]

Figure 4 displays the results of comparing the performances of the MLC and SQR losses in training classifiers with these parametrizations. The performances of the three parametrizations between the two losses are the same: in this case, the exponential parametrization performs remarkably better than the ReLU parametrization, and square parametrization performs the worst amongst all three.

### Generalized Loss Families

#### iii.3.1 Motivation

The MSE and SQR loss functionals are easily generalizable to a parametric family of loss functionals. While there are several possible parametrizations2 to choose from, we select the following for simplicity: for the MSE loss, we consider a power parameter \(p\in\mathbb{R}\), where \(p=2\) is the default value, and for the SQR loss, we consider a root parameter \(r\in\mathbb{R}\), where \(r=1\) is the default value. This yields the two families of losses presented in Table 2.

Footnote 2: For example, to enforce non-singular behavior at \(r=0\) for SQR, one could consider \(A(f)=(1-f^{-\frac{\pi}{2}})/|r|\) and \(B(f)=(1-f^{\frac{\pi}{2}})/|r|\). Another interesting parametrization is \(A(f)=(f^{q}-1)/q\) and \(B(f)=1-f^{(q+1)}/(q+1)\), which is minimized at \(q=1\).

Since the rescaling functions \(A\) and \(B\) have changed, the likelihood ratio recovered from \(f\) changes as well.

For the \(p\)-MSE losses, for \(p\notin(0,1)\),

\[\mathcal{L}(x) =-\frac{B^{\prime}(f)}{A^{\prime}(f)}=-\frac{-pf^{p-1}\cdot f^{ \prime}}{p(1-f)^{p-1}\cdot f^{\prime}} \tag{20}\] \[=\left(\frac{f}{1-f}\right)^{p-1}. \tag{21}\]

We exclude the case where \(p\in(0,1)\) since the corresponding loss functional is not convex, and as such the likelihood ratio trick no longer works.

And for the \(r\)-SQR losses,

\[\mathcal{L}(x) =-\frac{B^{\prime}(f)}{A^{\prime}(f)}=-\frac{\frac{r}{2}f^{\frac {\pi}{2}-1}\cdot f^{\prime}}{\frac{\pi}{2}f^{-\frac{\pi}{2}-1}\cdot f^{\prime}} \tag{22}\] \[=f^{r}. \tag{23}\]

A whole family of losses arises from both of the two original losses, each loss still maintaining the property that the function that minimizes the corresponding functional can recover the likelihood ratio. In addition to comparing how the four original losses performed against one another, we can compare among the losses in each of these two loss families.

#### iii.3.2 Methods

Since we were working over an uncountably infinite set of loss functionals, we decided to constrain our investigation to just the interval \([-2,2]\). We scanned along the interval \([-2,2]\); for each value \(p\) which we looked at, we trained 20 logistically-parametrized models on the \(p\)-MSE loss functional corresponding to that value of \(p\). Then we averaged the mean absolute errors of the 20 models together.

We did the same for values of \(r\) in the interval \([-2,2]\) as well; in that case, the models were parametrized with the exponential activation function instead.

We expect that near \(p^{*}=1\) and \(r^{*}=0\), where the generalized loss functionals will resemble the MAE loss, the figure-of-merit of MAE will likely be minimized, too. Due to this intrinsic relationship between the choice of loss functional and figure-of-merit, we also considered two additional figures-of-merit for evaluating these scans: the Mean Ratio and the Null Statistic, defined as:

\begin{table}
\begin{tabular}{c c c} Loss Name & \(A(f)\) & \(B(f)\) \\ \hline \hline \(p\)-MSE & \(-(1-f)^{p}\) & \(-f^{p}\) \\ \(r\)-SQR & \(-f^{-\frac{\pi}{2}}\) & \(-f^{\frac{\pi}{2}}\) \\ \end{tabular}
\end{table}
Table 2: The generalization of the MSE and SQR loss functionals to entire families of losses. Values of \(p=2\) and \(r=1\) correspond to the original definitions of the loss functionals.
Figure 4: Parametrizations of \(f\) for the MLC and SQR losses. (a) The average likelihood ratio fits of the ReLU, square, and exponential parametrizations for the MLC loss, with mean absolute errors 0.0148, 0.0684, and 0.0083, respectively. (b) The average likelihood ratio fits of the ReLU, square, and exponential parametrizations for the SQR loss, with mean absolute errors 0.0367, 0.6756, and 0.0075, respectively.

Figure 3: Parametrizations of \(f\) for the BCE and MSE losses. (a) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the BCE loss, with mean absolute errors 0.0080, 0.01240, and 0.0092, respectively. (b) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the MSE loss, with mean absolute errors 0.0084, 0.0127, and 0.0094, respectively.
\[\text{Mean Ratio}[\hat{\mathcal{L}}]=E[\hat{\mathcal{L}}(X)/\mathcal{L}(X)] \tag{24}\]

\[\text{Null Statistic}[\hat{\mathcal{L}}]=|E_{0}[\mathcal{L}(X)]-E_{0}[\hat{ \mathcal{L}}(X)]| \tag{25}\]

We found that the overall trends reported here using MAE were similar across these alternative figures-of-merit, though the trends were less dramatic when we used the Mean Ratio figure-of-merit.

#### iii.3.3 Results

In Figure 5, we show the performance of the classifiers trained by these losses when modifying their power and root parameters. The values \(p^{*}\) and \(r^{*}\) minimizing the MAE were \(p^{*}=1.08,1.24\) (with \(p^{*}=1.24\) having a similar performance to that of \(p^{*}=1.08\) while being more numerically stable) and \(r^{*}=0.018\).

#### iii.3.4 Discussion

In 5(a), we observe vertical features for \(p\in(0,1)\). This is to be expected, as the likelihood ratio trick does not apply in the range where the corresponding loss functional is non-convex. Similarly, the vertical feature in 5(b) is due to the fact that for \(r=0\), our loss functional is constant (\(L[f]=1\)), and thus it is not strictly convex; therefore the likelihood ratio again does not work.

Values of \(p\) slightly less than \(0\) or slightly greater than \(1\) resulted in the smallest mean absolute errors, while values of \(r\) close to \(0\) resulted in the smallest mean absolute errors.

This result was further investigated in Section III.5 in a simple, two-dimensional classifier model.

### Optimized Implementation

Altering the parametrization of the learned function \(f\) or using a more generalized loss functional yielded considerable increases in performance from the initial parametrizations and loss functionals.

In Figs. 6 and 7, we chose the best-performing parametrization for each loss (logistic for BCE and MSE; exponential for MLC and SQR), and, for the MSE and SQR, chose the best-performing loss functional from each loss family (\(p^{*}=1.24\) for MSE and \(r^{*}=0.018\) for SQR), and trained classifiers with each "optimized" parametrization and loss. This was done \(100\) times for each parametrization/loss, and the resulting likelihood ratio models were averaged.

In the naive implementation, the BCE and MSE models performed the best, while the SQR model had an

Figure 5: (a) The mean absolute errors averaged over models trained on the generalized MSE loss family for the logistic parametrization. The mean absolute error is minimized at \(p^{*}=1.08\), but we choose the second-lowest value (\(p^{*}=1.24\)) for stability, i.e. avoiding the steep increase in MAE near \(p=1\). The arrow indicates the typical choice of \(p=2\) for MSE loss. (b) The mean absolute errors averaged over models trained on the generalized SQR loss family for the exponential parametrization. The mean absolute error was smallest at \(r^{*}=0.018\). The arrow indicates the typical choice of \(r=1\) for SQR loss.
average error at least 0.015 larger than the other losses, even for large \(N\). In the optimized implementation with \(N=10^{6}\), all four loss functionals perform approximately equally, as shown in Figure 6. Figure 7 shows that for \(N>10^{5}\), the four optimized loss functionals continue to perform approximately equally well, but the new loss functionals \(p^{*}\)-MSE and \(r^{*}\)-SQR perform significantly better, reaching mean absolute errors about 2 to 4 times smaller than the other losses. The strong influence of the inductive bias of the activation function is also mitigated in the optimized implementation, as the losses are no longer grouped by activation function.

### Simple Classifiers

#### iii.5.1 Motivation

To better understand the behavior of the generalized loss models in Sec. III.3, we examined a much simpler classifier than the multi-layer fully-connected network used to train the models in this paper. This allows us to visualize the dynamics of each model, using numerical integration to compute the loss. The model is

\[f(x)=\phi(ax+b), \tag{26}\]

where \(a\) and \(b\) are the two weights of the model and \(\phi\) is its activation.

In the case of the \(p\)-MSE model with the logistic parametrization,

\[f_{\text{MSE}}(x) =\sigma(ax+b) \tag{27}\] \[\hat{\mathcal{L}}_{\text{MSE}}(x) =\left(\frac{\sigma(ax+b)}{1-\sigma(ax+b)}\right)^{p-1}=\left(e^{ ax+b}\right)^{p-1}. \tag{28}\]

The exponentially-parametrized \(r\)-SQR model is

\[f_{\text{SQR}}(x) =e^{ax+b}\] \[\hat{\mathcal{L}}_{\text{SQR}}(x) =\left(e^{ax+b}\right)^{r}.\]

As a result, only the analysis of one of the two models is necessary, since the resulting likelihood ratio model is the same for \(r=p-1\). In particular, we will analyze the \(r\)-SQR model, keeping in mind that for the model MAEs, the results will be identical for \(p=r+1\).

We continue working with \(X_{0}\sim\text{Normal}(+0.1,1)\) and \(X_{1}\sim\text{Normal}(-0.1,1)\). The exact likelihood ratio is given by

\[\mathcal{L}(x)=e^{0.2x}, \tag{29}\]

so the two-dimensional classifier will yield an exact solution at

\[a^{*}=\frac{0.2}{r},\qquad b^{*}=0. \tag{30}\]

#### iii.5.2 Methods

To better understand how the parameter \(r\) affects the optimization landscape, we first created a grid with fineness 0.005 of \((a,b)\) pairs in the box \([-1,1]^{2}\):

\[B=\frac{1}{200}\mathbb{Z}^{2}\cap[-1,1]^{2}. \tag{31}\]

Figure 6: Average likelihood ratio fits for the different loss categories. The MAEs are 0.0079, 0.0045, 0.0077, 0.0034, 0.0046, and 0.0034, for the BCE, MSE, MLC, SQR, \(p^{*}\)-MSE, and \(r^{*}\)-SQR likelihood ratio models, respectively.

Figure 7: Mean absolute errors computed for the different loss categories trained with increasingly larger samples.
The loss functional for a particular value of \(r\), \(L_{r}\) is given by

\[L_{r}[f]=\int\mathrm{d}x\bigg{(}p(x\mid\theta_{A})f(x)^{-\frac{r}{2}}+p(x\mid \theta_{B})f(x)^{\frac{r}{2}}\bigg{)} \tag{32}\]

Then, we visualized the loss landscape as the contour plot of \(L_{r}\) over the set of classifiers \(F=\{(e^{ax+b})^{r}:(a,b)\in[-1,1]^{2}\}\) for different values of \(r\). The loss functional \(L_{r}\) was computed via numerical integration.

#### iv.2.3 Results

Figure 8 displays in the first row the resulting contour plots for \(r\in\{0.1,0.25,0.5,1\}\). Drawn over each plot, in white, are level sets of the loss at increments of \(0.02\).

#### iv.2.4 Discussion

While the actual values of the losses are not comparable between different values of \(r\), since each value of \(r\) corresponds to a different loss functional, it is clear that the loss functional becomes increasingly steep as \(r\) increases. As expected, as \(r\to\infty\), \(a^{*}\to 0\), and as \(r\to 0\), \(a^{*}\to\infty\). In particular, the loss landscape of \(r=0.1\) is shaped like an extremely shallow pool, indicating that there is a large space of classifiers with close to optimal performance. The minimum value \(a^{*}=2\) is not visible in the box, since small values of \(r\) correspond to large values of \(a^{*}\). On the other hand, the loss landscape of \(r=1\) is much steeper, with a minimum at \(a^{*}=0.2\) around which the landscape quickly increases to high loss values.

However, since loss values between values of \(r\) are incomparable, it is unclear how the loss reflects the actual performance of the likelihood ratio model. In particular, given two classifiers \(f\) and \(g\) with \(L_{r}[f]<L_{s}[g]\), \(r\neq s\), we cannot be sure that \(f\) will yield a better likelihood ratio model than \(g\), since \(L_{r}\) and \(L_{s}\) are different loss functionals.

To this end, we visualized the error landscape as the contour plot of MAE over the same set of classifiers \(F=\{(e^{ax+b})^{r}:(a,b)\in[-1,1]^{2}\}\). Since the MAE is computed only from the expected absolute difference between a predicted likelihood ratio \(\hat{\mathcal{L}}\) and the true likelihood ratio \(\mathcal{L}\), we can compare across different values of \(r\) to see which values of \(r\) result in easily obtainable well-performing classifiers.

The second row of Figure 8 displays these error contour plots; indicated in white are the level sets of the error at increments of \(0.05\). For each of these, we can see that the error is zero at \((a^{*},0)\) and increases radially outwards from the minimum. The shape of the loss landscapes reflect the true nature of the performance of the classifiers; for \(r=0.1\), we still have a shallow pool of many well-performing classifiers, whereas for \(r=1\), there is a small set of well-performing classifiers around which the classifiers begin to perform much worse. That is to say, for small values of \(r\), there are many classifiers that perform well at modeling the likelihood ratio. It may be harder to find the true minimum, but most classifiers have comparable performance. On the other hand, for large values of \(r\), the loss landscape is steep with few

Figure 8: Contour plots of the losses and MAEs of the two-dimensional SQR classifier \(f(x)=\exp{(ax+b)}\) over \([-1,1]^{2}\). The first row plots the value of the loss functional \(L[f]\), obtained through numerical integration, on a grid of \((a,b)\) pairs over \([-1,1]^{2}\) for various values of \(r\), with contours curves at increments of \(0.02\). The second row plots an empirically computed absolute error, \(\mathrm{MAE}[f]\), over the same grid of points, for the same values of \(r\), with contour curves at increments of \(0.05\).
classifiers with decent performance. Slight perturbations around the minimum correspond to large errors.

## IV Multivariate Gaussians

### Parametrizing \(f\)

#### iv.1.1 Motivation

A natural extension from the univariate Gaussians analysis in the previous section would be to multivariate Gaussians, wherein the setting is complicated by the higher dimensions, but we still have knowledge of the true likelihood ratio. To this end, we first established five different case studies of different Gaussian arrangements to examine in our multivariate analysis.

The first case study, labeled "Vertical," corresponds to independent Gaussians with variance 1, and means at a distance of 0.2, as in the univariate case.

In this case, the background distribution is more likely over the right half-plane, whereas the signal distribution is more likely over the left half-plane.

\[X_{0} \sim\text{Normal}\left(\begin{bmatrix}+0.1\\ 0\end{bmatrix},\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\right) \tag{33}\] \[X_{1} \sim\text{Normal}\left(\begin{bmatrix}-0.1\\ 0\end{bmatrix},\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\right) \tag{34}\]

The next case study, "Slant," simply rotates the vertical case study by \(45^{\circ}\). This results in the same likelihood ratio as the vertical case, except rotated by \(45^{\circ}\).

\[X_{0} \sim\text{Normal}\left(\begin{bmatrix}+\frac{0.1}{\sqrt{2}}\\ -\frac{0.1}{\sqrt{2}}\end{bmatrix},\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\right) \tag{35}\] \[X_{1} \sim\text{Normal}\left(\begin{bmatrix}-\frac{0.1}{\sqrt{2}}\\ +\frac{0.1}{\sqrt{2}}\end{bmatrix},\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\right) \tag{36}\]

In "Circle," we consider the case where the background distribution has low variance in comparison to the signal distribution. As a result, values close to the origin are more likely to be from the background, whereas values far from the origin are more likely to be from the signal. This likelihood structure is visualized in Figure 9.

\[X_{0} \sim\text{Normal}\left(\begin{bmatrix}+0.1\\ 0\end{bmatrix},\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\right) \tag{37}\] \[X_{1} \sim\text{Normal}\left(\begin{bmatrix}-0.1\\ 0\end{bmatrix},\begin{bmatrix}2&0\\ 0&2\end{bmatrix}\right) \tag{38}\]

The "Hyperbola" case study looks at the case when both the background and the signal have different variances in each coordinate. This results in a hyperbola-like likelihood structure, as visualized in Figure 9.

#### iv.1.2 Methods

The methodology was similar to that done in Sec. III. For each case study, we implemented all four classifiers with each of the three parametrizations. Each resulting classifier architecture was trained 100 times to minimize the corresponding loss functional. We evaluated each classifier on the box \([-2,2]^{2}\), and we averaged the resulting 100 predictions for the likelihood ratio over that

Figure 9: Two of the five multivariate Gaussian cases we examined, as well as some of the likelihood ratio model fits. The first row corresponds to the circle case, while the second row corresponds to the hyperbola case. The first column plots the likelihood structure of each case; red regions are regions where \(\mathcal{L}(x)\leq 1\), and blue regions are regions where \(\mathcal{L}(x)>1\). The second and third columns display contour plots of the mean absolute error for some models trained with the various losses to learn the likelihood ratios. The plot is suggestively colored to show how the structure of the data corresponds to the structure in the likelihood ratio models.
box. We used the MAE as the performance metric, again as an empirical average over 100,000 samples.

#### iv.1.3 Results

The resulting MAEs are shown in Figure 10. Some contour plots of the mean absolute errors of some of the different parametrizations are presented in Figure 9; the remaining contour plots are provided in Appendix B.

#### iv.1.4 Discussion

In the univariate case, we found that the logistic and exponential parametrizations were uniformly the best parametrizations for the BCE/MSE and MLC/SQR losses, respectively. This trend bears out for the most part in this higher-dimensional case. In almost all of the case studies, the logistic and exponential parametrizations perform the best; in the cases where they don't have the smallest MAE, the difference between their MAE and the best MAE is not large.

Unlike in Section III, once the optimal parametrizations are chosen for each of the four loss functionals, some differences persist in the performance of each loss. Across all five cases, the SQR loss yields the largest errors. For the Vertical and Slant cases, all four optimized loss functionals perform equally well, overlapping within one standard deviation. For the remaining cases (Checker, Circle, and Hyperbola), the optimized MLC loss with exponential parametrization performs significantly better than the other three optimized losses.

It is striking to note that the MLC loss with exponential parametrization emerges as the best-performing loss configuration in some of the more complex datasets considered for these studies. The typical choice for a neural network classifier loss is arguably BCE. For the purposes of the likelihood ratio trick, however, we are interested in reinterpreting the classifier output to approximate the likelihood ratio, so it is possible that optimizing for raw classification performance alone is misguided. The MLC loss has the advantage of explicitly relating the signal and background probability distributions; in particular, the MLC loss can be intuitively understood to maximize the likelihood of \(\hat{\mathcal{L}}(x)\) with respect to \(p(x\mid\theta_{0})\) subject to the constraint that \(\hat{\mathcal{L}}(x)p(x\mid\theta_{1})\) is a probability distribution [10]. Therefore, it may be a more natural choice for this particular application than the default BCE loss.

### Generalized Loss Families

#### iv.2.1 Motivation

By treating the square in the MSE loss and the root in the SQR loss as parameters \(p\) and \(r\), respectively, we were able to generalize those loss functional to entire continuous parametric families of losses. We saw in the univariate case that we can optimize over \(p\) and \(r\), and were even able to see, through examining the landscapes of the different loss functionals in a simple case, what kinds of values of \(p\) and \(r\) will correspond to "better" loss functionals.

We now continue this investigation in the situation of multivariate Gaussians to get a sense of how much the trend we observe continues into more complex situations.

#### iv.2.2 Methods

We used the same methods as in Section III.3; in this case, however, we worked with the five different multivariate Gaussians cases rather than the single univariate Gaussians case.

#### iv.2.3 Results

In Table 3, we list the optimal values \(p^{*}\) and \(r^{*}\) in each of the five cases. An overall comparison of the four loss functionals with optimized parameterizations alongside \(p^{*}-\)MSE and \(r^{*}-\)SQR losses is shown in Figure 11. The plots of the MAE for the various values of \(p\) and \(r\) are presented in Appendix B.

#### iv.2.4 Discussion

The simpler multivariate cases considered (Vertical, Slant) result in very similar values to those found in the univariate Gaussian case: \(p^{*}\) and \(r^{*}\) are close to \(1.24\) and \(0.018\), respectively.

In the more complex multivariate cases (Circle, Hyperbola, and Checker), the optimal values of \(p^{*}\) are also between \(1\) and \(2\), with the expection of Hyperbola, for which \(p^{*}=-0.44\). It's possible that an equally-performing

\begin{table}
\begin{tabular}{l c c} Case & \(p^{*}\) & \(r^{*}\) \\ \hline \hline Vertical & \(1.12\) & \(0.018\) \\ Slant & \(1.16\) & \(0.018\) \\ Circle & \(1.28\) & \(-0.1\) \\ Hyperbola & \(-0.44\) & \(-0.2\) \\ Checker & \(1.6\) & \(-0.1\) \\ \end{tabular}
\end{table}
Table 3: The optimal values for \(p^{*}\) and \(r^{*}\) for the five different multivariate Gaussian cases. Note that in the univariate Gaussian case, the optimal values chosen were \(p^{*}=1.24\) and \(r^{*}=0.018\).
value of \(p\) larger than 2 could also exist, but our studies did not scan far enough to probe the asymptotic behavior in that direction. The optimal values \(r^{*}\) are all negative. However, it is worth noting that the MAE landscapes for the \(r\)-SQR are symmetric, and the corresponding MAEs for \(|r^{*}|\) are small, so the signs of these values is likely due to random chance. For these cases, the optimal values of \(r^{*}\) are less than 1, as in the univariate case, but very small values of \(r^{*}\) (\(|r^{*}|<0.01\)) are too numerically unstable to consistently yield useful outputs.

Overall, as shown in Figure 11, if one chooses only from the four loss functionals as defined in Table 1, but with optimized parametrizations, all four show equally good performance for the simpler cases (Vertical, Slant), but the MLC loss is significantly better than the other three choices in the more complex cases (Circle, Hyperbola, Checker). However, if one chooses \(p^{*}\) and \(r^{*}\) by scanning along these generalized loss families, the improvements are immense: for all cases except Hyperbola, the optimized \(r^{*}-\)SQR MAE is between 30% and 50% smaller than the optimized MLC MAE.

## V Physics data

### Parametrizing \(f\)

#### v.1.1 Motivation

In our final case study, we extended our comparison of classifier parametrizations and loss functionals to simulated high-energy particle physics datasets [39]. While there are a number of observables present in the datasets, in our analysis we considered only the leading jet transverse momentum (\(p_{T}\)), rapidity (\(y\)); azimuthal angle (\(\phi\)), and invariant mass (\(m\)).

The datasets consist of particle-level and detector-level simulated QCD jets originating from \(Z\) + jets events. \(Z\) + jets events from proton-proton collisions generated at \(\sqrt{s}=14\) TeV were simulated using Herwig 7.1.5 [40, 41, 42] with the default tune and Pythia 8.243 [43, 44, 45] tune 21 [46] (ATLAS A14 central tune with NNPDF2.3LO). We call the Pythia simulation "Monte Carlo" and Herwig "data". For the generated events, the \(p_{T}\) of the \(Z\) boson is required to be larger than 150 GeV. Events then

Figure 10: Mean absolute errors are compared for the four different losses considered (binary cross-entropy, mean squared error, maximum likelihood classifier, and square root), each with 3 different activation functions. For each loss, five different multivariate normal cases are studied: “Vertical”, “Slant”, “Circle”, “Checker”, and “Hyperbola”. For each case study, the best performing parametrization for each loss is shown in either red or blue. Errors represent the standard deviation across 100 independent model trainings.
Figure 11: Mean absolute errors are compared for the four different losses considered (binary cross-entropy, mean squared error, maximum likelihood classifier, and square root), each with their respective optimal parametrizations. For each loss, five different multivariate normal cases are studied: “Vertical”, “Slant”, “Circle”, “Checker”, and “Hyperbola”. Errors represent the standard deviation across 100 independent model trainings.

Figure 12: Histograms of the four jet features, transverse momentum \(p_{T}\) [GeV], rapidity \(y\), azimuthal angle \(\phi\), and mass \(m\) [GeV], for the Monte Carlo and data \(Z\) + jet events. Here, we treat the Monte Carlo as the signal and the data as the background.
are passed through the Delphes 3.4.2 fast detector simulation [47] of the CMS detector. The datasets consist of the highest-momentum jet from \(Z\) boson events with \(p_{T}\geq 200\) GeV. This process ultimately yields about 1.6 million jets for each simulation. Figures 12 and 13 display histograms of each of the four observables for both the "Monte Carlo" and the "data".

In this more complex setting, we no longer have access to the true likelihood ratio, as we do not know the underlying distributions generating these datasets. To allow for a more complete comparison of the different parametrizations' ability to model the "true" likelihood ratio, we therefore fit Normalizing Flows [48] to each sample. These flows estimate the generating distribution of the samples, and thus allow us to compute "true" likelihood ratios for these datasets.

#### iv.2.2 Methods

We first trained a FFJORD Normalizing Flow [48] for each of the "Monte Carlo" and "data" simulated samples. The models were tuned by comparing the performance of a classifier on distinguishing between data generated from

Figure 13: A corner plot of the four jet features. Blue contours correspond to the particle-level data and red contours correspond to the detector-level data.
the flow and the true data. We then used the flows as proxies for the underlying distributions of these datasets, creating new proxy datasets by sampling from the flows.

The methodology following this point was again similar to what has been established before in Sec. III. We implemented all four classifiers with each of the three parametrizations on the dataset, training 100 independent copies of each classifier architecture to minimize the corresponding loss functional. We used the MAE as the performance metric, computed as in Equation 15. In particular, it was computed with the true likelihood ratio \(\mathcal{L}(X)\) from the flows and the model likelihood ratio \(\hat{\mathcal{L}}(X)\) averaged over the 100 copies of each classifier. As before, the MAE was computed as an empirical average over 100,000 samples.

#### iv.1.3 Results

The distributions of the "data" and "Monte Carlo" learned by the flows are plotted along side the empirical distributions in Figure 12. To quantify the quality of the flows' learned distributions, we trained classifiers to try to distinguish between proxy datasets sampled from the flows and the original datasets; for the "Monte Carlo", the AUC was 0.54, and for the "data", the AUC was 0.56. These AUCs close to 0.5 indicate that the classifier has difficulty distinguishing between these two distributions, and therefore that the flows have performed reasonably well at reflecting the target distributions.

To visualize the performance of these classifiers, we performed a scan of the likelihood ratio along \(\phi\). We fixed the three other observables at values near their medians (\(p_{T}=221.8\), \(y=0.0\), \(m=16.0\)) and compared the flow likelihood ratio to the likelihood ratio modeled from the variously parametrized classifiers trained on the different loss functionals.

The scans for the logistic, arctangent, and hyperbolic tangent parametrizations of the likelihood ratio, trained with the BCE and MSE losses, are displayed in Figure 14. The analogous scans for the ReLU, square, and exponential parametrizations with the MLC and SQR losses are displayed in Figure 15.

The resulting MAEs are shown in Figure 16.

#### iv.1.4 Discussion

The observed trend in the Gaussian studies that the log-gistical and exponential parametrizations were the best for the BCE/MSE and MLC/SQR losses, respectively, also holds in the physics case. Of the four optimized loss functionals, the MLC loss with exponential parametrization performs distinctly better than the other three loss configurations, as shown in Figure 16.

### Generalized Loss Families

#### iv.2.1 Motivation

In the previous studies with univariate and multivariate Gaussians, we found that the performances of likelihood ratio models trained with losses from the generalized families of \(p\)-MSE and \(r\)-SQR losses followed a similar structure across various cases. In order to examine the robustness of this observed structure, we repeated the same study with the high energy particle physics dataset.

#### iv.2.2 Methods

The methodology for this study was similar to that of the previous studies. We scanned over values of \(p\) and \(r\) in the interval \([-2,2]\). For each increment of \(p\), we trained 20 models with the \(p\)-MSE loss functional defined by that value of \(p\) and averaged together their mean absolute errors. Likewise, for each increment of \(r\), we trained 20 models with the \(r\)-MSE loss functional defined by that value of \(r\) and averaged together their mean absolute errors. We parametrized the \(p\)-MSE classifiers with logistic activation functions and the \(r\)-SQR classifiers with exponential activation functions. All models were trained on the same set of one million samples from the flows fit to the distributions of the physics data.

#### iv.2.3 Results

The plots of the MAEs of the likelihood ratio models for the loss functionals are provided in Figure 17. As before, we observe vertical features in the plots when the loss functional is no longer strictly convex (\(p\in(0,1)\) and \(r=0\)). The MAE was minimized at \(p^{*}=-0.28\) and \(r^{*}=-1.3\). A comparison of the \(p^{*}-\)MSE and \(r^{*}-\)SQR losses with these chosen values alongsize the other four losses with optimized parametrizations is shown in Figure 16.

#### iv.2.4 Discussion

The shape of the \(p\) and \(r\) scans looks approximately similar to those observed in the previous case studies (e.g. Figure 5); however, since the MAE landscape is flat away from the non-convex regions (\(p\in(0,1)\) for \(p\)-MSE and \(r=0\) for \(r\)-SQR), the best choices \(p^{*}=-0.28\) and \(r^{*}=-1.3\) perform about the same as the unoptimized choices of \(p=2\) and \(r=1\). In this particular case, the evidence does not suggest that changing \(p^{*}\) or \(r^{*}\) from their default values of 2 and 1, respectively, would yield a significant benefit in reducing mean absolute error. It is possible that better values exist beyond the ranges \(r,p\in[-2,2]\) considered here. Overall, as shown
in Figure 16, the best-performing loss is MLC with exponential parameterization.

## VI Conclusions

The likelihood ratio \(\mathcal{L}(x)=\frac{p(x|\theta_{0})}{p(x|\theta_{1})}\) is a statistical quantity essential for characterizing whether an experimental dataset \(x\) better supports one of two hypotheses defined by sets of parameters \(\theta_{0}\) and \(\theta_{1}\). It is used beyond hypothesis testing, too, for applications such as reweight

Figure 14: The performance of different parametrizations of \(f\) for the BCE and MSE losses for the azimuthal angle \(\phi\). (a) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the BCE loss. (b) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the MSE loss.

Figure 15: The performance of different parametrizations of \(f\) for the MLC and SQR losses for the azimuthal angle \(\phi\). (a) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the BCE loss. (b) The average likelihood ratio fits of the logistic, hyperbolic tangent, and arctangent parametrizations for the MSE loss.
Figure 16: The MAEs are compared for the Pythia/Herwig + Delphes particle physics jet datasets [39] for the four different losses considered. Errors represent the standard deviation across 100 independent model trainings. In **(a)**, each loss is shown with 3 different parametrizations. In **(b)**, the best-performing parametrization is chosen for each loss, and these optimized losses are then directly compared.

Figure 17: (a) The mean absolute errors averaged over logistically-parametrized models trained with the generalized MSE loss family. The mean absolute error is minimized at \(p^{*}=-0.28\); however, the MAE landscape here is rather flat. The arrow indicates the typical choice of \(p=2\) for the standard MSE loss. (b) The mean absolute errors averaged over exponentially-parametrized models trained with the generalized SQR loss family. The mean absolute error was smallest at \(r^{*}=-1.3\). The arrow indicates the typical choice of \(r=1\) for the standard SQR loss.
ing high-dimensional distributions for background estimation and more. In contexts where calculating the likelihood ratio is impossible or very tedious, researchers can use the "likehood ratio trick", leveraging a neural network classifier to approximate the likelihood ratio.

Often, the likelihood ratio trick is implemented by minimizing a typical choice of loss functional for a classifier: the binary cross-entropy loss. However, many loss functionals satisfy the likelihood ratio trick setup.

In this paper, we presented detailed studies comparing four choices of loss functionals: binary cross-entropy (BCE), mean squared error (MSE), maximum likelihood classifier (MLC), and square root (SQR). For each of these four loss functionals, we also explored a suite of choices of final activation functions for parametrizing the neural network output. For the MSE and SQR losses, we performed a scan along the exponential parameter (replacing \(2\to p\) for MSE and replacing \(\frac{1}{2}\rightarrow\frac{r}{2}\) for SQR) to understand the behavior of these generalized families of loss functionals.

As a result of these studies, we present the following recommendations for optimized implementations of each of these loss functionals in the likelihood ratio trick:

\[\begin{split}\text{Loss}&\text{Activation}\\ \hline\hline\text{Binary Cross-Entropy (BCE)}&\sigma(z)\\ \text{Mean Squared Error (MSE)}&\sigma(z)\\ \text{Maximum Likelihood Classifier (MLC)}&\exp(z)\\ \text{Square Root (SQR)}&\exp(z)\\ \end{split}\]

For MLC and SQR losses, we find that choosing small, nonzero values of \(r\) (and, correspondingly, \(p=r+1\)) tend to result in smaller mean absolute errors than the default choices (\(r=1\) and \(p=2\)) for these loss functionals. As we illustrate by mapping the loss landscape of a simple neural network, this is because smaller values of \(r\) can yield shallower loss landscapes where many values are nearly optimal, while larger values of \(r\) have steeper landscapes for models to traverse, with a much smaller proportion of the phase space corresponding to optimum values of the loss.

The loss landscape will vary with each new application, so we recommend that future researchers perform a scan along \(p\) or \(r\) to find an optimum value as part of hyperparameter optimization. If a scan over \(p\) or \(r\) is not feasible, we recommend comparing the default selections (i.e. \(p=2\) and \(r=1\)) with our alternative recommendations derived from the average optimum values across our various trials \(p^{*}=1.15\) and \(r^{*}=0.1\), or:

\[L_{\text{MSE}^{*}}[f]=-\int\mathrm{d}x\bigg{(} (1-f)^{1.15}p(x\mid\theta_{0})\] \[+f^{1.15}p(x\mid\theta_{1})\bigg{)}\] \[L_{\text{SQR}^{*}}[f]=-\int\mathrm{d}x\bigg{(} f^{-0.05}p(x\mid\theta_{0})+f^{0.05}p(x\mid\theta_{1})\bigg{)}.\]

Across the majority of the various datasets we considered, these choices tend to have significantly smaller mean absolute errors than the default selections while maintaining good numerical stability across multiple trainings. An interesting future investigation would be to consider how to dynamically optimize \(p\) and \(r\) as learned parameters during training.

When tested on univariate Gaussians and simple multivariate Gaussians (Vertical and Slant cases), all four loss implementations with optimized parametrizations perform similarly when approximating the desired likelihood ratio. For larger datasets (\(N>10^{5}\)), choosing different exponents in the definitions of MSE and SQR loss functionals results in an additional \(\geq 50\%\) reduction in errors for these cases.

On more complex datasets, including multidimensional Gaussians (Checker, Hyperbola, Circle) as well as simulated high-energy physics data, the Maximum Likelihood Classifier (MLC) loss with exponential parametrization performs the best out of the four default losses considered. Choosing different exponents in the definitions of MSE and SQR loss functionals additionally results in between \(30\%\) and \(50\%\) smaller errors for the Checker and Circle cases. For the Hyperbola and simulated high-energy physics case, choosing alternate \(p^{*}\) and \(r^{*}\) values in the range \([-2,2]\) does not yield a significant performance improvement, though it is possible that better values could exist outside of this range.

While these configurations performed well in our chosen case studies, these results should not be read as a guarantee that these choices will result in optimal performance for any dataset. We therefore recommend that other researchers compare the results of several of the optimized losses described in this work to yield the most effective setup for a given dataset.

There remain several open questions in this line of inquiry. For instance, can an analytical analysis of these loss functionals explain some of the performance differences observed? How much can we further characterize the uncountably many possible loss functionals that satisfy this setup? How else can we generalize certain loss functionals? Pursuing these answers could help us achieve even better scientific measurements enabled by machine learning in the near future.

## Acknowledgements

We are grateful to Jesse Thaler for very helpful feedback about figures-of-merit and ways to generalize our loss functions. We thank Dag Gillberg for the suggestion to compare NNs with BDTs and Vinicius Mikuni for the idea of using normalizing flows to model the LR of the physics datasets. M.P. thanks Shirley Ho and the Flatiron Institute for their hospitality while preparing this paper. S.R., M.P., and B.N. are supported by the Department of Energy, Office of Science under contract number DE-AC02-05CH11231.