Evidence Networks: simple losses for fast, amortized, neural Bayesian model comparison

Bayesian model comparison aims to evaluate the relative posterior probability of different underlying models given some observed data x_O . The ratio of probabilities (the posterior odds) for two models, M₁ and M₀, can be related to the model evidence via:

$$\begin{align} \frac{p\left(M_1 | x_O\right)}{p\left(M_0 | x_O\right)} = \frac{p\left(x_O | M_1 \right)}{p\left( x_O | M_0 \right)} \ \frac{p\left( M_1 \right)}{p\left( M_0 \right)}, \end{align} \tag{ 1 }$$

where $p(x_O | M_1 )$ is the evidence for model M₁ and $p(M_1)$ is the prior probability of the model. The evidence ratio for the two models is known as the Bayes factor, K (see Jaynes 2003 for details).

Under the assumption of equal prior probability for both models, the Bayes factor (i.e. the ratio of model evidences) becomes equal to the ratio of posterior probabilities for models:

$$\begin{align} K = \frac{p\left(x_O | M_1 \right)}{p\left( x_O | M_0 \right)} = \frac{p\left(M_1 | x_O\right)}{p\left(M_0 | x_O\right)} \ , {\mathrm{if} \ } p\left( M_1 \right) = p\left( M_0 \right). \end{align} \tag{ 2 }$$

The evaluation of the Bayes factor, given some data and taking into account all sources of uncertainty, is typically the ultimate goal for model comparison. Bayes factors have been used or discussed in the fields of epidemiology, engineering, particle physics, cosmology, law, neuroscience, and many others (Wakefield 2009, Yuen 2010, Jasa and Xian 2012, Feroz 2013, Fenton et al 2016, Handley et al 2019, Keysers et al 2020, Massimi 2020). Though the simplest interpretation is that of an odds ratio under equal model priors, interpretations like the Jeffreys' scale (Jeffrey 1998) are also popular for Bayesian model comparison. In general these interpretations are concerned with the log Bayes factor, thought of as the relative magnitude of the odds.

In section 1.3, 'Motivation' we introduce the challenges faced in Bayesian model comparison that are solved by Evidence Networks: intractability of computationally challenging integration, unknown likelihood/prior, and lack of model parameterization. In section 1.4, 'Evidence Networks & comparison with existing methods', we consider existing methods in the context of these three motivations.

In section 2, 'Optimization objective', we show how to construct bespoke losses to estimate convenient functions of the model evidence, with focus on the log Bayes factor. Section 2.1 describes the first case of the 'Direct estimate of model posterior probability'. Section 2.2 introduces our primary case of 'Direct estimate of Bayes factors', which includes the description of the novel 'l-POP' transform and associated 'l-POP-Exponential' loss. We include a 'Summary of "symmetric losses"' in section 2.3.

In section 3, we discuss 'Pitfalls & validation techniques' with Evidence Networks. Section 3.1 describes the practical 'Challenge' in validating a trained Evidence Network in the standard case in which the true Bayes factor of the validation data is not known. Section 3.2 presents a 'Blind coverage testing' solution.

Section 4 presents 'Demonstration: Analytic time-series'. This demonstration uses a 'Generative model' (section 4.2) defined such that the Bayes factors can be calculated analytically using a closed-form expression. Section 4.3, 'Evidence Network results', describes the results of the Evidence Networks, demonstrates the blind coverage test, and explores the choice of loss functions. Alternatives are considered in section 4.4, 'Comparison to state-of-the-art: nested sampling with known likelihood', and section 4.5, 'Alternative likelihood-free method: density estimation'.

Section 5 presents 'Demonstration: Dark Energy Survey data', in which Evidence Networks are used to compare two alternative models of galaxy intrinsic alignments using gravitational lensing data.

Section 6, 'Pedagogical Example: Rastrigin Posterior', addresses a conceptual point: there are cases where estimation of the Bayes factor (e.g. with Evidence Networks) for model comparison can be done even in cases where solving the parameter inference problem is intractable, e.g. due to the complexity of the posterior probability distribution.

We conclude in section 7, with section 7.1 presenting a discussion of 'Simple extensions for further applications': absolute evidence (rather than Bayes factor ratios), the connection to frequentist hypothesis testing, and Evidence Networks for posterior predictive tests.

The challenge for Bayesian model comparison arises in the calculation of the model evidence. Assuming the statistical model is known, the model evidence can be expressed as the marginal likelihood:

$$\begin{align} p\left(x_O | M_1\right) = \int p\left(x_O | \theta, M_1 \right) \ p\left(\theta | M_1 \right) \ \mathrm{d} \theta . \end{align} \tag{ 3 }$$

This integral over the space of permissible model parameters θ is often intractable. There could be three distinct reasons for this:

• (i)  
  The integral is computationally challenging. This often renders the calculation effectively impossible for a high-dimensional parameter space.
• (ii)  
  The probability densities $p(x_O | \theta)$ and $p( \theta)$ are unknown.
• (iii)  
  There is no known underlying parameterization that generates the data. For example, real-world data can be drawn from different model classes, and there may be no evident underlying parameterization.

We expand upon these three reasons in section 1.4 below. However, for practitioners used to the marginal likelihood presentation of Bayesian evidence (equation (3)), the key points to consider first are:

• model evidence does not require an integration over model parameters.
• model evidence does not require knowledge of the unnormalized posterior distribution.

The marginal likelihood is a technique to marginalize the joint distribution $p(x_O, \theta | M_1)$ to recover $p(x_O | M_1)$, however a knowledge of $p(x_O, \theta | M_1)$ is not actually necessary to evaluate the model evidence $p(x_O | M_1)$. Evidence Networks provide an alternative approach that avoids equation (3) entirely.

Evidence Networks are fast amortized neural estimators of Bayes factors for model comparison. In section 2 we describe classes of loss functions that lead to the correct optimization objective for achieving this goal. In section 3 we describe the potential pitfalls and validation techniques that are applicable to real-world applications.

We will now address the three motivations highlighted above, comparing to existing methods and linking to our demonstrations that are presented within this paper.

1.4.1. Motivation (i) - intractable integral

1.4.2. Motivation (ii) - unknown likelihood/prior

1.4.3. Motivation (iii) - no parameterization

Though we favor direct estimation of functions of the Bayes factor, we will first review simple optimization objectives that result in neural network estimates of the posterior model probability $p(M_i|x_O)$.

The first of these cases, the Squared Loss, can be interpreted as simple special case of parameter inference (e.g. in the context of Moment Networks: Jeffrey and Wande 2020). The alternative case, the Cross-Entropy Loss, matches the $\mathcal{M}$-closed limit (i.e. all possible models are accounted for such that $\sum_i p(M_i) = 1$) as described in Radev et al (2020). In section 4 we will show that, in comparison to direct Bayes factor estimation, this approach of direct model posterior estimation is not ideal for model comparison.

2.1.1. Squared & Polynomial Loss

First let us consider the typical Squared Loss

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = \left(f\left(x\right) - m \right)^2 . \end{align} \tag{ 7 }$$

The functional to be optimized from equation (6) is then given by:

$$\begin{align} I\left[\,f\,\right] = \int \ \left(f\left(x\right) - 1 \right)^2 \ p\left(x, M_1\right) + f\left(x\right)^2 \ p\left(x, M_0\right) \ \mathrm{d}x . \end{align} \tag{ 8 }$$

If we minimize this objective using standard methods from variational calculus, i.e. $\delta I = 0$, we recover the solution for the global optimal function:

$$\begin{align} f^{\;*}\left(x\right) = p\left(m = 1 | \ x\right) = p\left(M_1 | \ x\right) \ \forall \ x , \end{align} \tag{ 9 }$$

under the assumption that there are two models (so that $p(M_1 | x) = 1 - p(M_0 | x)$). Therefore, the optimal function, evaluated for the observed data x_O , gives the model posterior for model M₁:

$$\begin{align} f^{\;*}\left(x_O\right) = p\left(M_1 | \ x_O\right). \end{align} \tag{ 10 }$$

We can further generalize to a Polynomial Loss with a power $\alpha \in \mathbb{R}$, α ≠ 1,

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = m\left(1-f\left(x\right)\right)^\alpha + \left(1-m\right) f\left(x\right)^\alpha , \end{align} \tag{ 11 }$$

we recover the same result as for the Squared Loss, but with the model posterior probability raised to a power. This results in:

$$\begin{align} f^{\;*}\left(x_O\right) = p\left(M_1 | \ x_O\right)^\frac{1}{\alpha-1} , \end{align} \tag{ 12 }$$

where α = 2 corresponds to the Squared Loss result, as expected. The choice of α leads to different loss landscapes during training, meaning α is a possible hyper-parameter that can be tuned depending on the use-case.

2.1.2. Cross-entropy loss

If we take the loss function

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = m \log \left(f\left(x\right) \right) + \left(1 - m\right) \log \left( 1 - f\left(x\right) \right) \ \ , \end{align} \tag{ 13 }$$

following the same procedure as before gives the result:

$$\begin{align} f^{\;*}\left(x_O\right) = p\left(M_1 | \ x_O\right) \ \ . \end{align} \tag{ 14 }$$

2.2.1. Exponential Loss

Taking the loss function

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) & = e^{\left(\frac{1}{2} - m \right)f\left(x\right)} \end{align} \tag{ 15 }$$

gives the global optimization objective:

$$\begin{align} I\left[\,f\,\right] & = \int \ e^{-\frac{1}{2}f\left(x\right)} \ p\left(x , M_1\right) \ + e^{\frac{1}{2}f\left(x\right)} \ p\left(x , M_0\right) \ \mathrm{d}x. \end{align} \tag{ 16 }$$

Following the same procedure as previously, i.e. minimizing this objective using variational calculus, recovers the optimal function:

$$\begin{align} f^{\;*}\left(x\right) = \log \left( \frac{p\left(M_1 | \ x\right)}{p\left(M_0 | \ x\right)} \right) \ \forall \ x. \end{align} \tag{ 17 }$$

Therefore, the optimal function evaluated for our observed data is:

$$\begin{align} f^{\;*}\left(x_O\right) = \log \left( K \ \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right) = \log K + \log \left( \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right). \end{align} \tag{ 18 }$$

If the model priors are equal (having the same number of examples of each model label in the training data for the neural network), then the second term vanishes.

2.2.2. Logistic loss

An alternative form gives the same final result as the Exponential Loss, with

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = \log \left( 1 + e^{\left(1-2m\right)f\left(x\right)} \right). \end{align} \tag{ 19 }$$

Following the same optimization procedure as before once again results in an estimate of the log Bayes factor with an additional model prior ratio term:

$$\begin{align} f^{\;*}\left(x_O\right) = \log K + \log \left( \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right). \end{align} \tag{ 20 }$$

2.2.3.  α-Exponential and log-α-Exponent

Analogously to the generalization from Squared Loss to Polynomial Loss with the α hyper-parameter, we can introduce a power of α to the losses that directly estimate the Bayes factor or log-Bayes factor.

The 'α-Exponential' Loss is given by:

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = \left( 1 + e^{\left(1 - 2m \right)f\left(x\right)} \right)^\alpha. \end{align} \tag{ 21 }$$

Following the same optimization procedure as before, we recover results the estimate of the rescaled log Bayes factor with an additional model prior ratio term:

$$\begin{align} \begin{split} f^{\;*}\left(x_O\right) & = \log K^\frac{1}{1+\alpha} + \log \left( \frac{p\left(M_1 \right)}{p\left(M_0 \right)}\right)^{\frac{1}{1+\alpha}} \ \ , \\ & = \frac{1}{1+\alpha} \log K + \frac{1}{1+\alpha} \log \left( \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right). \end{split} \end{align} \tag{ 22 }$$

Similarly, we can construct the 'α-log-Exponential' Loss:

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = f\left(x\right) ^{\left(\frac{1}{2}-m\right)\alpha} . \end{align} \tag{ 23 }$$

For this case, we analogously recover the Bayes factor to a power, multiplied by a model prior ratio factor:

$$\begin{align} f^{\;*}\left(x_O\right) = K^\frac{1}{1+\alpha} \times \left( \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right)^{\frac{1}{1+\alpha}} . \end{align} \tag{ 24 }$$

2.2.4. l-POP-Exponential Loss

Finally, we present the 'l-POP-Exponential Loss', which is our recommended choice of default loss for the Evidence Network. This loss applies a bespoke transformation to the network output followed by the usual Exponential Loss. This is completely equivalent to defining a final layer activation function; for this discussion we include the transformation within the loss to aid interpretability.

We first define the parity-odd power transformation, as:

$$\begin{align} \tilde{\mathcal{J}}_\alpha \left(x\right) : = x \ |x|^{\alpha-1} \ \ \forall \ x \in \mathbb{R} \ {\mathrm{where}} \ \alpha \in \mathbb{R} \ {\mathrm{and}} \ \alpha \unicode{x2A7E} 1 . \end{align} \tag{ 25 }$$

Using the sign function (sgn), another way of representing the parity-odd power transformation is:

$$\begin{align} \tilde{\mathcal{J}}_\alpha \left(x\right) = {\mathrm{sgn}}\left(x\right) \ |x|^{\alpha} . \end{align} \tag{ 26 }$$

This function is differentiable everywhere. Note, however, that the intermediate $\textrm{sgn}$ or absolute value functions are not differentiable at x = 0, so care should be taken in code implementation.

Though the parity-odd power transform $\tilde{\mathcal{J}}$ works well, its value and gradient remains close to zero, an effect which increases with increasing α > 1. We therefore use the leaky parity-odd power, or l-POP, transform, which we define as:

$$\begin{align} {\mathcal{J}}_\alpha \left(x\right) : = x + x \ |x|^{\alpha-1}. \end{align} \tag{ 27 }$$

We could choose to generalize this further, using $\beta x + x \ |x|^{\alpha-1}$, though we effectively choose β = 1 for all that follows. The resulting parity-odd power Exponential Loss, or 'l-POP Exponential Loss' is given by:

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = e^{\left(\frac{1}{2} - m \right)\mathcal{J}_\alpha \left(f\left(x\right)\right)} . \end{align} \tag{ 28 }$$

As discussed above, this is just equivalent to using $\mathcal{J}_\alpha$ as a novel choice of activation function in the network output layer.

Following the same optimization procedure as before, we recover a rescaled estimate of the log Bayes factor plus an additional model prior ratio term:

$$\begin{align} \begin{split} f^{\;*}\left(x_O\right) & = \mathcal{J}^{-1}_\alpha \left( \log K + \log \frac{p\left(M_1 \right)}{p\left(M_0 \right)} \right) \\ & = {\mathrm{sgn}} \left( \log K + \log p\left(M_1 \right) - \log p\left(M_0\right) \right) |\log K + \log p\left(M_1 \right) - \log p\left(M_0 \right) |^\frac{1}{\alpha} \ \ . \end{split} \end{align} \tag{ 29 }$$

For equal model priors, $p(M_1 ) = p(M_0 )$, we expect that

$$\begin{align} \log K = \mathcal{J}_\alpha \left(f\left(^*x_O\right)\right), \end{align} \tag{ 30 }$$

for the optimal trained network $f^{\;*}$.

As with many of the losses that we have previously described, the l-POP-Exponential Loss results in a (transformed) estimate of the log Bayes factor, which is ideal for Bayesian model comparison (e.g. the Jeffreys Scale).

We recommend a default value of α = 2. This provides the reweighting of errors for high Bayes factor values, ensuring accuracy across many orders of magnitude, while being sufficiently simple that the network trains consistently well.

The l-POP transformation $\mathcal{J}_\alpha$ to act on f(x) is not the only available choice that could be used to construct useful losses. Any strictly monotonic function $\mathcal{J}: \mathbb{R} \rightarrow \mathbb{R}$ would fulfil the condition that $\mathcal{J}(f(^*x_O)) = \log K$. Nevertheless, the choice is not arbitrary. For example, we have considered the sinh function in the place of $\mathcal{J}_\alpha$, but find the exponential scaling for high K leads to unreliable performance. Since we have empirically found our l-POP transform $\mathcal{J}_\alpha$ to give accurate results for several test cases presented in this paper we recommend it as the default choice.

In table 1, we summarize both the losses presented in the previous section and the corresponding optimal network outputs (the latter shown in terms of model posteriors or Bayes factors).

Table 1. A set of simple loss functions that can be used to estimate useful functions of the model posterior probability (in particular the Bayes factor K for Bayesian model comparison).

	$\mathcal{V}(f(x), m)$	$f^{\;*}(x_0)$
Loss name	using model label $m \in \{0,1\}$	if $p(M_1) = p(M_0)$
Polynomial	$m(1-f(x))^\alpha + (1-m) f(x)^\alpha$	$p(M_1 | \ x_O) ^{\frac{1}{\alpha - 1}}$
Cross Entropy	$- m \log \big(f(x) \big) - (1 - m) \log \big( 1 - f(x) \big)$	$p(M_1 | \ x_O)$
Exponential	$e^{(\frac{1}{2} - m )f(x)}$	$\log K$
Logistic	$\log \big( 1 + e^{(1-2m)f(x)} \big)$	$\log K$
α-Exponential	$\big( 1 + e^{(1 - 2m )f(x)} \big)^{\alpha-1}$	$\frac{1}{\alpha} \log K$
α-log-Exponent	$f(x) ^{(\frac{1}{2}-m)\alpha}$	$K^{\frac{1}{\alpha}}$
l-POP-Exponential	$e^{(\frac{1}{2} - m )\mathcal{J}_\alpha (f(x))}$	$\mathcal{J}_\alpha^{-1}(\log K)$

This list is not at all comprehensive; we have not attempted to study the family of all loss functions that can be used to construct an Evidence Network. Nor are all of these loss functions novel to this paper. The Squared Loss (a special case of the Polynomial), Cross-Entropy, Exponential, and Logistic have long been known in terms of Bayes risk for classification (Masnadi-Shiraz and Vasconcelos 2008). The '$\alpha-$'type losses are new variations on existing losses and the POP or l-POP Exponential losses are novel to this paper. The main aim, however, is to use these losses in such a way that we can perform quantitative Bayesian model comparison via an Evidence Network.

Furthermore, we have included only 'symmetric losses', i.e. those that are symmetric with respect to relabelling M₀ and M₁. There are many alternative loss functions without this property; such functions might be useful for certain applications (e.g. when it is known a priori that the data will prefer one model far more than another).

For a network to calculate $p(M_1|x)$, such an 'asymmetric' loss might have the form:

$$\begin{align} \mathcal{V}\left(f\left(x\right), m\right) = m \mathcal{A} \left(f\left(x\right)\right) + \left(1-m\right) \mathcal{B} \left(f\left(x\right)\right) , \end{align} \tag{ 31 }$$

for a choice of functions $\mathcal{A}$ and $\mathcal{B}$ where $\mathcal{A} (f(x)) \neq \mathcal{B} (1 - f(x))$. Some of these losses have been explored for computing likelihood-ratios using classifiers in Rizvi et al (2023).

Even if theoretically there exists an optimal function $f^{\;*}$ that minimizes the global optimization objective equation (6), in practice three factors (typical of neural methods) prevent us from finding it exactly. First, we must approximate the space of functions with a set of neural network architectures. Second, we have access to only a finite number of samples from p(x) and hence we only have a (Monte Carlo) approximation to the optimization objective, via:

$$\begin{align} I\left[\,f\,\right] \approx \frac{1}{N} \sum_{j = 1}^N \ \mathcal{V}\left(f\left(x_j\right), m_j\right). \end{align} \tag{ 32 }$$

Third, we have no means of guaranteeing that a local optimum discovered during training is in fact a global optimum. It is, therefore, inevitable that there will be some error associated with our estimate of the Bayes factor.

In the next section 4, we will demonstrate that the Evidence Network does indeed recover robust estimates of the Bayes factor. For such synthetic tests, we can fully validate our approach by choosing a model for which the true Bayes factors can be calculated analytically (i.e. in closed-form).

In actual applications, the likelihood and or prior are intractable or implicitly defined only through labeled training data (e.g. consisting only of data realizations and their associated model labels from simulations or a generative model). In this setting a closed-form solution will neither be available nor necessary since the true Bayes factors or model posteriors are not required to train an Evidence Network. Despite this, it is still necessary to validate the trained network. We will now describe a method for validating the Bayes factor computation by the Evidence Network that is generally applicable and requires only labeled training data of the same type that is needed to train the Evidence Network itself.

The solution to this unique validation problem is blind coverage testing. We can compare the estimate of $p(M_1 | x)$ to the fraction of validation data with the model label M₁.

Let us assume that the trained neural network is optimal, such that the network recovers the function $f^{\;*}$. For example, if the network is trained using the l-POP-Exponential Loss (equation (28)), for equal model priors we expect that:

$$\begin{align} \mathcal{J}_\alpha \left(f^{\;*}\left(x_O\right)\right) = \ \log K. \end{align} \tag{ 33 }$$

For simplicity, consider the case of only two models M₁ and M₀ (this can be easily generalized for any number of models). We can then rearrange the previous expression to give the model posterior for M₁:

$$\begin{align} p\left(M_1 | \ x\right) = \frac{\exp \left(\mathcal{J}_\alpha \left(f^{\;*}\left(x_O\right)\right) \right)}{1+ \exp \left(\mathcal{J}_\alpha \left(f^{\;*}\left(x_O\right)\right) \right)} . \end{align} \tag{ 34 }$$

Such an expression can be constructed for any of the losses in table 1 that estimate some function of K.

Fix P and let ε > 0. Consider validation data $x^\textrm{v} \sim p(x)$ for which the posterior probability $p(M_1 |x^\textrm{v})$ lies in $\left[ P, P +\epsilon \right]$. As $\epsilon \rightarrow 0$, the expected fraction of such data having label M₁ converges to P. By comparing the true fraction of model labels in some small probability range with that expected from the estimated K, we can validate the trained network without knowing the true model posterior probabilities.

Put another way, consider a set of N validation data samples that all have the same posterior probability $p(M_1 |x)$. In this set, if the number of data samples that are drawn from model M₁ is n₁, then the expected fraction of data from model M₁ is given by $\langle \frac{n_{1}}{N} \rangle = p(M_1 |x)$, following the usual definition of posterior probability. In practice, no two samples in the validation data have exactly equal posterior probabilities, so we use small bins of probability to calculate this fraction (corresponding to a small ε). To quantify the statistical significance of the scatter in the probability bins, we can use a binomial distribution to model the count statistics (see section 5 for an example).

This procedure resembles the re-calibration of class probabilities in classification tasks (Zadrozny and Elkan 2002, Niculescu-Mizil and Carua 2005). In this work, it is used as a blind validation of our Bayes factor estimates for Bayesian model comparison.

We construct two generative models to sample time series data x. These are defined such that the model evidences and associated Bayes factors can be calculated analytically using a closed-form expression. In section 4.2, this is compared to the estimated log Bayes factors from our simple Evidence Network.

We compare the Evidence Network to likely alternative methods. Even if the true likelihood/prior were known, with the use of the most advanced nested sampling methods (Handley et al 2015) the problem can still be too arduous for high parameter-space dimensionality $\mathrm{dim}(\theta) \gtrsim 10^2$; we demonstrate this in section 4.4. Assuming that we only have data samples and model labels, we would have no knowledge of the distributions $p(x | \theta)$ or $p(\theta)$ to evaluate the marginal likelihood (equation (3)). As previously discussed, one could try to estimate those probability densities using implicit inference, but with high data dimensionality $\mathrm{dim}(x) \gtrsim 10^2$ or high parameter-space dimensionality $\mathrm{dim}(\theta) \gtrsim 10^2$, density estimation becomes intractable or prone to significant error; we show this in section 4.5.

We construct two nested models that generate time series data with heteroschedastic random noise. Both models are linear Gaussian models with a linear operation that is non-trivial with respect to time.

The prior for the model parameters θ is a simple normal distribution for each parameter θ_i :

$$\begin{align} p\left(\theta_i\right) = \mathcal{N}\left(\theta; 0,1\right) . \end{align} \tag{ 35 }$$

Per data (vector) realization x, each element x_j is given by:

$$\begin{align} x_j = \sum\limits_{i} A_{ji} \theta_j + n_j = \mu^x_j + n_j , \end{align} \tag{ 36 }$$

where the heteroscedastic noise is drawn from a multivariate normal distribution, such that:

$$\begin{align} p\left(n\right) = \mathcal{N} \left(n;0, {\Sigma}\right). \end{align} \tag{ 37 }$$

The noise covariance Σ is a diagonal matrix with diagonal elements ${\Sigma}_{k k}$ given by:

$$\begin{align} \sqrt{{\Sigma}_{k k}} = \sqrt{\frac{N}{100}} \left( \left(k+\frac{5}{2} \right) \frac{8}{5}\right)^2 . \end{align} \tag{ 38 }$$

For the first model, M₁, the linear operation is given by:

$$\begin{align} A_{ji} = \cos \left(\left(i-\frac{1}{2}\right)t_j \right) \end{align} \tag{ 39 }$$

for i > 0. For i = 0, we have

$$\begin{align} \ A_{j0} = 2 t_j . \end{align} \tag{ 40 }$$

The i = 0 term corresponds to a linear growth over time. For the second model, M₀, this first term is removed, so there is no linear growth term over time. This is equivalent to setting $\theta_0 = 0$.

An example data sample is shown in figure 1.

Zoom In Zoom Out Reset image size

In both of these cases we can evaluate the marginal evidence integral (equation (3)) analytically to recover a closed-form solution for the Bayesian model evidence, which is equivalent to evaluating a normal distribution:

$$\begin{align} p\left(x | \mathcal{M}\right) = Z\left(x| \mathcal{M}\right) = \mathcal{N}\left(x ; \mu_Z, C_Z\right), \end{align} \tag{ 41 }$$

where

$$\begin{align} \mu_Z = A \cdot \mu^x , \end{align} \tag{ 42 }$$

and

$$\begin{align} C_Z = \Sigma + A A^T . \end{align} \tag{ 43 }$$

This can be simply evaluated for each model to calculate the ratio to give the Bayes factor K:

$$\begin{align} K = \frac{Z\left(x | M_1\right)}{Z\left(x | M_0\right)} . \end{align} \tag{ 44 }$$

Using data drawn from the generative model described in the previous Subsection, we train the Evidence Network using the l-POP-Exponential Loss with α = 2 to recover a simple transformation of $\log K$ (see table 1). This choice of loss function is our default recommendation:

$$\begin{align} {\mathrm{l}}\text{-}{\mathrm{POP}}\text{-}{\mathrm{Exponential} \mathrm\,{Loss}: } e^{\mathcal{J}_{\alpha = 2}\left(f\left(x\right)\right)} = e^{\left(\frac{1}{2} - m \right)\left(f\left(x\right)+f\left(x\right)|f\left(x\right)|\right)} \implies f^{\;*}\left(x_O\right)\left(1 + |f^{\;*}\left(x_O\right)|\right) = \log K . \end{align} \tag{ 45 }$$

Using the model labelling described above, a Bayes factor greater than one corresponds to a preference for model M₁ (i.e. $\theta_0 = 0$ is not preferred in this case).

We find that taking the average output of an ensemble of networks, each with independent random weight initialization, significantly improves the accuracy of our log Bayes factor estimation. For these results, each network has the identical architecture.

We have chosen a simple network with little tuning that uses 6 dense layers (see appendix for details). The network was purposefully chosen to be simple to demonstrate that the network does not need to be tuned for the Bayes factor estimates to be reliable.

Though we have chosen a simple dense architecture, Evidence Networks can use any type of network architecture, provided they use the correct optimization objective and are validated using our proposed coverage test. Depending on the application, they could use geometric architectures (e.g. deep sets, graph networks; Battagli et al 2018), recurrent neural networks (e.g. long short-term memory—LSTM; Hochreiter and Schmidhuber 1997), or any appropriate architecture.

Our main result uses the generative model described above with $\mathrm{dim}(\theta) = 10^2$ and $\mathrm{dim}(x) = 10^2$. As the data samples for this generative model are symmetric under a sign flip, we make use of this standard data augmentation during training.

The left panel of figure 2 shows our main result with the ensemble Evidence Network. We recover estimates of the Bayes factor over many orders of magnitude. When compared with the analytic (true) Bayes factors, we calculate the root-mean-square error (RMSE) and find a $\log K$ error ${\approx}0.02$.

Zoom In Zoom Out Reset image size

This is an excellent estimation of the Bayes factor for such a high dimensional problem. The high dimensionality of the data and parameters render this task either intractable or error-prone for alternative methods as we demonstrate in the following subsections.

For most applications, the true Bayes factor value will typically not be available for validation; our training data consists only of data realizations with their associated model label. We can, however, leverage our blind coverage test to validate the model posterior probabilities are recovered. This is shown in the right panel of figure 2 for the same trained network and data.

The plotted error-bars in the right panel of figure 2 are the expected standard error for binomial distributed samples: $\sigma_\textrm{err} = \sqrt{p(M_1|x_O) (1 - p(M_1|x_O))/n_p}$, where n_p are the number of data samples in a given p probability bin. We find the measured standard deviation of the residual error (between the expected and the predicted probabilities) does indeed match σ_err, further validating the method.

For a general problem where the true Bayes factors are not known, the result of such a blind coverage test can be used to construct error bounds for the Bayes factors estimates.

Figure 3 shows different choice of loss functions. The left panel does not use a loss function to estimate the Bayes factor directly, but uses the Cross Entropy Loss to individually estimate $p(M_1 | x_O)$ and $p(M_0 | x_O)$, whose ratio gives the Bayes factor K. This approach leads to systematic inaccuracies for moderately high values of K. These inaccuracies are ameliorated with the choice of the Exponential Loss (center panel), from which the Evidence Network directly estimates $\log K$. Using the l-POP-Exponential Loss (right panel) the systematic error at high K is no longer apparent.

Zoom In Zoom Out Reset image size

In this subsection we compare the Evidence Network with nested sampling; we assume the correct likelihood for nested sampling, while restricting the Evidence Network to the likelihood-free problem, and find that the Evidence Network is still able to outperform nested sampling for high-dimensional parameter spaces.

Let us assume we know the likelihoods, $p(x | \theta, M_1)$ and $p(x | \theta, M_0)$, and priors, $p( \theta | M_1$ and $p( \theta | M_0)$. We assume that we are able to either evaluate a given likelihood, $\mathcal{L}_1(\theta) = p(x_O | \theta, M_1)$, or sample data realizations from it, $x_i \sim p(x | \theta, M_1)$. For a Gaussian likelihood, these two operations are of equal computational complexity.

We use the PolyChord (Handley et al 2015) code to perform the nested sampling; this uses an advanced algorithm that combines both cluster detection for multi-modal distributions and slice sampling using affine 'whitening' transformations. We use the default number of live points $= 25\times \textrm{dim} (\theta)$, which then leads to a certain number of $p(x|\theta)$ evaluations. We use a maximum cut-off of 10⁸ likelihood evaluations for reasons of computation time.

We choose the number of $p(x|\theta)$ samples for the Evidence Network to match the same polynomial scaling with $\textrm{dim} (\theta)$ as PolyChord evaluation, but with only approximately 1 per cent of that of PolyChord. Despite significantly fewer $p(x|\theta)$ samples than nested sampling $p(x|\theta)$ evaluations (figure 4, left panel), we consistently recover more accurate estimates of $\log K$ with the Evidence Network (figure 4, right panel).

Zoom In Zoom Out Reset image size

This is an indication that it may be more efficient to train an Evidence Network to estimate Bayes Factors even when the likelihood is known. For many closed-form likelihoods, it is a relatively simple code change to replace the $p(x|\theta)$ evaluation with a $p(x|\theta)$ sample of x.

Let us now consider alternative approaches that are available for likelihood-free (e.g. simulation-based) inference problems.

As previously discussed in section 1.4, we could first learn the form of $p(x | \theta )$ from the training data; this is known as NLE. We would still need to perform the high-dimensional marginal evidence integration and using nested sampling (or a similar method). This step is unnecessary: the intermediate density estimation step can only induce potential errors and we nevertheless return to the difficult procedure discussed in section 4.4.

An alternative density estimation approach would be to learn the form of $p(x | \mathcal{M})$ from the training data. So far, this method has not been widely considered for scientific model comparison, even though it would also avoid marginal evidence integration. Therefore, we compare Evidence Networks with this density estimation approach.

We use an ensemble of neural density estimators to estimate $p(x | M_1)$ and $p(x | M_0)$ and use an ensemble of neural networks (matching the set-up used so far) to form the Evidence Network to estimate $\log K$. In both cases we use 10⁶ samples of training data x. Density estimation suffers particularly badly from increasing dimensionality. We will, therefore, reduce from a dimensionality of $\textrm{dim} ( x ) = 100$ (as in figure 2) to a more manageable $\textrm{dim} (x) = \textrm{dim} (\theta) = 20$ for this comparison.

We use Normalizing Flows for neural density estimation; we use Neural Spline Flows (Durkan et al 2019), which are a popular and flexible neural density estimation method. As an implementation, we use the PZFlow package (Crenshaw et al 2021) with the default settings.

Figure 5 shows the comparison between the two methods. The left panel shows the analytic $\log K$ value against the estimated value using both the Evidence Network and the Normalizing flow ratio (using $p(x | M_1)/p(x | M_0)$). The relatively large residual scatter for the Normalizing Flow can also be seen in the right panel, which shows the histogram of the residual errors.

Zoom In Zoom Out Reset image size

The RMSE (error) in the $\log K$ using the Normalizing Flow is more than a factor of 10 larger than the Evidence Network. The reason for this is twofold: (i) density estimation suffers exponentially from increasing dimensionality (which is even a problem for this restricted 20 dimensional data), and (ii) the errors in each $p(x | \mathcal{M})$ combine in the ratio that gives $\log K$.

7.1.1. Absolute evidence

7.1.2. Connection to frequentist hypothesis testing

7.1.3. Using Evidence Networks to do posterior predictive tests

While the interpretation of the Bayes factor is clear, the strong dependence on prior choice can be seen as a weakness. A common alternative is posterior predictive testing (PPT, see Gelman et al 2013 for overview). PPT partitions the data x into subsets, x₀ and x₁, say, and evaluates the relative probabilities of two models having predicted x₁ conditional on x₀. This typically involves evaluating $p(\theta|x_0)$ (posterior inference), computing the integral

$$\begin{align} p\left(x_1|M_i,x_0\right) = \int p\left(x_1|\theta,x_0\right) p\left(\theta|x_0\right) \ {\mathrm{d}}\theta \end{align} \tag{ 48 }$$

for M₀ and M₁, and then computing the posterior odds ratio

$$\begin{align} K_\mathrm{PPT} = \frac{p\left(x_1|M_1,x_0\right)}{p\left(x_1|M_0,x_0\right)}. \end{align} \tag{ 49 }$$

The difficulty is that all the challenges of computing the evidence ratio with conventional methods, which motivated the present study, also apply to $K_\mathrm{PPT}$.

But it is straightforward estimate $K_\mathrm{PPT}$ using our methodology. The simplest way to see this is to write

$$\begin{align} K_\mathrm{PPT} = \frac{p\left(x_1|M_1,x_0\right)}{p\left(x_1|M_0,x_0\right)} = \frac{p\left(x_0,x_1|M_1\right)}{p\left(x_0,x_1|M_0\right)} \frac{p\left(x_0|M_0\right)}{p\left(x_0|M_1\right)} = {K\left(x\right)}/{K\left(x_0\right)}, \end{align} \tag{ 50 }$$

which can be computed simply by training two Evidence Networks, one to compute the Bayes factor for the full data x and one for the subset x₀.

In conclusion, we have presented a novel approach for Bayesian model comparison with a wide range of applications; it overcomes the limitations of even state-of-the-art classical computational methods.