# Bayes factors from MCMC samples

I'm working to implement Bayesian model selection among models whose posteriors have already been sampled via MCMC. After reviewing some discussions of Bayes factors, I understand that they are sensitive to choice of prior, that their computation via the harmonic mean of marginal likelihood samples is treacherous, and that alternate information criteria are available.

[ $$y$$=data, $$M$$=model, $$\vec{\theta}$$=parameters, $$V_{\vec{\theta}}$$=volume of parameter space, $$N$$=# MCMC samples ]

Question 1)

The Bayes factor is the ratio of marginal likelihoods of two models, $$\frac{\text{p}(y|M_0)}{\text{p}(y|M_1)}.$$ For a given model, is it OK to think of the marginal likelihood as an integral over the parameter posterior? I.e. $$\begin{split} \text{p}(y|M_i) & = \int \text{p}(y|\vec{\theta},M_i) \text{p}(\vec{\theta},M_i)\text{d}\vec{\theta} \\ & \propto \int \text{p}(\vec{\theta}|y,M_i)\text{d}\vec{\theta} \text{ ?} \end{split}$$

Question 2)

If my understanding in Question 1 is correct, when would it be safe to calculate $$\text{p}(y|M_i)$$ using MC approximation of $$\int \text{p}(\vec{\theta}|y,M_i) \text{d}\vec{\theta}$$ ? If my MCMC procedure has already converged upon the full parameter posterior, and each parameter's marginal posterior appears well-sampled, can I reliably calculate $$\text{p}(y|M_i)$$ using $$\text{p}(y|M_i) \propto \int \text{p}(\vec{\theta}|y,M_i) \text{d}\vec{\theta} \approx \frac{V_{\vec{\theta}}}{N} \sum_{j=1}^{N} \text{p}(\vec{\theta_j}|y,M_i) \text{ ,}$$ where $$\text{p}(\vec{\theta_j}|y,M_i)$$ are the MCMC samples?

Computation of the marginal likelihood from MCMC samples seems like a detailed response, but I'd like to know if I'm on track with my initial concept of what goes into a calculation for $$\text{p}(y|M_i)$$ from MCMC samples.

## 1 Answer

Regarding 1, your expression needs to be adjusted a little: $$\begin{split} \text{p}(y|M_i) & = \int \text{p}(y,\vec{\theta} | M_i)\text{d}\vec{\theta} \\ & =\int \text{p}(y|\vec{\theta},M_i) \text{p}(\vec{\theta} |M_i)\text{d}\vec{\theta} \\ &= E_{\vec{\theta} |M_i}[\text{p}(y|\vec{\theta},M_i)]. \end{split}$$ Your $$\int \text{p}(\vec{\theta}|y,M_i)\text{d}\vec{\theta}$$ is actually just equal to $$1$$ (because it's the integral of a density), and your $$\int \text{p}(y|\vec{\theta},M_i) \text{p}(\vec{\theta},M_i)\text{d}\vec{\theta}$$ is equal to $$\text{p}(y,M_i)$$ (not $$\text{p}(y|M_i)$$).

This does not use draws from the posterior, though. For that you would need a different expression. This is explained in the thread you linked to.

Regarding 2, yes, that's quite a meaty thread. I'd start there. You might be trying to write $$E_{\vec{\theta} |M_i}[\text{p}(y|\vec{\theta},M_i)] \approx N^{-1}\sum_{j=1}^N\text{p}(y|\vec{\theta}^j,M_i) \tag{1}$$ where $$\vec{\theta}^j$$ are samples from the prior. That's one version of a Monte-Carlo or importance sampling estimate.

Edit: as was mentioned in the comments, you might consider using another importance sampling estimator, one with a different proposal/instrumental density: $$N^{-1}\sum_{j=1}^N\frac{\text{p}(y|\vec{\theta}^j,M_i) \text{p}(\vec{\theta} |M_i)}{ \text{q}(\vec{\theta})} \tag{2}.$$

• Yes, with an importance weight, another distribution works better. The only one that does not work is the posterior itself (harmonic mean estimator warning!). – Xi'an Mar 5 at 19:10
• @Taylor, each step in my MCMC chain represents a $\vec{\theta}^j$ (drawn from joint prior) for which $p(y|\vec{\theta}^j,M_i)$ was calculated, so I should be able to use those existing samples to approximate $E_{\vec{\theta}|M_i}[\text{p}(y|\vec{\theta},M_i)]$. Regarding my second question, if my MCMC procedure has already converged upon the joint parameter posterior, can I be confident in going ahead and using those converged $p(y|\vec{\theta}^j,M_i)$ samples in calculating $E_{\vec{\theta}|M_i}$? – curiousStudent Mar 5 at 19:51
• @Xi'an I apologize I was thinking $\text{p}(\vec{\theta} |M_i)$ is the posterior instead of the prior(!) – Taylor Mar 5 at 20:23
• @curiousStudent samples are from the prior in estimate (1). In (2), you could set $q$ equal to the posterior, and evaluate that fraction for each sample, but that would probably give you a high-variance estimate. – Taylor Mar 5 at 20:25