PyMC – Calculate Evidence (Bayes' Theorem denominator)

Question

I have one simple model and one more complex model (both are PDEs that describe displacement, when different kind of forces and thermal effects take place). From those, given some data, I calculate through PyMC the posterior distribution of a parameter (elasticity). Now I want to calculate the Bayes' factor, and I am supposed to calculate it by directly calculating the two evidences, but I cannot find any information about how to do this. I use a uniform prior for my parameter (=a parameter of my partial differential equations system), and a normal likelihood with mean that depends on the physical parameter of interest, and precision Uniformly distributed.

Is there a way to calculate the denominator (evidence) of Bayes' formula?

Edit: In general my probabilistic Model-1 looks like this: (parameter of interest: E)

$E = Uniform$ #Prior for Hyperparameter

$\tau = Normal$ #Prior for parameter

$Deterministic: u_1(E)$

$Observed: DataLikelihood: Normal(\mu=u_1(E), \tau=\tau, value=data)$

For Model-2, everything is the same, except for the deterministic u. i.e.

$Observed: DataLikelihood: Normal(\mu=u_2(E), \tau=\tau, value=data)$

EDIT:

I used kernel density estimators scipy.stats.kde.gaussian_kde on the posteriors of the two different models. Now that I have the densities, I can evaluate them at the data points (i.e. find the likelihood for each data point for each model). So by taking the mean of those likelihoods I approximate the marginal likelihoods, right? So by dividing them I get the Bayes' factor. Please correct me if there is a flaw in this procedure.

Edit: I'm starting to believe that this is just model fitness comparison and not pure evidence, since no information about the prior were being used. Any suggestions are more than welcome.

Answer 1

The Bayes factor is the ratio of two Bayesian evidences and must be computed as an integral over the entire parameter space (see e.g. this review).

For some scientific applications I was using MultiNest and its python interface that you may implement in your code, see the examples in the github repository. Even better, since it is a more efficient algorithm, would be to use PolyChord, but the original code is in Fortran. I just discovered that it also has a python interface, but I never used it.

Answer 2

I am not familiar with PyMC but this should solve your problem.

http://stronginference.com/post/bayes-factors-pymc

I would have posted this as a comment but I don't have enough reputation to do that yet!

---

Or you could calculate the bayes factor by hand by calculating

$$\text{BF}_{i,j} = \frac{L(Y \mid M_i)}{L(Y \mid M_j)} = \frac{\int L(Y \mid M_i,\theta_i)p(\theta_i \mid M_i)d\theta}{\int L(Y \mid M_j,\theta_j)p(\theta_j \mid M_j)d\theta}$$

where the numerator and denominator is the marginal likelihood under each model (i.e., the evidence or denominator in Bayes' Theorem).

---

So, without filling in too much of the details, what you are going to need to do is calculate the following:

$$\int_{-\infty}^\infty\int_0^\infty \text{Normal}(\mu(\alpha),\tau^2)\text{Uniform}(\alpha;a,b)\text{Uniform}(\tau^2;c,d)d\tau d\mu$$ but because you are dealing with two uniforms it is going to simply to $$\frac{1}{d-c}\frac{1}{b-a}\int_{-\infty}^\infty\int_0^\infty \text{Normal}(\mu(\alpha),\tau^2)d\tau d\mu$$ So you just need to be able to compute that integral.

Answer 3

I used kernel density estimators scipy.stats.kde.gaussian_kde on the posteriors of the two different models. Now that I have the densities, I can evaluate them at the data points (i.e. find the likelihood for each data point for each model). So by taking the mean of those likelihoods I approximate the marginal likelihoods, right? So by dividing them I get the Bayes' factor. Please correct me if there is a flaw in this procedure.

Edit: I'm starting to believe that this is just model fitness comparison and not pure evidence, since no information about the prior were being used. Any suggestions are more than welcome.