Bayes Factor approximation

Question

A brute force method to approximate the Bayes Factor (the ratio of the denominators (normalizing constants) in the Bayes formula) is to do the following for the two models of interest:

repeat multiple times

• draw parameters from the prior density
• compute the likelihood given those parameters.

If you then take the average of all those likelihoods for both models and compute their ratio, you get an approximation of the BF.

My questions are:

1. Is this the correct way to approximate the BF (disregarding the fact that there are more efficient ways such as Importance Sampling)?
2. Imagine I want to repeatedly compute the BF as the sample size increases so that my initial priors are replaced by the posteriors after each updating step. Is it still a correct approximation if I use posterior samples as the "prior" density. Put differently, if I don't sample the parameters from a proper prior/posterior density function but rather from a vector with samples from the prior/posterior.

Answer 1

The answer to your first question is "yes": Using notation similar to Wikipedia, the Bayes factor is $$ K = \frac{P(D|M_1)}{P(D|M_2)} = \frac{\int P(D|\theta_1,M_1) P(\theta_1|M_1) \,d\theta_1} {\int P(D|\theta_2,M_2) P(\theta_2|M_2) \,d\theta_2}, $$ and the procedure you describe just replaces both numerator and denominator with Monte Carlo estimates: for the numerator we get $$ \int P(D|\theta_1,M_1) P(\theta_1|M_1) \,d\theta_1 = E_{M_1} \bigl( P(D|\theta, M_1) \bigr) \approx \frac{1}{N} \sum_{j=1}^N P(D|\theta_1^{(j)}, M_1), $$ where the subscript $M_1$ at the expectation indicates that $\theta$ in the expectation is distributed according to (the prior of) model 1 and, correspondingly, the $\theta_1^{(j)}$ on the right-hand side are independent and identically distributed samples from the prior of model 1. Using the same approximation for the denominator, we get $$ K \approx \frac{\frac{1}{N} \sum_{j=1}^N P(D|\theta_1^{(j)}, M_1)} {\frac{1}{N} \sum_{j=1}^N P(D|\theta_2^{(j)}, M_2)}. $$ The law of large numbers guarantees that this estimate converges to the correct value $K$ as the sample size $N$ increases.

There are a few things to consider:

1. If $P(D|M_2)$ is small, the estimate will likely perform badly: small fluctuations in the estimate for the denominator will turn into large fluctuations of the estimate for $K$.

2. Variance reduction methods (like importance sampling) can be used to make the estimates for the denominator and numerator more efficient.

3. The estimates for the denominator and numerator are unbiased, but the estimate for $K$ is not.

I believe that more detail about the adaptive version is required before a definite answer to your second question can be given.

Answer 2

The answer to your second question is yes as well if you consider the sequence of Bayes factors: $$\mathfrak{B}_{12}^t=\dfrac{\int f_1(x_1,\ldots,x_t|\theta_1)\pi_1(\theta_1)\text{d}\theta_1}{\int f_2(x_1,\ldots,x_t|\theta_2)\pi_2(\theta_2)\text{d}\theta_2}\qquad t=1,\ldots,T$$ Then, assuming the observations are iid, \begin{align}\mathfrak{B}_{12}^{t+1}&=\dfrac{\int f_1(x_{t+1}|\theta_1)f_1(x_1,\ldots,x_t|\theta_1)\pi_1(\theta_1)\text{d}\theta_1}{\int f_2(x_{t+1}|\theta_2)f_2(x_1,\ldots,x_t|\theta_2)\pi_2(\theta_2)\text{d}\theta_2}\\ &=\dfrac{\int f_1(x_{t+1}|\theta_1)\dfrac{f_1(x_1,\ldots,x_t|\theta_1)\pi_1(\theta_1)}{\int f_1(x_1,\ldots,x_t|\theta_1)\pi_1(\theta_1)\text{d}\theta_1}\text{d}\theta_1}{\int f_2(x_{t+1}|\theta_2)\dfrac{f_2(x_1,\ldots,x_t|\theta_2)\pi_2(\theta_2)}{\int f_2(x_1,\ldots,x_t|\theta_2)\pi_2(\theta_2)\text{d}\theta_2}\text{d}\theta_2}\times \mathfrak{B}_{12}^t\\ &=\dfrac{\int f_1(x_{t+1}|\theta_1)\pi_1(\theta_1|x_1,\ldots,x_t)\text{d}\theta_1}{\int f_2(x_{t+1}|\theta_2)\pi_2(\theta_2|x_1,\ldots,x_t)\text{d}\theta_2}\times \mathfrak{B}_{12}^t\\ \end{align} This representation means that if

1. You already have an approximation of $\mathfrak{B}_{12}^t$, $\hat{\mathfrak{B}}_{12}^t$;
2. You have samples from both posteriors $\pi_1(\theta_1|x_1,\ldots,x_t|\theta_1)$ and $\pi_2(\theta_2|x_1,\ldots,x_t)$;

you can derive an approximation of the first ratio using those samples and hence deduce an approximation of $\mathfrak{B}_{12}^{t+1}$. However you have to update your samples at each step $t$ in order to get samples from both posteriors $\pi_1(\theta_1|x_1,\ldots,x_{t+1})$ and $\pi_2(\theta_2|x_1,\ldots,x_{t+1})$. Which is why I do not believe this form of update is particularly useful, as indicated in a question of mine's.