Targeting Bayes factors with direct-path non-equilibrium thermodynamic integration

Thermodynamic integration (TI) for computing marginal likelihoods is based on an inverse annealing path from the prior to the posterior distribution. In many cases, the resulting estimator suffers from high variability, which particularly stems from the prior regime. When comparing complex models with differences in a comparatively small number of parameters, intrinsic errors from sampling fluctuations may outweigh the differences in the log marginal likelihood estimates. In the present article, we propose a TI scheme that directly targets the log Bayes factor. The method is based on a modified annealing path between the posterior distributions of the two models compared, which systematically avoids the high variance prior regime. We combine this scheme with the concept of non-equilibrium TI to minimise discretisation errors from numerical integration. Results obtained on Bayesian regression models applied to standard benchmark data, and a complex hierarchical model applied to biopathway inference, demonstrate a significant reduction in estimator variance over state-of-the-art TI methods.


Introduction
A central quantity in Bayesian statistics is the marginal likelihood where D are the data, and M represents a given statistical model with parameter vector θ . The difficulty in practically computing the marginal likelihood is exemplified by considering the Monte Carlo sum where {θ i } is an iid sample from p(θ|M). Under fairly general regularity conditions the estimator X converges almost surely to p(D|M), by the strong law of large numbers, and is asymptotically efficient, with asymptotic variance C/ √ N (where N is the size of D), by the central limit theorem. However, even for modestly complex systems, the constant in the numerator, C, can reach exorbitant magnitudes, rendering the scheme not viable for practical applications. The practical shortcomings of a variety of alternative numerical methods, like the harmonic mean estimator (Gelfand and Dey 1994), bridge sampling (Gelman and Meng 1998), or Chib's method (Chib and Jeliazkov 2001), have been discussed in the statistics and machine learning literature (e.g. Murphy (2012)). The most widely used and robust method appears to be thermodynamic integration (TI). This method was originally proposed by Kirkwood (1935) and further developed in statistical physics for the mathematically equivalent problem of computing free energies; see e.g. Schlitter (1991) and Schlitter and Husmeier (1992). Gelman and Meng (1998) adapted TI to the computation of marginal likelihoods, Lartillot and Philippe (2006) demonstrated the application of TI to complex systems, and Friel and Pettitt (2008) and Calderhead and Girolami (2009) popularised TI more widely in the statistics community by demonstrating a computationally powerful combination with parallel tempering (Earl and Deem 2005).
Thermodynamic integration is based on an integral of the expected log likelihood along an inverse annealing path from the prior to the posterior distribution. The resulting estimator typically suffers from high variability, which particularly stems from the parameter prior regime. When comparing complex models with differences in a comparatively small number of parameters, these intrinsic errors from sampling fluctuations may outweigh the differences in the log marginal likelihood estimates. The objective of the present study is to explore the scope for variance reduction by directly targeting the log Bayes factor via a modified transition path between the two models such that the high-variance prior domain is avoided. This idea is not new. In statistical physics it is well known (Schlitter 1991;Schlitter and Husmeier 1992) that applying TI to the computation of a reaction free energy, which is mathematically equivalent to the log Bayes factor, is computationally more efficient than the separate computation of the standard free energies for the two reaction states involved (educt versus product states); the latter is mathematically equivalent to the difference of the log marginal likelihoods of two statistical models to be compared. Also in the statistics literature, the direct targeting of the log Bayes factor has been discussed before. For instance, path sampling (Gelman and Meng 1998) and annealed importance sampling (Neal 2001) have been conceived in a way to allow the direct computation of the ratio of two partition functions, Z 1 and Z 2 , associated with two models M 1 and M 2 . However, in the work of Neal (2001), Z 1 is set to the normalisation factor of the prior distribution, and the method thus reduces to the computation of the log marginal likelihood. 1 Gelman and Meng (1998) do consider a direct comparison between two alternative models: a homoscedastic versus a heteroscedastic linear regression model. Rather than computing the Bayes factor, the authors apply their path sampling approach to infer the posterior distribution of the entire spectrum of intermediate models. While this is a more ambitious approach than model selection with Bayes factors, it will be computationally onerous beyond the one-dimensional regime considered in their example.
To the best of our knowledge, the present article presents the first systematic study of the variance reduction that can be achieved with a thermodynamic integration path that directly targets the log Bayes factor by transiting between the posterior distributions of the two models involved. The mathematical exposition and implementation of this scheme is combined with a comprehensive comparative performance assessment based on a set of standard benchmark data to quantify the improvement in variance reduction, accuracy and computationally efficiency that can be achieved over stateof-the-art established TI methods, in particular the recent improvement proposed by Friel et al. (2014).
This article is organised as follows. In Sect. 2 we provide a brief rationale for targeting Bayes factors directly rather than indirectly via the marginal likelihood. Section 3.1 reviews standard thermodynamic integration. In Sect. 3.2 we discuss a modified numerical integration and sampling scheme from statistical physics, termed non-equilibrium TI (NETI), to reduce numerical discretisation errors. Section 3.3 describes NETI-DIFF, the proposed new TI scheme along an alternative integration path between two posterior distributions. Sections 3.4, 3.5 describe practical numerical implementations based on Metropolis-Hastings and Gibbs sampling, and Sect. 3.6 proposes a new improved inverse temperature ladder. Section 4 provides an overview of a set of benchmark problems on which we have evaluated the methods, and Sect. 5 presents our empirical findings. We conclude this article in Sect. 6 with a discussion, a comparison with the controlled thermodynamic integral of Oates et al. (2016), and an outlook on future work.

Rationale
Consider two alternative models, M 1 and M 2 , and define E i = − log p(D|θ, M i ), the negative log likelihood of model i. Further, define the log likelihood ratio ΔE = 720 M. Grzegorczyk et al.
E 2 − E 1 , the negative unnormalised log posteriorẼ i = E i +log p(θ |M i ), the negative log posterior ratio ΔẼ =Ẽ 2 −Ẽ 1 , and let . . . i denote the posterior average with respect to the posterior distribution (θ |D, M i ). We can then adapt Jarzynski's theorem from statistical physics (Híjar and Zárate 2010) to show that A proof is given in the Appendix. In real applications with non-trivial models, the negative log likelihood is typically in the order of a two to three digit figure, which when put into the argument of the exponential function will lead to an astronomically large number. An estimator aiming to approximate p(D|M i ) = exp(E i [θ]) −1 i from a limited sample drawn from the posterior distribution will inevitably suffer from substantial variation. For nested models or models with sufficient parameter overlap, on the other hand, ΔẼ(θ) will typically be small, |ΔẼ(θ)| min{|E 1 (θ )|, |E 2 (θ )|}. We can therefore reduce the intrinsic estimation uncertainty considerably by computing the Bayes factor directly rather than indirectly via two separate marginal likelihood estimations.

Thermodynamic integration for marginal likelihoods
Thermodynamic integration is based on an inverse annealing path from the prior to the posterior distribution, and computing the expectation of the log likelihood with respect to the following annealed posterior distributions at inverse temperatures τ ∈ [0, 1]: Taking the derivative of log Z (D|τ, M) gives: From Eq. (5) we get: This one-dimensional integral can be solved numerically, e.g. with the trapezoid rule: Some care has to be taken with respect to the choice of discretisation points τ k , k = {0, 1, 2, . . . , K }, as the major contributions to the integral usually come from a small region around τ → 0. This motivates the form for α > 1. Theoretical results for the optimal choice of α can be found in Schlitter (1991), but require knowledge that is usually not available in practice (like the functional dependence of E τ [log p(D|θ, M)] on τ ). In practice, α = 5 is widely used, as e.g. in Friel et al. (2014), and we have used this value in the present study. A potentially numerically more stable alternative was proposed in Friel et al. (2014). The authors show that: where V τ (.) is the variance w.r.t. the power posterior in Eq. (4  Friel et al. (2014) then employ the corrected trapezoid rule 2 to compute each subintegral

Nonequilibrium thermodynamic integration
The computation of the expectation values E τ k log p(D|θ, M) is expensive and limits the number of discretisation points K that can be practically applied. An alternative scheme we use in the present work is to approximate where θ (τ ) is a single draw from the power posterior defined in Eq. (4), and 0 = τ 1 < τ 2 < · · · < τ K = 1. The computational resources gained are used to choose K orders of magnitude larger than in equilibrium TI, 3 with the implication that (τ k − τ k−1 ) → 0 and discretisation errors in numerical integration are avoided. This scheme was originally proposed in statistical physics (Schlitter and Husmeier 1992) under the name non-equilibrium thermodynamic integration (NETI), and is conceptionally similar to annealed importance sampling (Neal 2001). The underlying rationale is as follows: rather than use the computational resources for the computation of the expectation value at a limited number of discretisation points-and incur discretisation errorsspread the computational resources over the whole "temperature" range and use as fine a discretisation as possible. This avoids the problem that had to be addressed in Friel et al. (2014): how to select the inverse temperatures and minimise the numerical integration error in standard TI. The price to pay is a relaxation error as a consequence of the non-equilibrium nature of the method, as discussed by Schlitter and Husmeier (1992). The authors proposed a scheme for correcting this relaxation error, by running simulations over different simulation lengths N iter , regressing the estimates against an approximate upper bound on the relaxation error R, and then extrapolating for R → 0.
In preliminary investigations omitted from the present article, we found that a single simulation with an increased value of N iter matching the total computational costs of the extrapolation scheme achieved similar results, and we used this conceptionally simpler approach in all our studies. 4

Novel thermodynamic integration for Bayes factors
When comparing two models, we are typically interested in the Bayes factor p(D|M 2 )/ p(D|M 1 ). The standard approach is to apply thermodynamic integration to both models M 1 and M 2 separately, by independently inversely annealing the prior distributions to the respective posterior distributions. This approach ignores the fact that both models usually have many aspects in common and share certain parameters. This applies particularly to nested models, where all the parameters of the less complex model are also included in the more complex model. One would expect to reduce the estimation uncertainty by following a direct transition path from the posterior distribution of the less complex model to that of the more complex model, rather than transiting through the uninformative prior distribution twice. Consider two models M 1 and M 2 with joint parameter vector θ and a joint parameter prior p(θ |M 1 , M 2 ) defined such that it reduces to the parameter priors for the separate models by marginalisation: (12) where M 2 /M 1 is the subset of parameters contained in model M 2 , but not in model M 1 , and M 1 /M 2 is the subset of parameters contained in model M 1 , but not in model M 2 . A mathematically more accurate notation would split θ into three subsets, θ = {θ 1 , θ 2 , θ 12 } such that θ 1 ∈ M 1 /M 2 , θ 2 ∈ M 2 /M 1 and θ 12 ∈ M 1 ∩ M 2 . Eq. (12) implies that p(θ|M 1 ) = p(θ 1 , θ 12 |M 1 ) and p(θ |M 2 ) = p(θ 2 , θ 12 |M 2 ). For that reason we can use a mathematically redundant but less opaque notation that does not make the partition θ = {θ 1 , θ 2 , θ 12 } explicit. Define the tempered posterior distribution where From Eq. (12) we get: Taking the derivative of the partition function in Eq. (14) gives: Combining Eqs. (15,16) gives the following thermodynamic integral for the log Bayes factor: Again, we follow the idea of non-equilibrium thermodynamic integration and make the approximation where θ (τ ) is a single draw from the tempered posterior distribution defined in Eq. (13), 0 = τ 1 < τ 2 < · · · < τ K = 1, K 1, and (τ k − τ k−1 ) 1. In comparison with statistical physics, the proposed scheme corresponds to the direct computation of a free energy difference (Schlitter 1991;Schlitter and Husmeier 1992), which is more efficient, in terms of reduced estimation variance for given computational costs, than computing the difference of two separately computed standard free energies. The analogy from classical statistics is model comparison via a paired test, which is known to have higher power than an unpaired test.
In what follows, we refer to the estimator defined by Eq. (18) as NETI-DIFF. We describe how to compute the variance of this estimator in the Appendix 7.2.

Metropolis-Hastings scheme
The implementation of a Metropolis-Hastings scheme to target the distribution in (13) is straightforward. Given the current parameters θ, sample new parametersθ from a proposal distribution q(θ |θ), and accept these new parameters with the following acceptance probability: Otherwise, setθ = θ, and follow this scheme iteratively.

Gibbs sampling for linear models
Consider a standard linear model with parameter vector θ, design matrix D, and prior distribution The data, D = {y 1 , . . . , y T } or Y = (y 1 , . . . , y T ) T , are assumed to be obtained under the assumption of independent and identically distributed normal noise, with variance σ 2 : p(y|θ , σ 2 ) = N (Dθ , σ 2 I) We want to compare two competing models M 1 and M 2 , represented by two alternative design matrices D (1) and D (2) : For notational compactness we choose a representation that leaves the dimension of θ invariant with respect to changing model dimensions by padding obsolete entries in the design matrix with zeros. For instance, to compare the models M 1 :y = θ 1 x 1 + θ 2 x 2 , and M 2 :y = θ 1 x 1 + θ 3 x 3 + θ 4 x 4 based on a data set of n observations {y t , x 1,t , x 2,t , x 3,t , x 4,t }, t = 1, . . . , n, we get the following design matrices: where the factor C(y) does not depend on θ . Comparing this with the identity we get: where Hence, we can directly sample θ from the tempered conditional distributions in a Gibbs sampling scheme without having to resort to Metropolis-Hastings. For linear models where the variance σ 2 is not known and has to be sampled from the tempered posterior distribution too, we refer to Appendix 7.6.

Sigmoid inverse temperature ladder
Given a single model M, conventional TI follows an inverse annealing path from the prior p(θ|M) to the posterior p(θ|M, D), symbolically p(θ |M) → p(θ |M, D). Unlike TI, NETI-DIFF is based on a direct transition from the posterior of one model M 1 to the posterior of another model M 2 , p(θ|M 1 , D) → p(θ |M 2 , D). For nested models, e.g. M 1 ⊂ M 2 , we start at the less complex model M 1 and move towards the more complex model M 2 , e.g. using the power-law inverse temperature ladder, defined in Eq. (8). For a power α > 1 the distances τ i+1 − τ i between neighbouring discretisation points τ k and τ k+1 increase in k and the discretisation points will be concentrated around the nested model, M 1 (τ = 0), and fewer points will be set near M 2 (τ = 1). However, in many applications non-nested models have to be compared, and it is then not clear which of the two models should be used as starting point. Imbalances can be avoided by choosing a sigmoid inverse temperature ladder, such that the discretisation points are mirrored at the midpoint τ = 0.5 of the interval [0, 1]. Every discretisation point τ < 0.5 closer to M 1 then has its counterpart τ = 1 − τ with the same distance τ to M 2 , and vice-versa.
To obtain a sigmoid inverse temperature ladder for NETI-DIFF we apply the following procedure. We first specify 50% of the discretisation points τ 1 < · · · < τ N iter 2 within the interval [0, 0.5], and then we mirror the ladder at the midpoint τ = 0.5. 5 This yields the remaining 50% of the discretisation points, τ N iter . . , N iter 2 ). As we want the first 50% of the discretisation points to get as close as possible to the midpoint τ = 0.5 subject to a power law with power α, we determine the minimal integer N such that The solution is:

Benchmark problems and data
We have evaluated the proposed method on four benchmark data sets. Given data D the goal is to estimate the log Bayes factor B between two models M 1 and M 2 . We assume the models to be equally likely a priori, p(M 1 ) = p(M 2 ), so that the Bayes factor is the ratio of marginal likelihoods: For nonuniform prior distributions, p(M 1 ) = p(M 2 ), it is straightforward to add the correction factor log p(M 2 )/ p(M 1 ), which is computationally cheap compared to the marginal likelihood ratio. For method evaluation, we need to compare with a ground truth. For a linear model, we have a proper ground truth, as the Bayes factor can be computed analytically. This applies to the Radiata pine data (Sect. 4.1) and the Radiocarbon data (Sect. 4.3). For the Pima Indian data (Sect. 4.2), we use a generalised linear model, and for the biopathway data (Sect. 4.4), we use a nonlinear model. In these cases, a closed-form solution of the Bayes factor does not exist. For the Pima Indian data, we follow the 5 For uneven N iter , we fix one point at τ = 0.5 and apply the procedure to the remaining N iter − 1 points. method suggested in Friel et al. (2014) and use the numerical result from a very long MCMC run as an approximate gold standard. For the biopathway data, we use the knowledge of the true interaction structure of the system as a surrogate gold standard and assess the performance in terms of network reconstruction accuracy. We think this provides an adequate balance between using linear models, for which a strong ground truth exists, and generalised linear/non-linear models, for which a strong ground truth is intrinsically unavailable, and a weaker surrogate ground truth has to be used instead.

Radiata pine
The Radiata pine data have been used in Friel et al. (2014) and were originally published in Williams (1959). Like Friel et al. (2014) we focus on the log Bayes factor between two competing non-nested linear regression models for explaining the 'maximum compression strength' y of n = 42 Radiata pine specimens. Both linear models contain an intercept and one single covariate. The first model (M 1 ) uses the 'density' x 1 and the second model (M 2 ) the 'adjusted density' x 2 of the specimen. After standardizing the observation vectors x 1 and x 2 of the two covariates to mean 0, the likelihood of model M k (k = 1, 2) is: where y is the vector of the observed 'maximum compression strengths', D (k) = (1, x k ) is the n-by-2 design matrix and θ (k) is the 2-dimensional vector of regression coefficients of model M k . Both models share the intercept parameter θ 0 , but differ w.r.t. the second parameter, i.e. θ (k) = (θ 0 , θ k ) . For comparability we use exactly the same Bayesian model as in Friel et al. (2014), where an inverse Gamma prior is imposed on the noise variance: p(σ −2 ) = G AM(3, 2 · 300 2 ) and Gaussian priors are used for the regression coefficient vectors: 6 This is a model with fully conjugate priors, so that the marginal likelihoods p(y|M k ) can be computed in closed form (Friel et al. 2014). With Eq. (29) we obtain for the log Bayes factor B(M 1 , M 2 ) = 8.8571. Like Friel et al. (2014) we apply Gibbs sampling and re-sample the model parameters iteratively from their full conditional distributions: p(σ 2 |y, θ (k) ) and p(θ (k) |y, σ 2 ).

Pima Indians
The Pima Indians data have also been used in Friel et al. (2014) and were originally published in Smith et al. (1988). Like Friel et al. (2014) we focus on the log Bayes factor between two nested logistic regression models for explaining the binary 'diabetes disease status' y of n = 532 female Pima Indians. The first model (M 1 ) contains an intercept and 4 covariates, namely 'the number of pregnancies', 'the plasma glucose concentration', 'the body mass index', and 'the diabetes pedigree function', while the second model (M 2 ) extends model M 1 by including one additional covariable 'age'. After standardizing all covariates to mean 0 and variance 1, the likelihood of model M k (k = 1, 2) is: where the i-th element of y, y i ∈ {0, 1}, is the diabetes status of female i, x i,k is the corresponding vector of covariates, including an initial 1 for the intercept, and θ (k) is the vector of regression coefficients of dimension m = 5 (M 1 ) or m = 6 (M 2 ). Again we follow Friel et al. (2014) and impose the following Gaussian priors on the regression coefficient vectors: gives rather uninformative priors. For the logistic regression neither the marginal likelihoods nor the full conditional distributions can be computed in closed form. We therefore use the Metropolis Hastings based Markov chain Monte Carlo (MCMC) sampling scheme from Friel et al. (2014), which employs the following proposal mechanism: In each iteration a new candidate regression coefficient vector is obtained by adding a sample u from an m-dimensional multivariate Gaussian distribution to the current vector θ (k) . The Gaussian distribution of u has a zero mean vector and a diagonal covariance matrix, whose diagonal entries d 1 , . . . , d m depend on the inverse temperature τ ∈ [0, 1] of the power posterior. For the TI approaches we set: d i = min{0.01τ −1 , 100}, as in Friel et al. (2014). For the proposed NETI-DIFF approach we use d 6 = min{0.01τ −1 , 100}, while we fix the first five diagonal entries d 1 , . . . , d 5 = 0.01. This modification is required, as the first five regression coefficients appear in both models M 1 and M 2 . That is, they effectively appear constantly with inverse temperature τ = 1 throughout NETI-DIFF simulations. The marginal likelihoods cannot be computed in closed-form. We therefore use those values reported in Friel et al. (2014), which were obtained from long TI simulations, as gold-standard: log{ p(y|M 1 )} = −257.2342 and log{ p(y|M 2 )} = −259.8519. Equation (29) yields the log Bayes factor: B(M 1 , M 2 ) = −2.6177.

Radiocarbon dating
We use the Radiocarbon data from Pearson and Qua (1993) to compute the Bayes factors among 10 nested linear regression models. For predicting the 'true calendar age' y of n = 343 Irish oaks from one single covariable: 'the Radiocarbon dating process' x, we fit polynomial calibration curves M i (i = 1, . . . , 10) of the following type: The likelihood of model M i is then where y is the vector of calendar ages, θ (i) = (θ 0 , θ 1 , . . . , θ i ) is the vector of regression coefficients, and D (i) is the n-by-(i + 1) design matrix. The first column of the design matrix consists entirely of ones (for the intercept), and the subsequent columns are built from the observation vector x, D (i) = (1, x 1 , . . . , x i ), where x j denotes an element-wise power operation on x. We impose conjugate priors on the parameters. For σ 2 we use an inverse Gamma distribution: p(σ −2 ) = G AM( a 2 , b 2 ), and on θ (i) we impose Gaussian priors: For fixed hyperparameters a, b, and δ 2 the marginal likelihood for a model M with design matrix D is given by: so that the log Bayes factors B(M i , M l ) for two models M i and M l can be computed in closed form with Eq. (29). For the Radiocarbon data we fix a = b = 0.2, δ 2 = 1, and we sample the parameters iteratively from their conditional distributions p(σ 2 |y, θ (i) ) and p(θ (i) |y, σ 2 ) with Gibbs sampling.

Biopathway
The objective of the last application is model selection with respect to two alternative candidate interaction structures of ten genes in the circadian gene regulatory network of Arabidopsis thaliana, shown in Fig 2017)). Let x i (t) denote the mRNA concentration of gene i at time t, and π i the set of its regulators. For instance, in the gene network of Fig. 1a, the regulators of gene PRR9 are two other genes, TOC1 and LHY. So if i = P R R9, then π i = {T OC1, L HY }. A regulator can either act as activator or as repressor, and we represent that with the binary variable I u,i , with I u,i = 1 indicating that gene u is an activator for gene i, and I u,i = 0 indicating that gene u is an inhibitor for gene i. For the example above, LHY is an activator for PRR9, hence I u,i = 1, while TOC1 is an inhibitor for PRR9, hence I u,i = 0. From the fundamental equation of transcriptional regulation based on Michaelis-Menten kinetics we have for the gradient of x i (Barenco et al. 2006): where the sum is over all genes u that are in the regulator set of π i of gene i. The first term, −v 0,i x i (t ), takes the degradation of x i into account, while v u,i and k u,i are the maximum reaction rate and Michaelis-Menten parameters for the regulatory effect of gene u ∈ π i on gene i, respectively. See the supplementary material of Pokhilko et al. (2010Pokhilko et al. ( , 2012 for similar examples in the mathematical biology literature. Without loss of generality, we now assume that π i is given by π i = {x 1 , . . . , x s }. Equation (35) can then be written in vector notation: where is the vector of the maximum reaction rate parameters, and the vector D i,t depends on the measured concentrations x u (t ) and the Michaelis-Menten parameters k u,i (u ∈ π i ) via Eq. (35): We combine the s Michaelis-Menten parameters k u,i in a vector K i = (k 1,i . . . , k s,i ) . For n time points t ∈ {t 1 , . . . , t n } we obtain n row vectors from Eq. (37), and we can arrange them successively in an n-by-(|π i | + 1) design matrix The corresponding gradient vector is given by y i := (y i,1 , . . . , y i,n ) , where y i, j is the gradient of x i at time point t j . With y i being the response vector the likelihood is: where ν > 0 is a hyperparameter, and the subscript, {K i ≥ 0}, indicates the truncation condition, i.e. that each element of K i has to be non-negative. For the maximum reaction rates, we use a truncated ridge regression prior: where δ 2 i is a hyperparameter that regulates the prior strength. For σ 2 i and δ 2 i we use inverse Gamma priors, σ 2 i ∼ I G(a σ , b σ ) and δ 2 i ∼ I G(a δ , b δ ). A graphical model representation can be found in Fig. 2.
The posterior distribution of the parameters and hyperparameters has no closedform solution, and we therefore resort to an MCMC scheme to sample from it. From the graphical model in Fig. 2 it can be seen that with the sole exception of the Michaelis-Menten parameters K i , the conditional distribution of each parameter conditional on its Markov blanket 7 is of standard form (due to conjugacy) and can be sampled from directly. The MCMC scheme is therefore of the form of a Gibbs sampler, in which all parameters are sampled directly from their conditional distributions, except for K i , which is sampled via a Metropolis-Hastings within Gibbs step. The conditional distribution of the maximum rate parameter vector V i is obtained from Eqs. (25-26) by replacing θ by V i , and adding an index i, for association with gene i, to all other quantities except for the identity matrix I and the inverse temperature τ . The derivation of the other conditional distributions is straightforward. Pseudo code of the standard MCMC algorithm can be found in Aderhold et al. (2017). Pseudo code of the modified MCMC algorithm integrated into the proposed NETI-DIFF scheme is provided in the Appendix, Table 2.
The data used for inference were obtained from Aderhold et al. (2014). These are synthetic gene expression time series, which were generated from a biologically realistic simulation of the molecular interactions in these networks, using the mathematical framework described in Guerriero et al. (2012) and implemented in the Biopepa software package (Ciocchetta and Hillston 2009). These time series correspond to gene expression measurements in 2h intervals over 24 h, repeated 11 times for different experimental conditions related to various gene knockouts. We repeated the simulations twice, for both of the two networks shown in Fig. 1. Hence, the true interaction network is known, which can be used to evaluate the accuracy of Bayesian model selection based on the modelling framework described above.

Results
In this section, we compare the efficiency and accuracy of three algorithms: standard thermodynamic integration (TI-standard) and optimal thermodynamic integration (TIoptimal) for computing the log marginal likelihood, and the proposed non-equilibrium thermodynamic integration for directly targeting the difference of the log marginal likelihood (NETI-DIFF).
In TI-standard we compute, based on Eq. (5) TI-optimal uses the two improvements proposed in Friel et al. (2014): the log marginal likelihood is computed with the improved numerical integration (Eq. 10), and the inverse temperatures are set iteratively according to an optimality criterion that aims to minimise the expected uncertainty; see Friel et al. (2014) for details. 8 Finally, NETI-DIFF is the algorithm proposed in the present article. For each inverse temperature τ in TI-standard and TI-optimal, we discarded the first 20% of the MCMC steps as burn-in (following Friel et al. (2014)). For NETI-DIFF, we discarded the first 1000 MCMC steps with the inverse temperature kept fixed at τ = 0, as burn-in. 9 We recorded the total number of non-burn-in MCMC steps for all three algorithms, N iter . As discussed in Appendix 7.7 this is a measure of the total computational complexity.
We repeated the MCMC simulations N simu = 5 times from different initialisations. Let B i denote the log Bayes factor obtained from the ith MCMC simulation, and B true the 'true' log Bayes factor. For the Bayesian linear regression models applied to the 8 Note that there is a typo in Eq. (17) of Friel et al. (2014) . 9 Due to the non-equilibrium nature of NETI, not discarding any burn-in phase made little difference to the results. Radiata and Radiocarbon data, a closed-form expression for B true is available. For the Bayesian logistic regression model applied to the Pima Indians data, and the hierarchical Bayesian model from Fig. 2 for biopathway data, the log Bayes factor is not analytically tractable, and B true was obtained from a very long simulation, as in Friel et al. (2014). We assessed the intrinsic estimation uncertainty in terms of the variance: and the accuracy in terms of the mean absolute error:

Radiata pine and Pima Indians
We start our empirical evaluation study with the analysis of the Radiate pine data (Sect. 4.1) and the Pima Indians data (Sect. 4.2). These two data sets have been used in the literature before for the evaluation of the TI method proposed by Friel et al. (2014), and in both cases the goal is to estimate the Bayes factor between two competing Bayesian regression models. For the Radiata pine data we compare two non-nested linear regression models. For the Pima Indians data we compare two logistic regression models, where the first model, M 1 , is nested in the second, M 2 . We apply the NETI-DIFF approach with a sigmoid inverse temperature ladder, defined in Sect. 3.6, and we instantiate NETI-DIFF such that in both applications the transition path runs from the first model, M 1 (τ = 0), to the second, M 2 (τ = 1). For the Pima Indians data this is the natural path, as M 1 is nested within M 2 . Figures 3 and 4 show the average absolute deviations (Eq. 41) between the analytically computed log Bayes factors and the estimated log Bayes factors for different total iteration numbers N iter . Figure 6 compares the variance of the log Bayes factor estimates for the three different methods: TI-standard, TI-optimal, and NETI-DIFF. Figure 7 shows ratios of the variances obtained with TI-optimal and NETI-DIFF.
For the Radiata data, NETI-DIFF only achieves a slight reduction in the absolute deviation (Fig. 3) and the variance (Figs. 6, 7) for the lowest number of iterations, N iter = 64k; otherwise NETI-DIFF and TI-optimal are on a par. Note that the two alternative linear regression models applied to the Radiata data only share the intercept, while their sets of covariables are disjunct. This lack of model overlap presents the least favourable scenario for NETI-DIFF, and our results confirm that there is little room for improvement over standard TI.
For the Pima Indians data, NETI-DIFF achieves a significant reduction in the absolute deviation (Fig. 4) and the variance (Figs. 6, 7) and clearly outperforms both TI methods: TI-standard and TI-optimal. The variance reduction ranges between ratios of 5 and 50. As opposed to the models applied to the Radiata data, the alternative logistic regression models applied to the Pima Indians data are nested, with the parameters of the less complex model forming a subset of those of the more complex one. Our Fig. 3 Average absolute error on the Radiata pine data: TI versus NETI-DIFF. The figure shows the average absolute deviation between the estimated and the true log Bayes factor in dependence on the total number of MCMC iterations N iter . In each panel the same NETI-DIFF results are shown, while the two TI approaches (TI-standard and TI-optimal) were applied with different numbers of discretisation points (10, 20, 50 and 100). The error bars represent standard deviations. The horizontal axes give the total number of (power posterior) MCMC iterations, N iter results demonstrate that in this scenario, the new thermodynamic integration path of NETI-DIFF has potential for significant improvement over the established TI methods.
We also investigated the effect of the inverse temperature ladder ('sigmoid' vs. 'power 5') and the starting point (M 1 vs. M 2 ). To this end, we systematically applied the proposed NETI-DIFF approach with all four combinations (two inverse temperature ladders times two starting points) to the two data sets: Radiata pine and Pima Indians. The results can be found in Fig. 5. First, consider the Pima Indians data, where the two alternative models are nested, and the power inverse temperature scheme of Eq. (8) has been applied. There is a clear advantage of starting the thermodynamic integration at the less complex model over starting at the more complex model: the absolute errors are significantly higher in the latter case. This is not surprising. It is well known from standard TI for computing marginal likelihoods that for the power law of Eq. (8), the optimal transition path is from the prior to the posterior, with the majority of the inverse temperature points at the prior end. Applying this principle to NETI-DIFF, starting the transition path for the differential parameters (i.e. the parameters that are only in the more complex model) at the prior, implies that the overall inverse temperature transition path has to lead from the less complex to the more complex model, in confirmation of our findings. Interestingly, for the sigmoid temperature ladder from Sect. 3.6, the difference between the two directions is substantially reduced, which is a natural consequence of the symmetry inherent in this scheme. There is no significant performance difference between the sigmoid and the power law inverse temperature paths when the models are nested (Pima Indians data, top row in Fig. 5). For the Radiata pine data on the other hand (bottom row in Fig. 5), where the alternative models are not nested, the power law of Eq. (8) is intrinsically suboptimal, 10 and the sigmoidal inverse temperature path of Sect. 3.6 is to be preferred.

Radiocarbon dating
Next, we consider model selection amongst different polynomial orders for polynomial regression on the Radiocarbon data. Since this is a linear model, the log Bayes factor is known and can be used for evaluating the accuracy of the different thermodynamic integration schemes. Besides comparing the proposed NETI-DIFF scheme with the established TI methods, we investigate the influence of the inverse temperature ladder and the transition path. Due to the comparatively low computational costs, we have increased the number of discretisation points from K ∈ {10, 20, 50, 100} to K ∈ {20, 50, 100, 200}. Figure 8 shows the absolute error (see Eq. 41) for NETI-DIFF and the better of the two established TI methods: TI-optimal. The task is to compute the log Bayes factor for the pairwise comparison of various polynomial orders, as indicated by the horizontal axis of each panel. It turns out that for TI-optimal, the accuracy of the estimate deteriorates with increasing difference of the model orders (black bars in the top panels of Fig. 8), while NETI-DIFF is unaffected by model choice. 11 In addition, NETI-DIFF considerably outperforms TI-optimal for the lower iteration numbers, as again seen from the top row in Fig. 8.
The right column of Fig. 6 compares the variances between NETI-DIFF and TIoptimal, and the right column of Fig. 7 shows the corresponding variance ratios. It is seen that NETI-DIFF consistently outperforms TI-optimal, with the variance ratios ranging between 5 and 2000. It appears that for low iteration numbers N iter , the improvement is most pronounced when the alternative models differ substantially (polynomial order 1 vs. 9), while for high iteration numbers N iter , the clearest improvement is achieved when the alternative models are more similar (polynomial orders 4 vs. 6).
The left panel of Fig. 9 compares the two inverse temperature ladders: the power law of Eq. (8) versus the sigmoidal form of Sect. 3.6. Since the models are nested, we would expect the polynomial scheme to perform well, like for the Pima Indians data discussed above. Interestingly, the sigmoidal scheme achieves a better stabilization of the results w.r.t. model order, and a slightly better performance for the largest difference between the polynomial orders of the two alternative models considered. To shed more light on this trend, we have investigated the evolution of the standard deviation of the thermodynamic integral up to a given inverse temperature τ . The results are shown in Fig. 10. While the power law indeed achieves a lower standard deviation than the sigmoidal scheme at the low-inverse-temperature end (near the low-complex model), it contributes a larger proportion to the standard deviation at the high-inverse-  40), for the TI-standard, TIoptimal and NETI-DIFF estimators of the log Bayes factor. For the Radiata pine and the Pima Indians data we varied the number of total MCMC iteration (horizontal axes). For the Radiocarbon data we performed N iter = 1024k iterations and considered four different pairwise model comparisons (horizontal axis). The rows represent different numbers of discretisation points for TI (NETI-DIFF is unaffected). The three columns refer to the four panels in Fig. 3 (right), Fig. 4 (center) and Fig. 8 (right). The corresponding ratios of the variances are shown in Fig. 7 temperature end (near the high-complex model). This suggests that the sparsity of inverse temperatures at the high-inverse-temperature end can be counterproductive due to insufficient sample size.
We finally investigated different model transition paths, with a comparison of three alternative schemes: (1) a staggered path from the low-complexity to the highcomplexity model via a series of all intermediate models; (2) Fig. 6). The horizontal reference line at value 1 indicates equal performance; values above 1 indicate that NETI-DIFF achieves a variance reduction over the established TI schemes. For the Radiata pine and the Pima Indians data we varied the number of total MCMC iterations N iter (horizontal axes). For the Radiocarbon data we performed N iter = 1024k iterations and considered four different pairwise model comparisons (horizontal axis). The rows represent different numbers of discretisation points for TI (NETI-DIFF is unaffected). The three columns refer to the four panels in Fig. 3 (left), Fig. 4 (center) and Fig. 8 (right) shown in the right panel of Fig. 9. The differences are small without a clear trend. This suggests that NETI-DIFF is remarkably robust w.r.t. the choice of the model transition path.

Biopathway
For the biopathway example, we considered two types of data. The first type was obtained from the wild type gene regulatory network shown in Fig. 1a; the second type was obtained from the mutant network shown in Fig. 1b. As we do not have a closed-form expression of the log Bayes factor we chose, as a proxy, the average of the log Bayes factors obtained with the longest TI and NETI-DIFF simulations, which tended to be in reasonably good agreement. Table 1 shows the values of the log Bayes factor thus obtained, which confirms that Bayesian model selection based on the hierarchical model of Fig. 2 consistently identifies the true gene network.
In a preliminary study, we compared the two inverse temperature ladders for NETI-DIFF: power law (see Eq. (8)) with power 5, as in Friel et al. (2014), versus the sigmoid transfer function of Sect. 3.6. We repeated the simulations on the 5 data sets of Table 1. From these data sets, we computed the mean of the variance V, Eq. (40), and the mean absolute error A, Eq. (41). The results are shown in Fig. 11. The trend is not as clear as in Fig. 5. However, the sigmoid inverse temperature ladder achieves more often a  performance improvement over the power law (in terms of lower mean absolute error A and average variance V) than the other way round, and we therefore adopted it for all subsequent studies. The main question of interest is to compare TI and NETI-DIFF with respect to accuracy, estimation uncertainty and computational efficiency. To improve the clarity of the presentation, we only show the comparison between NETI-DIFF and TI-optimal, i.e. the TI scheme with the improvements proposed by Friel et al. (2014). In what follows, we refer to "TI-optimal" simply as "TI". The simulations were repeated for different total iteration lengths, N iter , ranging from N iter = 10,000 to N iter = 6,400,000 MCMC steps. We repeated TI for different numbers of inverse temperatures, K , ranging from K = 10 to K = 100 [(the same values as used in Friel et al. (2014)]. Figure 12 shows the distribution of estimated log Bayes factors obtained from N simu = 5 independent MCMC runs. 12 The two columns refer to the different data types (from the wild type network, left column, and the mutant network, right column), and the rows (Panels 12a-d) to the number of inverse temperatures used for TI (from K = 10 to K = 100; note that NETI-DIFF is unaffected by that choice). The horizontal dashed lines show the 'true' value, as described above. As expected, the distribution width tends to decrease with increasing computational costs, N iter , and for the highest value, TI and NETI-DIFF tend to be in close agreement, with distributions tightly focused on the 'true' values. However, for lower computational costs, N iter ≤ 400k, bias and uncertainty tend to be considerably lower for NETI-DIFF than for TI, irrespective of the number of inverse temperatures used for TI.
For a more systematic investigation, we repeated the MCMC simulations on ten independent data instantiations, for the ten data sets used in Table 1. Five data sets were obtained from the biopathway of Fig. 1a (wildtype), and five data sets were obtained from the biopathway of Fig. 1b (PRR7/PRR9 mutant). For each data set, we computed the mean absolute deviation A, defined in Eq. (41), and the variance V, as defined in Eq. (40).
The top row in Fig. 13 shows the average variance V, averaged over all data instantiations. The second row shows the ratio of the average variance obtained with TI, divided by the average variance obtained with NETI-DIFF, averaged over all five data instantiations: V(TI)/V(NETI − DIFF). The third and fourth rows show the distribution of the variance ratios V(TI)/V(NETI − DIFF) over the five different data instantiations, for different numbers of inverse temperatures (for TI), and different total interation numbers N iter . For all ratios, values above 1 indicate a performance improvement with NETI-DIFF over TI. Our results indicate that NETI-DIFF consistently achieves a considerable variance reduction over TI. This reduction is particularly pronounced for small numbers of inverse temperatures, where it reaches up to three orders of magnitude. However, even for the highest number of inverse temperatures the variance reduction NETI-DIFF achieves over TI still varies between one and two orders of magnitude. This clear reduction in estimation uncertainty is matched by a consistent reduction in the estimation error, as quantified in terms of A and shown in Fig. 14. The reduction becomes stronger with decreasing iteration numbers N iter and decreasing numbers of inverse temperatures, which indicates that the performance improvement of NETI-DIFF over TI is particularly relevant in the regime of limited computational resources.  Fig. 1a (wildtype), and M 2 is the biopathway from Fig. 1b (PRR7/PRR9 mutant). NETI-DIFF is the same for all four rows. Left column: data generated from M 1 ; negative log Bayes factors select the correct model. Right column: data generated from M 2 ; positive log Bayes factors select the correct model. The horizontal line shows the 'true' value of the log Bayes factor (in the sense defined in the text). The box plots show distributions over 5 independent MCMC runs. The horizontal axis shows N iter , the total number of iterations, ranging from 10 to 6400 k a TI with K = 10 inverse temperatures. b TI with K = 20 inverse temperatures. c TI with K = 50 inverse temperatures. d TI with K = 100 inverse temperatures

Discussion
The objective of our work has been the direct targeting of the log Bayes factor via a modified thermodynamic integration path. This has been motivated by statistical physics, where the computation of a reaction free energy (mathematically equivalent to the log Bayes factor) is computationally more efficient than the computation of the difference of standard free energies (equivalent to the difference of log marginal likelihoods). The modified transition path directly connects the posterior distributions of the two models involved. In this way, the high variance prior regime is avoided. We have carried out a comparative evaluation with the state-of-the-art TI method of Friel et al. (2014). Our study confirms that a substantial variance reduction can be achieved when the models to be compared are nested. There is little room for improvement when comparing non-nested models with non-overlapping parameter sets. However, even in this least favourable case, the performance achieved with the proposed method, referred to as NET-DIFF in the present manuscript, is still on a par with established TI methods. For inference in a complex systems described by coupled nonlinear differential equations (biopathway), we found that NETI-DIFF reduces the variance by up to two orders of magnitude over state-of-the-art TI methods. Our work has also revealed that NETI-DIFF achieves a considerable performance stabilisation with respect to a variation of the parameter prior.
When the task is model selection out of a set of cardinality m, carrying out direct pairwise comparisons is of computational complexity m 2 and may not be viable in practice. However, rather than reverting to the standard TI scheme and computing the marginal likelihoods p(D|M 1 ), . . . , p(D|M m ) it appears more sensible to compute the Bayes factors where M 0 is a typical or representative model chosen from the set of models compared. The results for the Radiocarbon data, reported in Sect. 5.2, have demonstrated a remarkable robustness of the proposed method w.r.t. a variation of the model transition path, meaning that there is no significant difference in efficiency and accuracy between the direct computation of log p(D|M 1 ) p(D|M 2 ) , and the indirect computation via log p(D|M 1 ) and log p(D|M 2 ) p(D|M 0 ) . This suggests that 1-out-of-m model selection can also be improved with the method we have proposed. It is beyond the scope of this article to investigate this conjecture at greater depth, but it appears plausible that targeting Bayes factors along an annealing path starting from a reference posterior distribution associated with a reference model should give smaller posterior variance than conventionally targeting marginal likelihoods along an annealing path starting from the prior distribution. If where the two Bayes factors on the right are known from Eq. (43), we get: (44) is formally equivalent to Eq. (4) in Berger and Delampady (1987). We have m models with discrete prior probabilities π i = P(M i ) > 0 and m i=1 π i = 1. We get, e.g., for model M 1 : where B is the Bayes factor: and the hypotheses stand for: H 0 :M = M 1 and H 1 :M ∈ {M 2 , . . . , M m } which are assumed to be true with the prior probabilities π 1 and 1 − π 1 , respectively. Equation (45) corresponds to Eq. (2) in Berger and Delampady (1987). 13 One of the referees raised the interesting question of how the proposed method is applied to graphical Gaussian models and mixture models.
We have included an additional section in the Appendix 7.4 where we discuss in detail how the proposed method can be applied to Graphical Gaussian models. We have also carried out an additional simulation study to illustrate the application of our method to Graphical Gaussian models. The key idea is to not apply the method to the configuration space of precision matrices directly, which would be cumbersome due to the constrained topology of this space (restriction to positive definite matrices). Instead, we make use of the theorem that every multivariate normal density can be represented by a Gaussian belief network, and vice versa; see Geiger and Heckerman (1994). This effectively defines an isomorphism between the space of Gaussian graphical models and the space of Gaussian belief networks. We exploit this isomorphism by defining the proposed NETI scheme in the space of Gaussian belief networks, as discussed in detail in Appendix 7.4.
For mixture models, the proposed NETI method will not achieve any improvement over the standard thermodynamic integration scheme. The reason is that according to Eq. (18), the modified thermodynamic integration path that we have proposed has the potential for a variance reduction if the two model likelihoods in the numerator and denominator share a substantial number of parameters. For mixture models, this is not the case, due to the intrinsic identifiability problem. In Appendix 7.5, we demonstrate on an empirical simulation study that for a mixture model, the proposed new method and the established thermodynamic integration scheme are on a par.
The focus of our study has been a comparison with the improved TI method proposed in Friel et al. (2014). Recently, a powerful new method for variance reduction in thermodynamic integration based on control variates, termed CTI (controlled thermodynamic integral), has been proposed (Oates et al. 2016). The idea is to add a zero-mean function from a given function family (e.g. a polynomial) to the integrand and then apply variational calculus to minimise the variance of the estimator. The resulting optimality equations depend on expectation values w.r.t. the unknown posterior distribution, which the authors approximate with samples from initial MCMC simulations.
On the Radiata data, CTI outperforms NETI-DIFF, due to the fact that NETI-DIFF offers little room for improvement on non-nested models with disjunct parameter sets, as discussed above. On the Pima Indians data, both NETI-DIFF and CTI achieve a significant variance reduction over the state-of-the-art TI method of Friel et al. (2014). 13 Berger and Delampady (1987) study the Bayesian test problem: H 0 : θ = θ 0 vs. H 1 : θ = θ 0 , where θ is a continuous parameter. In Berger and Delampady (1987)  = m j=2 p(D|M j ) · g(M j ). Oates et al. (2016) applied their method with the standard trapezoid sum of Eq. (7), CTI-1, and with the improved trapezoid sum of Eq. (10), CTI-2. A comparison between Fig. 7 in the present paper and Fig. 3 in Oates et al. (2016) shows that the performance of NETI-DIFF, which reduces the variance over state-of-the-art TI by a whole order of magnitude, lies between CTI-1 and CTI-2. Oates et al. (2016) argue that the linear curvature sum of Eq. (7) is known to be biased, and the quadratic curvature rule of Eq. (10) should be used. However, in Aderhold et al. (2017) it was demonstrated that quadratic curvature can lead to an increase in the estimation error when vague prior distributions are used, and it is therefore not always the automatic method of choice. Current work in statistics is increasingly aiming to tackle more complex models, e.g. based on coupled nonlinear differential equations, like the biopathway model discussed in Sect. 4.4. For data generated from an ordinary differential equation model of circadian regulation (Goodwin oscillator), Oates et al. (2016) found that CTI achieved little improvement over state-of-the-art TI. The authors discuss that a potential problem CTI faces for complex models is multimodality of the posterior distributions, rendering the approximation of the posterior expectation values, which enter the optimality equations from variational calculus, less reliable. NETI-DIFF, on the other hand, does not rely on such estimates. In fact, our results, presented in Fig. 13, suggest that NETI-DIFF achieves the most substantial variance reduction over state-of-the-art TI for the most complex, nonlinear biopathway model, reaching up to and exceeding two orders of magnitude.
We conclude that CTI and NETI-DIFF are not competing methods, but rather conceptionally different approaches with the potential to complement each other. CTI aims to achieve variance reduction by adding control variates to the integrand; it requires a reliable estimation of posterior averages of quantities related to these control variates from initial MCMC runs. NETI-DIFF aims to achieve variance reduction by modifying the thermodynamic integration path; it works best for models with substantial parameter overlap. Both approaches can be combined, that is, the natural next step is to add control variates and change the integration path, i.e. to target the log Bayes factor with the principles of CTI. This combination of NETI-DIFF and CTI has the potential to further extend the feasibility of Bayesian model selection to ever more complex models, and a closer investigation of such a hybrid approach poses a promising avenue for future research.

Proof of Jarzynski's theorem
Using the definitions from Sect. 2, we get:

Uncertainty quantification
From Eq. (16) we have: For the first term we get: Combining Eqs. (46) and (47), we get: Define the following shorthand notation: which is an estimator of E τ log p(D|θ,M 2 ) p(D|θ,M 1 ) with sample size 1. We can rewrite Eq. (18) as: For the variance we get: Table 2 shows the NETI-DIFF pseudocode for the Bayesian hierarchical model of Fig. 2. Pseudocode for standard MCMC, following a Metropolis-Hastings within Gibbs scheme, was provided in Table 1 of Aderhold et al. (2017). Table 2 shows the modification required to sample with the NETI-DIFF scheme from the tempered posterior distribution in Eq. (13).

Application to Gaussian Graphical Models
In this Appendix we show how the new method (NETI-DIFF) can be used to infer the Bayes factor between Gaussian graphical models (GGMs). We propose an indirect procedure which exploits that multivariate Gaussians can be represented as 'Gaussian belief networks' (Geiger and Heckerman 1994). A Gaussian graphical model corresponds to an M-dimensional multivariate Gaussian distribution with mean vector m and covariance matrix Σ so that the density (PDF) is given by where x = (x 1 , . . . , x M ) T and W = Σ −1 is called the precision matrix. Each 0 element of W indicates that the partial correlation between the corresponding variables is zero, e.g. W i, j = 0 if the partial correlation between x i and x j is zero. We follow Geiger and Heckerman (1994) and identify this Gaussian distribution with a 'Gaussian belief network', i.e. we factorise the density in Eq. (52) with the chain rule: where the conditional distributions are univariate Gaussians The differences to the standard MCMC scheme from Table 1 in Aderhold et al. (2017) The most convenient way to compute the Bayes factor between two competing GGMs is to work with their Gaussian belief network representations. 14 For a GGM with precision matrix W, we impose a Wishart prior onto W, and we represent the GGM in terms of the parameters m = (m 1 , . . . , m M ) T , σ 2 = (σ 2 1 , . . . , σ 2 M ) T , and Given two GGMs M 1 and M 2 with precision matrices W 1 and W 2 we represent both as Gaussian belief networks with the regression coefficient matrices B k whose elements are given by β k j,i (k = 1, 2). We have β k j,i = 0 if W k j,i = 0 (k = 1, 2) and β 1 j,i = β 2 j,i if β 1 j,i , β 2 j,i = 0, so that all shared non-zero regression coefficients are equal. We assume that both GGMs share the mean vector m, which we assume to be known, and the conditional variances (σ 2 i ) k = σ 2 i . Let B denote the matrix of all regression coefficients which are non-zero in at last one of the GGMs. The elements of B are: Given n data points x 1 , . . . , x n the tempered posteriors take the form: where τ ∈ [0, 1] and the three precision matrices W 1 , W 2 , and W can be computed with Eq. (55) from the conditional variances σ 2 i and the regression parameters in B 1 , B 2 and B. Sampling from the tempered posterior can be done with Metropolis-Hastings (MH) MCMC moves which we define in the space of the non-zero regression parameters in B and in the space of the logarithms of the conditional variances σ 2 i . We obtain a new candidate state B and σ 2 i, by adding randomly sampled numbers to the non-zero elements of B and to log(σ 2 i ). 15 From the new candidate matrix B we extract the matrices B k, (k = 1, 2) as follows: β k, j,i = 0 if W k i, j is restricted to be zero and β k, j,i = β j,i otherwise. The new precision matrices W , W 1, , and W 2, can then be computed from B , B 1, , and B 2, with Eq. (55), and the MH acceptance probability depends on the ratio of the tempered posteriors of the new precision matrix W and the old precision matrix W.
For a proof of concept we perform a simulation study: We consider the M = 7 genes (1=LHY, 2=TOC1, 3=PRR9, 4=PRR7, 5=GI, 6=Y, and 7=TOC1) of the Arabidopsis networks, shown in Fig. 1, and we parametrize both graphs M 1 and M 2 as Gaussian belief networks. We set: m = 0 and σ 2 i = 1 for all i, and the non-zero regression coefficients appearing in both graphs are set to β j,i = 1, while the regression coefficients appearing only in the wildtype (M 1 ) are set to β ∈ R. The latter coefficients correspond to the edges 'PRR9-PRR7' (β 2,3 ) and 'PRR7-NI' (β 3,4 ) in Fig. 1. β is a tuning parameter for the strength of the two additional partial correlations in M 1 . For β = 0 the partial correlations are zero and the nested mutant network M 2 is the correct model. As prior on W we use a Wishart distribution with df = 10 degrees of freedom and the identity matrix as precision matrix P = I 7 . We generate data sets with n = 100 data points from M 1 , and we use NETI-DIFF (with 100k iterations, a sigmoidal temperature ladder and ε = 0.1) to compute the Bayes-factors. Figure 15 shows the results. The Bayes factors are in favour of the true wildtype network (M 1 ) if the additional regression coefficients have a sufficient size (β ≥ 0.3). For low values (β ≤ 0.2) the Bayes factor is in favour of the mutant network (M 2 ), which is actually the true network for β = 0. Only for low positive values (β = 0.1 and β = 0.2) the wrong model is favoured over the true model. The latter can be explained by the prior. The Wishart prior with hyperparameters df = 10 and P = I 7 corresponds to 10 pseudo data points from a GGM without any non-zero partial correlations, and hence, yield a higher penalty for the wildtype netwerk M 1 than for the sparser mutant network M 2 .

Application to mixture models (galaxy data)
In this Appendix we show that the new method (NETI-DIFF) can also be used to compute the Bayes factor between mixture models with different numbers of mixture components. Like Friel et al. (2014) we consider the Galaxy data from Richardson and Green (1997), which contain n = 82 measurements y 1 , . . . , y 82 of galaxy velocities, and we compute the Bayes factor between two Bayesian Gaussian mixture models M 3 with K = 3 components and M 4 with K = 4 components. For our study we use exactly the same mixture model as Friel et al. (2014) with the same prior distributions and the same hyperparameters. A latent allocation vector z = (z 1 , . . . , z 82 ) T allocates the individual data points to the K mixture components, where z i = k if data point y i has been allocated to component k (k = 1, . . . , K ; i = 1, . . . , 82). On the mixture weights w k := P(z i = k) we impose a Dirichlet prior: The data points within each component k (1 ≤ k ≤ K ) are assumed to stem from a univariate Gaussian distribution with mean μ k and variance σ 2 k , so that and for μ k and σ 2 k we use a Gaussian prior and an Inverse-Gamma prior: We define θ K to be the set of all parameters of the mixture model M K with K components: In the absence of limiting conditions, mixture models with different numbers of components (here: M 3 and M 4 ) do not share any parameters, and the tempered NETI-DIFF posteriors take the form p τ (θ 3 , θ 4 |y 1 , . . . , y m , M 3 , M 4 ) ∝ p(y 1 , . . . , y n |θ 3 ) τ p(y 1 , . . . , y n |θ 4 ) 1−τ Because of this modular form, the parameters in the sets θ 3 and θ 4 can be sampled by disjunct MCMC sampling steps, which either re-sample subsets of the parameters θ 3 ⊂ θ 3 (orθ 4 ⊂ θ 4 ) from their full conditional distributions: p τ θ 3 |θ 3 , θ 4 , y 1 , . . . , y m , M 3 , M 4 ∝ p (y 1 , . . . , y n |θ 3 ) τ · p (θ 3 |M 3 ) p τ θ 4 |θ 4 , θ 3 , y 1 , . . . , y m , M 3 , M 4 ∝ p (y 1 , . . . , y n |θ 4 ) 1−τ · p (θ 4 |M 4 ) whereθ K ∪θ K = θ K , or via Metropolis Hastings sampling steps, whose acceptance probabilities are: where HR is the move-specific Hastings ratio and the symbol indicates a new candidate parameter set. Since these are the standard equations for power posterior sampling, as used by the thermodynamic integration (TI) approach, the adaptation of the Metropolis-Hastings and Gibbs sampling steps of the power posterior sampling scheme for TI (Friel et al. (2014)) is straightforward. At each temperature τ ∈ [0, 1] NETI-DIFF updates the parameters in θ 3 and in θ 4 independently by performing the corresponding steps of the MCMC sampling scheme. The only difference is that the parameters in θ 4 are subject to the complementary temperature 1 − τ rather than τ , and we therefore implement NETI-DIFF with the sigmoid inverse temperature ladder from Sect. 3.6. Moreover, we also take into account that NETI-DIFF has to perform twice as many sampling steps as TI, since NETI-DIFF re-samples the parameters of both models M 3 and M 4 within each iteration. Thus, NETI-DIFF iterations are approximately double as expensive as TI iterations, and we can perform only 50% of the total number of iterations N iter with NETI-DIFF. In our empirical study we compare the performance of NETI-DIFF with TI-standard and TI-optimal, and we implement both TI approaches with 100 discretisation points. We compute the Bayes factor between the mixture models M 3 and M 4 based on N iter = 1000k and N iter = 2000k iterations. 16 The results of our study are shown in Fig. 16. It can be seen that there are no significant differences between the performances. The NETI-DIFF estimates appear to be minimally less biased than the TI estimates, but on the other hand the NETI-DIFF estimates have a slightly increased standard deviation. This finding, that NETI-DIFF does not lead to any improvement over the standard TI approach, is not surprising: Due to the fact that the two mixture models do not have any parameters in common, targeting the Bayes factor directly cannot have any advantages. For models with disjunct parameter spaces NETI-DIFF effectively just corresponds to two simultaneously performed but independent non-equilibrium thermodynamic integration (NETI) approaches, where one model is subject to the complementary temperature transition from τ = 1 to τ = 0. Targeting the Bayes factor directly, as described in Sect. 3.3, can only lead to an improvement if the two models share parameters. In the direct transition paths between the two model posteriors, only those shared parameters constantly appear with the inverse temperature 1 and do not undergo any temperature transitions (i.e. they are excluded from the annealing process). All non-shared parameters have to undergo the transitions from τ = 0 to τ = 1 or from τ = 1 to τ = 0, respectively.

Full conditional distributions of variance parameters
For linear models where the variance parameter σ 2 in Eq. (20) in Sect. 3.5 is not known, a prior distribution has to be imposed on σ 2 . A common choice is the conjugate Inverse-Gamma distribution with hyperparameters a/2 and b/2, symbolically σ −2 ∼ GAM( a 2 , b 2 ). The tempered full conditional distribution of σ −2 is then of closed-form and can be derived as follows: p τ (σ −2 |D, θ, M 1 , M 2 ) ∝ p(y|θ , σ 2 , M 2 ) τ · p(y|θ , σ 2 , M 1 ) 1−τ · p(θ|σ 2 , M 1 , M 2 ) · p(σ 2 ) ∝ N n (D (2) θ, σ 2 I) τ · N n (D (1) θ, σ 2 I) 1−τ · N p (μ 0 , σ 2 δ 2 I) · G AM(σ −2 ) where p is the length of the regression coefficient vector θ and The data size is given by n observations. The number of variables can differ depending on the model. Time values are averages over 10 runs and different models on an Intel Core i7 6700 HQ processor. The standard deviations are indicated in brackets. The runtime for Biopepa is for a single response out of seven possible responsesã Comparing this with the identity: we get the full conditional distribution Hence, we can also sample σ −2 directly from the tempered full conditional distribution in a Gibbs sampling scheme, and σ 2 = 1/σ −2 .

Computational run times and convergence diagnostics
It is important to assess the convergence of the NETI simulations accurately. However, conventional convergence diagnostics for MCMC, like the Gelman-Rubin potential scale reduction factor, are not applicable here. The reason is that the combination of the NETI scheme, described in Sect. 3.2, and the new thermodynamic integration path, described in Sect. 3.3, continuously transform one model into another via a series of non-equilibrium configurations. We need to point out that any samples taken during this transformation are of no interest in themselves; the only quantity of interest is the log Bayes factor, computed according to Eq. (18). The estimate of the log Bayes factor from Eq. (18) is a random variable that is subject to the intrinsic stochasticity of the MCMC sampler. A natural convergence diagnostic is the variance of this estimator: for an infinite simulation time, the variance should go to zero as the estimate should not (a) (b) Fig. 17 Convergence of NETI-DIFF for the Biopepa data. This figure gives the variance of marginal log likelihood estimators for NETI-DIFF and the two ladder types (panel a and b). The variance is calculated from five repetitions and for five different data instances shown in the legend. This figure complements Fig. 11, where the average variances, averaged over all five data instantiations, are shown. a NETI-DIFF with power law ladder. b NETI-DIFF with sigmoidal ladder depend on the particular idiosyncrasies of any MCMC trajectory. We have investigated this conjecture in Figs. 6, 7, 11c, d and 13a. Since Figs. 11c, d and 13a provide average variances over five independent data instantiations, we have included Fig. 17 that shows the individual variances for each data set separately. All these figures demonstrate that the variance approaches zero as the simulation time, regarding the number of MCMC steps, is increased. Figure 6 quantifies the improvement in convergence that the proposed method achieves over the established schemes, in the form of a faster decrease of the variance with increasing simulation times. The figures mentioned above, e.g. Figs. 6 and 13a, monitor convergence in terms of iteration numbers. For a fair comparison between different methods, we also need to take into consideration the computational costs per iteration shown in Table 3: The computational run times of the three algorithms compared are approximately equal; if there is any difference at all, it appears to be in favour of the proposed NETI scheme. From this, we can conclude that monitoring inference uncertainty as a function of MCMC iteration numbers, as carried out throughout our paper, provides an appropriate quantification of computational complexity.