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.
  • 1
    $\begingroup$ Your second question is not clear. Specifically what Bayes Factor (BF) are you trying to approximate? As worded, it seems as though you have a sequence of BFs. $\endgroup$
    – jaradniemi
    Commented Mar 29, 2016 at 14:53
  • $\begingroup$ I hope this is more clear: Is it still valid If I can't sample the parameters from a known density function but only have samples available from Gibbs/MCMC sampling (stored in a vector). So in order to draw a paramter I cannot use the r... function in R such as rnorm() etc but instead have to pick the parameters randomly from this vector. $\endgroup$
    – beginneR
    Commented Mar 29, 2016 at 14:59
  • 1
    $\begingroup$ Actually according to this $p$-values are consistent with BF, so $p$-values could serve as an approximation ;) $\endgroup$
    – Tim
    Commented Nov 2, 2016 at 8:15

2 Answers 2


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.

  • 2
    $\begingroup$ (+1) While this particular importance sampler works in principle, there exist more efficient approaches like bridge sampling (Gelman and Meng, 1998). $\endgroup$
    – Xi'an
    Commented Oct 16, 2016 at 17:09

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.