# Bayesian Justification of Cross-validation

If I understand correctly, K-fold cross-validation is supposed to approximate expected log predictive density (ELPD), which is defined as $$\mathop{\mathbb{E}}_{D_{new}\sim P(.|M_{true})}\log P(D_{new}|D_{observed},M)$$, for a model $$M$$, and assuming there is a "true" data generating process $$M_{true}$$ (that we do not know what it is). This is used as a measure of how good model $$M$$ is. Regarding this, I have two questions:

1. In this expectation, why do we use the logarithm of the likelihood of the new data, i.e. $$\log P(D_{new}|D_{observed},M)$$ and not the new posterior, i.e. $$\log P(M|D_{new},D_{observed})$$. How is the former "better", whatever that's supposed to mean? Is there any source that has used the latter?

2. More generally, what is the mathematical justification that ELPD, or its approximation using cross-validation, is a better measure for model selection than pure posterior calculated using all of our data? Bayesian posterior is uniquely derived from a set of coherency criteria and any other measure is strictly inferior to it (at least when we are only concerned with those coherency criteria). I suppose the usefulness of ELPD and cross-validation is usually manifested when we cannot quantify our priors well enough. But this understanding I have is purely based on intuition, and I do not even trust it that much. Is there a more rigorous treatment of this? For example, is it possible to show, mathematically, that ELPD is better, on some scale, than posterior calculated with an inaccurate prior?

Edit: Here is my failed attempt at answering this question.

• Regarding (1): what do you propose as the prior $p(M)$? If it's anything other than uniform, I'd be curious to hear it; however, if it's uniform, consider the implications for the relationship between the two expectations. Commented Apr 22 at 21:12
• Regarding (2): the posterior is often intractable, so we approximate it, and this is one avenue for doing so. Commented Apr 22 at 21:22
• Regarding your first comment, regardless of the fact that sometimes we do have useful priors and we will use that, even if we use a flat prior the implication is weighting the likelihood of the new data with the likelihood of the old data. It would not make them equal.
– Feri
Commented Apr 22 at 21:44
• Regarding your second comment, I don't think intractability is the only reason this method is used. For one thing, even in cross-validation we will still need to calculate $P(D_{new}|D_{obs},M)$ and that would also be intractable so we are not solving anything.
– Feri
Commented Apr 22 at 21:50
• You may want to read "The fallacy of sample re-use" in Jaynes' book "Probability Theory, the Logic of Science" (p.264). Commented Apr 26 at 7:40

## 1 Answer

Several ideas could serve to justify model selection based on the ELPD instead of the model-posterior.

1. Direct relation to frequentist-like MSE comparison for regression models with Gaussian likelihood and Gaussian prior. Consider the following regression models ($$\sigma^2$$ known):

Model $$M_1$$: \begin{align} y_i &\overset{iid}{\sim} N(\mu,\sigma^2), \quad i=1,\cdots,n \\ \mu &\sim N(\mu_0,\sigma^2_0). \end{align}

Then, the log of the predictive distribution is: $$\log p(\tilde{y}\mid D,M_1) = -\frac{1}{2}\log\left[2\pi(\sigma^2+\sigma^2_1) \right] -\frac{1}{2}\frac{(\tilde{y}-\mu_1)^2}{\sigma^2+\sigma^2_1},$$

where $$\mu_1$$ and $$\sigma^2_1$$ are the posterior mean and posterior variance of $$\mu$$.

Model $$M_2$$, with $$x_i\in\mathbb{R}^d$$:

\begin{align} y_i\mid x_i &\overset{iid}{\sim} N(x_i^\top\beta,\sigma^2), \quad i=1,\cdots,n \\ \beta &\sim N_d(\beta_0,\Sigma_0). \end{align}

In this case, the log of the predictive distribution is: $$\log p(\tilde{y}\mid \tilde{x},D,M_2) = -\frac{1}{2}\log\left[2\pi(\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}) \right] -\frac{1}{2}\frac{(\tilde{y}-\tilde{x}^\top\beta_2)^2}{\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}},$$

where $$\beta_2$$ and $$\Sigma_2$$ are the posterior mean and posterior covariance matrix of $$\beta$$.

When $$n$$ is large, $$\sigma^2+\sigma^2_1\approx\sigma^2$$ and $$\sigma^2+\tilde{x}^\top\Sigma_2\tilde{x}\approx\sigma^2$$. This is because, under no misspecification, the posterior variances $$\sigma_1^2$$ and $$\Sigma_2$$ shrink to zero with $$n$$. Thus, the ELPD of each model is approximately:

\begin{aligned} \mathbb{E}_{\tilde{y}} \log p(\tilde{y}\mid D,M_1) &\approx -\frac{1}{2}\log\left(2\pi\sigma^2 \right) -\frac{1}{2\sigma^2}\mathbb{E}_{\tilde{y}}(\tilde{y}-\mu_1)^2 = K_1-K_2^2\,\text{MSE}_1,\\ \mathbb{E}_{\tilde{y},\tilde{x}} \log p(\tilde{y}\mid\tilde{x}, D,M_2) &\approx -\frac{1}{2}\log\left(2\pi\sigma^2 \right) -\frac{1}{2\sigma^2}\mathbb{E}_{\tilde{y},\tilde{x}}(\tilde{y}-\tilde{x}^\top\beta_2)^2 = K_1-K_2^2\,\text{MSE}_2. \end{aligned}

So, in this context, and for a large sample size $$n$$, choosing the model with highest ELPD would mostly be equivalent to choosing the model with the lowest MSE. This is, selecting/comparing Bayesian models using the ELPD finds a direct analogue in the frequentist framework using the MSE. For other models beyond the Gaussian likelihood + Gaussian prior, the connection can be given in terms of different loss functions.

1. The ELPD is directly related to several criteria from Information Theory and, in particular, to those based on Kullback-Leibler (KL) information. The KL divergence from the true data-generating distribution $$p^*$$ to the (marginal) predictive distribution in model $$M$$ is:

$$D_{KL}(p^*\,|\!|\,p_{\text{pred}})=-\mathbb{E}_{\tilde{y}} \log p(\tilde{y}\mid D,M)-\mathbb{H}(p^*)\geq 0,$$

where $$\mathbb{H}(p^*)$$ is the entropy of $$p^*$$. Most information criteria measures can be expressed as -2 times an estimated total ELDP, $$\hat{L}=\sum_{i=1}^n\hat{L}_i$$, in combination with a penalization term, $$2\hat{\Lambda}(M)$$. The latter is required to account for overfitting when no holdout data is available for evaluation.

For instance, for a finite-dimensional parametric model $$M=M_\theta,\,\theta\in\mathbb{R}^k,\, k\in\mathbb{N}$$, with $$m$$ samples from the posterior $$p(\theta\mid D,M)$$, one can compute the Watanabe-Akaike information criterion (Gelman's definition) as $$\text{WAIC} =-2\cdot(\hat{L}_W-2\hat{\Lambda}_W)$$, with:

$$\hat{L}_W = \sum_{i=1}^n\log\left(\frac{1}{m}\sum_{j=1}^m p(y_i\mid \theta_j) \right),\quad \text{ and }\quad \hat{\Lambda}_W=\hat{L}_W-\sum_{i=1}^n\frac{1}{m}\sum_{j=1}^m\log p(y_i\mid \theta_j) .$$

The WAIC employs a posterior-sample-based approximation of the total ELPD. For this reason, it is considered a more fully Bayesian approach than the Akaike IC (AIC) or the Deviance IC (DIC). Model selection using the WAIC is asymptotically equivalent to model selection using LOO-CV. Furthermore, if the true model is also contained in the family of candidate models $$\mathcal{M}$$, both are also asymptotically equivalent to model selection using the model-posteriors:

$$p(M\mid D) = \frac{p(M)\,\int_\theta p(\theta\mid M)\,\prod_{i=1}^n p(y_i\mid\theta,M)\, \text{d}\theta}{{\sum_{\tilde{M}\in\mathcal{M}} p(\tilde{M})\,\int_\theta p(\theta\mid \tilde{M})\,\prod_{i=1}^n p(y_i\mid\theta,\tilde{M})\, \text{d}\theta}},\,\quad M\in\mathcal{M}.$$

Yet, when the true model is not in the candidate model space $$\mathcal{M}$$, these approaches might not agree. This is because, in such a case, the candidate model with the lowest KL divergence from the true generating distribution might not be the same model with the highest posterior within $$\mathcal{M}$$.

So, the ELDP also helps bridging Bayesian analysis with model selection based on Information Theory and Information Geometry, and it helps defining "the best" candidate model as the one in $$\mathcal{M}$$ with the lowest KL divergence from the true data-generating distribution. For some purposes, this optimal projection is just enough.

1. The issues with prior specification and high-dimensional model spaces.

Computing $$p(M\mid D)$$ for all $$M\in\mathcal{M}$$ has many practical and computational drawbacks. To start, it requires the specification of model-priors $$p(M)$$. If $$\mathcal{M}$$ is composed of a finite number of candidates, then an uninformative model-prior can be employed: $$p(M)=|\mathcal{M}|^{-1},\ \forall M\in\mathcal{M}$$. Other model-prior specifications can be more difficult to motivate.

More importantly, the model-posteriors depend sensibly on how the specification of parameter-priors $$p(\theta\mid M_k)$$ change across models $$M_k\in\mathcal{M}$$. If models are highly heterogeneous, then the specification of $$p(\theta\mid M_k)$$ should reflect model-adaptive prior-uncertainty. But, $$\theta$$ is presented as model-shared parameters, and so different specifications may lead to different model-posteriors $$p(M\mid D)$$.

Finally, as $$p(M\mid D)$$ is probably not analytically available, a Monte Carlo approach would be required. This, in itself, presents several computational challenges. It is particularly tedious when model space $$\mathcal{M}$$ is discrete, large, and with a high number of posterior modes. For instance, take the case of $$M$$ representing a directed acyclic graph (DAG) on a set of $$d\in\mathbb{N}$$ variables in the data, with no structural constrains (such as prior time-ordering). For a large $$d$$, computing the MAP model is NP-hard, as the number of DAGs is super-exponential, $$O\left(2^{d^2}\right)$$, and thus exploring the entire model space for sampling would be practically unfeasible.

Clearly, some of these troubles would be avoided by employing ELPD-based model comparison.

I share some references on the topic:

Burnham, K. P., & Anderson, D. R. (Eds.). (2002). Model selection and multimodel inference: a practical information-theoretic approach. New York, NY: Springer New York.

Watanabe, S., & Opper, M. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of machine learning research, 11(12).

Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and computing, 24, 997-1016.

Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using Stacking to Average Bayesian Predictive Distributions (with Discussion). Bayesian Analysis, 13(3), 917-1007.

• Thanks! This is an elaborate answer and the bounty is ended so I award the points to you. But I cant accept it because my question remains: consider the scenario where the posterior is analytically tractable (your point 3 wouldnt apply). Cox theorem says Bayes posterior is the only reasonable way to assign plausibility to models. We shouldnt care about compatibility to any frequentist or information theoretical framework (points 1 and 2). But even in this situation people do crossvalidation e.g. to reduce overfitting, sth that intuitively makes sense but I cant find Bayesian justification for.
– Feri
Commented May 2 at 21:52
• I doubt you will. Cross-validation's advantage is that it does not assume your model specification is correct, whereas if you knew with 100% certainty that it was correct, the Bayes posterior would contain all the relevant information and cross-validation would be pointless, as Cox implies. Commented Jun 5 at 15:11
• @jbowman That is my understanding, as I also said in the question. I am looking to formalize this understanding, e.g. to quote the from question: "to show, mathematically, that ELPD is better, on some scale, than posterior calculated with an inaccurate prior". Also, "correct" is a vague word to me. In a way, all models are incorrect. You are looking for an explanation of the data within a model. If using half the data cannot train a model that works on the rest, why does it imply that the model is "incorrect", and not that half your data just isn't representative enough?
– Feri
Commented Jun 5 at 15:56
• Well, you can always construct a sufficiently inaccurate prior to wreck any analysis! But when you are working within a particular framework, e.g., Bayesian, naturally the mathematical results will be predicated on the correctness of the framework. If you want to be a "Bayesian with a badly misspecified likelihood function", you can do better if you don't make your decisions purely within the context of the Bayesian framework. This is addressed in point 2 of the answer above. There is a big element of risk management here. Also... models are correct or incorrect regardless of the data. Commented Jun 5 at 16:24
• @jbowman correctness of a model is independent of the data but our inference about it's correctness obviously depends on the data!
– Feri
Commented Jun 5 at 16:53