On the marginal likelihood and cross-validation

In Bayesian statistics, the marginal likelihood, also known as the evidence, is used to evaluate model fit as it quantifies the joint probability of the data under the prior. In contrast, non-Bayesian models are typically compared using cross-validation on held-out data, either through $k$-fold partitioning or leave-$p$-out subsampling. We show that the marginal likelihood is formally equivalent to exhaustive leave-$p$-out cross-validation averaged over all values of $p$ and all held-out test sets when using the log posterior predictive probability as the scoring rule. Moreover, the log posterior predictive is the only coherent scoring rule under data exchangeability. This offers new insight into the marginal likelihood and cross-validation and highlights the potential sensitivity of the marginal likelihood to the setting of the prior. We suggest an alternative approach using aggregate cross-validation following a preparatory training phase. Our work has connections to prequential analysis and intrinsic Bayes factors but is motivated through a different course.


Introduction
Probabilistic model evaluation and selection is an important task in statistics and machine learning, particularly when multiple models are under initial consideration. In the non-Bayesian literature, models are typically compared using out-of-sample performance criteria such as cross-validation [Geisser and Eddy, 1979, Shao, 1993, Vehtari and Lampinen, 2002, or predictive information [Watanabe, 2010]. Computing the leave-p-out cross-validation score requires n-choose-p test set evaluations for n data points, which in most cases is computationally unviable and hence approximations such as k-fold cross-validation are often used instead [Geisser, 1975]. A survey is provided by Arlot and Celisse [2010], and a Bayesian perspective on cross-validation by Vehtari and Ojanen [2012], Gelman et al. [2014].
In Bayesian statistics, the marginal likelihood or model evidence is the natural measure of model fit. For a model M with likelihood function or sampling distribution {f θ (y) : θ ∈ Θ} parameterised by θ, a prior π(θ), and observations y 1:n ∈ Y n , the marginal likelihood or the prior predictive is defined as p M (y 1:n ) = f θ (y 1:n ) dπ(θ) .
The marginal likelihood can be used to calculate the posterior probability of the model given the data, p(M | y 1:n ) ∝ p M (y 1:n ) p(M), as it is the probability of the data being generated under the prior when the model is correctly specified [Robert, 2007, Chapter 7]. The ratio of marginal likelihoods between models is known as the Bayes factor that quantifies the prior to posterior odds on observing the data. The marginal likelihood can be difficult to compute if the likelihood is peaked with respect to the prior, although Monte Carlo solutions exist; see Robert and Wraith [2009] for a survey. Under vague priors, the marginal likelihood may also be highly sensitive to the prior dispersion even if the posterior is not; a well known example is Lindley's paradox [Lindley, 1957, O'Hagan and Forster, 2004, Robert, 2013. As a result, its approximations such as the Bayesian information criterion [Schwarz, 1978] or the deviance information criterion [Spiegelhalter et al., 2002] are widely used, see also Gelman et al. [2014]. For our work, it is useful to note from the property of probability distributions that the log marginal likelihood can be written as the sum of log conditionals, where p M (y i | y 1:i−1 ) = f θ (y i ) dπ(θ | y 1:i−1 ) is the posterior predictive for i > 1, p M (y 1 | y 1:0 ) = f θ (y 1 ) dπ(θ) , and this representation is true for any permutation of the data indices. While Bayesian inference formally assumes that the model space captures the truth, in the model misspecified or so called M -open scenario [Bernardo and Smith, 2009, Chapter 6] the log marginal likelihood can be simply interpreted as a predictive sequential, or prequential [Dawid, 1984], scoring rule of the form S(y 1:n ) = i s(y i | y 1:i−1 ) with score function s(y i | y 1:i−1 ) = log p M (y i | y 1:i−1 ). This interpretation of the log marginal likelihood as a predictive score [Kass and Raftery, 1995, Gneiting and Raftery, 2007, Bernardo and Smith, 2009 has resulted in alternative scoring functions for Bayesian model selection [Dawid and Musio, 2014, 2015, Watson and Holmes, 2016, Shao et al., 2018, and provides insight into the relationship between the marginal likelihood and posterior predictive methods [Vehtari and Ojanen, 2012]. Key et al. [1999] considered cross-validation from an M -open perspective and introduced a mixture utility for model selection that trades off fidelity to data with predictive power.

Uniqueness of the marginal likelihood under coherent scoring
To begin, we prove that under an assumption of data exchangeability, the log posterior predictive is the only prequential scoring rule that guarantees coherent model evaluation. The coherence property under exchangeability, where the indices of the data points carry no information, refers to the principle that identical models on seeing the same data should be scored equally irrespective of data ordering.
In demonstrating the uniqueness of the log posterior predictive, it is useful to introduce the notion of a general Bayesian model [Bissiri et al., 2016], which is a framework for Bayesian updating without the requirement of a true model. Define a parameter of interest by where F 0 (y) is the unknown true sampling distribution giving rise to the data, and l : Θ × Y → [0, ∞) is a loss function linking an observation y to the parameter θ. Bissiri et al. [2016] argue that after observing y 1:n , a coherent update of beliefs about θ 0 from a prior π G (θ) to the posterior π G (θ | y 1:n ) exists and must take on the form where l(θ, y 1:n ) = i l(θ, y i ) is an additive loss function and w > 0 is a loss scale parameter; see Holmes and Walker [2017], Lyddon et al. [2019] on the selection of w. For w = 1 and l(θ, y) = − log f θ (y), we obtain traditional Bayesian updating without assuming the model f θ (y) is true for some value of θ. From (3), M -open Bayesian inference is simply targeting the value of θ that minimizes the Kullback-Leibler divergence between dF 0 (y) and f θ (y). The form (4) is uniquely implied by the assumptions in Theorem 1 of Bissiri et al. [2016], and we now focus on the coherence property of the update rule. An update function ψ(l(θ, y), π G (θ)) = π G (θ | y) is coherent if, for some inputs y 1:2 , it satisfies ψ[l(θ, y 2 ), ψ{l(θ, y 1 ), π G (θ)}] = ψ{l(θ, y 1 ) + l(θ, y 2 ), π G (θ)} This coherence condition is natural under an assumption of exchangeability as we expect posterior inferences about θ 0 to be unchanged whether we observe y 1:2 in any order or all at once, as it is in traditional Bayesian updating. We now extend this coherence condition to general Bayesian model choice, where the goal is to evaluate the fit of the observed data under the general Bayesian model class M G = {l(θ, y) : θ ∈ Θ} with a prior π G (θ). We treat w as a parameter outside of the model specification, as there are principled methods to select it from the model, prior and data. We define the log posterior predictive score as s G (ỹ | y 1:n ) = log g{l(θ,ỹ)}dπ G (θ | y 1:n ) where g : [0, ∞) → [0, ∞) is a continuous monotonically decreasing scoring function that transforms l(θ, y) into a predictive score for a test pointỹ. We define the cumulative prequential log score as where s G (y 1 | y 1:0 ) = log g{l(θ, y 1 )}dπ G (θ). The cumulative prequential log score sums the log posterior predictive score of each consecutive data point in a prequential manner, where a large score indicates that the model is predicting well. An intuitive choice for the scoring function might be the negative loss g(l) = −l, but we will see that this violates coherency, as defined below.
Definition 1. The model scoring function g(l) is coherent if it satisfies n i=1 s G (y i | y 1:i−1 ) = log g{l(θ, y 1:n )}dπ G (θ) for all Θ, π(θ) and n > 0, such that S G (y 1:n ) is invariant to the ordering or partitioning of the observations.
We now present our main result on the uniqueness of the choice of g. where w is the loss-scale in the general Bayesian posterior.
Proof. The proof is given in A.1 of the Appendix.
This holds irrespective of whether the model is true or not. More importantly for us is the corollary below.
Corollary 1. The marginal likelihood is the unique coherent marginal score for Bayesian inference.
The marginal likelihood arises naturally as the unique prequential scoring rule under coherent belief updating in the Bayesian framework. The coherence of the marginal likelihood implies an invariance to the permutation of the observations y 1:n under exchangeability, including independent and identically distributed data, a property that is not shared by other prequential scoring rules, such as Dawid and Musio [2014], Grünwald and van Ommen [2017], Shao et al. [2018].

Equivalence of the marginal likelihood and aggregate cross-validation
The leave-p-out cross-validation score is defined as 1:p denotes the t'th of n-choose-p possible held-out test sets, with y (t) 1:n−p the corresponding training set, such that y 1:n = ỹ (t) , y (t) , and S CV records the average predictive score per datum. Although leave-one-out cross-validation is a popular choice, it was shown in Shao [1993] that it is asymptotically inconsistent for a linear model selection problem, and requires {p/n} → 1 as n → ∞ for consistency. We will not go into further detail here but instead refer the reader to Arlot and Celisse [2010]. Selecting a larger p has the interpretation of penalizing complexity [Vehtari and Ojanen, 2012], as complex models will tend to over-fit to a small training set. However, the number of test set evaluations grows rapidly with p and hence k-fold cross-validation is often adopted for computational convenience.
From a Bayesian perspective it is natural to consider the log posterior predictive as the scoring function, s(ỹ | y) = log f θ (ỹ)dπ(θ | y), particularly as we have now shown that it is the only coherent scoring mechanism, which leads us to the following result.
Proposition 2. The Bayesian marginal likelihood is equivalent to the aggregate leave-p-out crossvalidation score using the log posterior predictive as the scoring rule, such that Proof. This follows from the invariance of the marginal likelihood under arbitrary permutation of the sequence y 1:n in (2). We provide more detailed explanations in A.2, A.3 of the Appendix.
The Bayesian marginal likelihood is simply n times the average leave-p-out cross-validation score, n × (1/n) n p=1 S CV (y 1:n ; p), where the scaling by n is due to (6) being a per datum score. Bayesian models are evaluated through out-of-sample predictions on all (2 n − 1) possible heldout test sets whereas cross-validation with fixed p only captures a snapshot of model performance. Evaluating the predictive performance on (2 n − 1) test sets would appear intractable for most applications, but we see through (7) and (1) that it is computable as a single integral.

Sensitivity to the prior and preparatory training
The representation of the marginal likelihood as an aggregate cross-validation score (7) provides insight into the sensitivity to the prior. The last term in the right hand side of (7) involves no training data, S CV (y 1:n ; , which scores the model entirely on how well the analyst is able to specify the prior. In many situations, the analyst may not want this term to contribute to model evaluation. Moreover, there is tension between any desire to specify vague priors to safeguard their influence and the fact that diffuse priors can lead to an arbitrarily large and negative model score for real valued parameters from (7). It may seem inappropriate to penalize a model based on the subjective ability to specify the prior, or to compare models using a score that includes contributions from predictions made using only a handful of training points even with informative priors. For example, we see that 10% of terms contributing to the marginal likelihood come from out-of-sample predictions using, on average, less than 5% of available training data. This is related to the start-up problem in prequential analysis [Dawid, 1992].
A natural and obvious solution is to begin evaluating the model performance after a preparatory phase, for example using 10% or 50% of the data as preparatory training prior to testing. This leads to a Bayesian aggregate leave-P -out cross-validation score defined as S ACV (y 1:n ; P ) = P p=1 S CV (y 1:n ; p) with a matched preparatory cross-validation score S P CV (y 1:n ; P ) = n p=P +1 S CV (y 1:n ; p), for S ACV (y 1:n ; P ) = 1 This equivalence is derived in A.4 of the Appendix in a similar fashion to Proposition 2. This has precisely the form of the the log geometric intrinsic Bayes factor of Berger and Pericchi [1996] but motivated by a different route. The intrinsic Bayes factor was developed in an objective Bayesian setting [Berger and Pericchi, 2001], where improper priors cause indeterminacies in the evaluation of the marginal likelihood. The intrinsic Bayes factor remedies this with a partition of the data into y 1:l , y l+1:n , where y 1:l is the minimum training sample used to convert an improper prior π(θ) into a proper prior π(θ | y 1:l ). In contrast, we set n − P to provide preparatory training and π(θ) can be subjective. Moreover, in modern applications we often have d ≫ n where intrinsic Bayes factors cannot be applied in their original form.
We can approximate (9) through Monte Carlo where the training data sets y (t) 1:n−P are drawn uniformly at random, and for non-conjugate models the inner term must also be estimated, for example throughŜ where samples θ are obtained via Markov chain Monte Carlo. In A.5 of the Appendix we provide some further details on efficient computation of (10).

Illustration for the normal linear model
We illustrate the use of Bayesian aggregate cross-validation in a polynomial regression example, where the r-th polynomial model is defined as We observe the data {y 1:n , x 1:n }, and we place a fixed vague prior on the intercept term, θ 0 ∼ N (θ 0 ; 0, 100 2 ), and θ d ∼ N (θ d ; 0, s 2 ) for d ∈ {1, . . . , r} on the remaining coefficients. In our example, we have n = 100 and the true model is r = 1, θ = 1 0.5 T with known σ 2 = 1. For our prior, we vary the value of s 2 ∈ 10 −1 , 10 0 , 10 4 to investigate the impact of the prior tails. For each prior setting, we calculate log p M (y 1:n ) and S ACV (y 1:n ; P ) for models r ∈ {0, 1, 2}. In this example, log p M (y 1:n ) is tractable, whereas S ACV requires a Monte Carlo average over tractable log posterior predictives. We report the mean over 10 runs of estimating S ACV with T = 10 6 random training/test splits. We calculate the Monte Carlo standard error over the 10 runs and report the maximum for each setting of P . The results are shown in Table 1, where the highest scoring model is in bold for each setting of s 2 , andŜ ACV is normalized to be on the same scale as log p r (y 1:n ). Under the strong prior s 2 = 10 −1 and the moderate prior s 2 = 10 0 , the marginal likelihood correctly identifies the true model, but when we increase s 2 to 10 4 it heavily over-penalizes the more complex models and prefers r = 0. In fact, the magnitude of the marginal likelihood and the discrepancy just described can be made arbitrarily large by simply increasing s 2 , which should be guarded against when a modeller has weak prior beliefs. This issue is not observed withŜ ACV for the values of P we consider. The vague prior does not impede the ability ofŜ ACV to correctly identify the true model r = 1 and the scores are stable within each column of P .
In A.6 of the Appendix, we present graphical tools for exploring the aggregate cross-validation and the effect of the choice of P on S ACV . In A.7 of the Appendix we provide an additional example using probit regression on the Pima Indian data set . Model log p r (y 1:n )Ŝ ACV (y 1:n ; P ) × n/P r P = 0.9n P = 0.5n M = 0.

Discussion
We have shown that for coherence, the unique scoring rule for Bayesian model evaluation in either M -open or M -closed is provided by the log posterior predictive probability, and that the marginal likelihood is equivalent to an aggregate cross-validation score over all training-test data partitions. The coherence flows from the fact that the scoring rule and the Bayesian update both use the same information, namely the likelihood function, which is appropriate as the alternative would be to learn and score under different criteria. If we are interested in an alternative loss function to the log likelihood, we advocate a general Bayesian update [Bissiri et al., 2016, Lyddon et al., 2019] that targets the parameters minimising the expected loss, with models evaluated using the corresponding coherent aggregate cross-validation score.
Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal
As g is continuous and monotonically decreasing, to satisfy (5) it must take on the form for λ ≥ 0. We now explicitly write out the form of p 1 Expanding, cancelling terms, and simplifying we obtain and so we must have λ = 0 or λ = w, where only the latter solution is non-trivial. We have thus shown that for n = 2, |Θ| = 2, the unique non-trivial solution to (5) is The remainder of the proof involves showing that this choice of g satisfies (5) for all n > 0 and all Θ and π(θ). Subbing (A.1.4) into (5), we obtain where for convenience we denote l(θ, y 1:0 ) = 0.

A.2 Proof of Proposition 2
Proof. Consider the (n! × n) matrix Z with elements (Z) ti = log p M (y 1:i−1 ), such that the t'th row of Z records the prequential sequence of log posterior predictives under the t'th of n! permutations of y 1:n . By the property of conditional probabilities, we have that the row sums of Z are equal, i (Z) ti = i (Z) t ′ i for all t, t ′ , and hence Within each column of Z, the values (Z) ti are invariant to the permutation of y 1:i−1 in the preceding i − 1 columns under exchangeability. There are thus n-choose-(i − 1) distinct training sets and n − i + 1 choices for y i given the training set. For each column i ∈ {1, . . . , n}, we can then write We have the result if we let p = n − i + 1.

A.3 Alternative proof of Proposition 2
We first begin by showing the following proposition.
Proposition 3. For a preparatory cross-validation score, S P CV (y 1:n ; P ), defined as the sum of cross-validation terms from leave-(P + 1)-out to leave-n-out, S P CV (y 1:n ; P ) = n p=P +1 S CV (y 1:n ; p), we have the following equivalence relationship S P CV (y 1:n ; P ) = 1 which states that S P CV is the average log marginal likelihood over all choices of the training set.
Proposition 2 then follows trivially by setting P = 0 in Proposition 3.

A.4 Derivation of S ACV for Bayesian models
The following corollary follows easily from Propositions 2 and 3. Proof. We note that log p M (y 1:n ) = S ACV (y 1:n ; P ) + S P CV (y 1:n ; P ) from their definitions and Proposition 2. From the permutation invariance of the marginal likelihood, we can write

A.5 Computing S ACV
To avoid the need for T Markov chain Monte Carlo chains in (10), we can take advantage of the fact that the partial posteriors for different training sets will be similar, and utilize importance sampling [Bhattacharya et al., 2007, Vehtari et al., 2017 or sequential Monte Carlo [Bornn et al., 2010] to estimate the posterior predictives for computational savings. We also note thatŜ ACV in (10) is a biased estimate, and Rischard et al. [2018] provides unbiased estimators of log p M (ỹ 1:P | y 1:n−P ) directly through unbiased Markov chain Monte Carlo and path sampling methods. The arithmetic averaging over training/test splitsŜ ACV may also be inherently unstable, as demonstrated by the following example. Suppose that y is a binary random variable which takes on either 0 or 1 with equal probability, and we are attempting to estimate S ACV (y 1:n ; n/2). For large n, it is likely that approximately half of the values in y 1:n are equal to 0 and the other half to 1. There will thus exist a permutation of the sequence y 1:n such that almost all the first n/2 values are equal to 0, with the remaining almost all equal to 1. The model will then be certain that y = 0 after observing the training set, and score the remaining n/2 points very poorly, giving a large negative log posterior predictive. This suggests that an arithmetic average may be unstable; the median or robust trimmed mean over permutations may be stabler alternatives.
The form in (9) relies on the conditional coherency of Bayesian updating and scoring. Without this, S ACV still exists as defined in (8), and can be directly estimated for example througĥ where p (t) ∼ U{1, P } and the training set y Normalized aggregate cross-validation score SACV Normalized aggregate cross-validation score S ACV against preparatory set size, s 2 = 1 r = 0 r = 1 r = 2 Figure A.6.1: Leave-p-out cross-validation scoreŜ CV (y 1:n ; p) against n − p (left) and normalized aggregate cross-validation scoreŜ ACV (y 1:n ; P ) × n/P against n − P (right) for s 2 = 1 and p, P ∈ {1, . . . , 99}; the maximum standard error is 0.001 forŜ CV and 0.005 forŜ ACV A visualization of the effects of the training/preparatory data size is shown in Figure A.6.1 for s 2 = 1. We omit S CV (y 1:n ; n) and S ACV (y 1:n ; n) for clarity of the plot, as both are significantly more negative than the other values. On the left we see that the individual cross-validation term S CV (y 1:n ; p) prefers the simplest r = 0 model when the training set is very small as over-fitting is penalized, but as n − p increases, the true r = 1 model overtakes it. The r = 2 model eventually overtakes the r = 0 model too, and we see the discrepancy between r = 2 and r = 1 decrease as over-fitting is penalized less and less. This latter effect is demonstrative of how leave-one-out crossvalidation under-penalizes complex models as argued in Shao [1993], and why a value of P > 1 should be preferred. On the right, we observe a similar effect for the aggregate cross-validation score S ACV , but the discrepancy between r = 2 and r = 1 remains more noticeable for moderate n − P as a cumulative sum of S CV terms is being taken.

A.7 Illustration for the probit model
To demonstrate the aggregate cross-validation score in an intractable example, we carry out model selection in the Pima Indian benchmark model with a probit model. We observe binary random variables y 1:n with associated r-dimensional covariates x 1:n , and the probit model is defined as where Φ is the standard normal cumulative distribution function andx = 1 x T T . As suggested in Marin and Robert [2018], we elicit a g-prior π(θ) = N θ; 0 r+1 , g(X T X) −1 where 0 r+1 is a r + 1 vector of 0s and X is the n by r + 1 matrix with rowsx T i . The dataset consists of n = 332 data points and we consider r = 3 covariates consisting of glu, bp and ped, which correspond to plasma glucose concentration from an oral glucose test, diastolic blood pressure and diabetes pedigree function respectively. We compare the full model M 0 : (glu,bp,ped) with M 1 : (glu,bp) through log p M (y 1:n ) and S ACV (y 1:n ; P ) to test for significance of ped. We standardize all covariates to have 0 mean and variance 1. We calculate log p M (y 1:n ) using importance sampling with a Gaussian proposal with 10 3 samples. The proposal mean is set to the maximum likelihood estimate of θ and proposal covariance to the estimated covariance matrix of the maximum likelihood estimate as suggested in Marin and Robert [2018]. For S ACV (y 1:n ; P ), we estimate each posterior predictive in (10) with the same importance sampling scheme where we temper the proposal such that its covariance matrix is divided by (n − P )/n. We also use 10 3 proposal samples and average over T = 10 5 random train/test splits. We carry out 10 runs of each and report the mean and maximum standard error as before.
We see in Table A.7.1 that for g = n, the simpler model with ped omitted performs worse for both scores, and there is thus strong evidence for ped. However, when we set g = 10n, we see that comparing models via the marginal likelihood suggests that ped is no longer significant, while the aggregate cross-validation score changes little with this increased variance of the prior. As a sanity check, we run a Gibbs sampler targeting π(θ | y 1:n , x 1:n ) for the two prior settings within the full model M 0 , and plot the marginal posterior of θ ped in Figure A.7.1. For reference, the posterior means of θ glu , θ bp are 0.70 and 0.12 respectively. The posteriors of θ ped are indistinguishable for the two prior settings, with a significant mean for θ ped . This agrees well with the aggregate crossvalidation scoreŜ ACV which is clearly robust to vague priors.