PyMC – Calculate Evidence (Bayes' Theorem denominator) 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.
 A: 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.
A: 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.
A: 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.
