# 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.

• My prior leads me to believe this probably can be done. Can you add more detail? – gung - Reinstate Monica Aug 21 '14 at 18:21
• I have one simple model and one more complex model (both are pde's 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 Uniform prior for my parameter and Normal likelihood. – Low Yield Bond Aug 21 '14 at 18:29
• I am not 100% sure, but I think there may be problems with using uniform priors. Not from a computational point of view but more from a is this a good choice when calculating Bayes Factors. – Dan Aug 21 '14 at 19:12
• You already have the prior, the likelihood, and the posterior. Can't you just divide your posterior by the prior and the likelihood to obtain the normalization constant / evidence / denominator? The whole point of posterior sampling is that the denominator is hard/impossible to compute. But that means there's more than one way to obtain an (approximate) posterior. As far as I can tell you can just plug that into the equation. This is as much a question for @gung as for anyone else. – shadowtalker Aug 24 '14 at 12:41

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.

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.

• Thank you. I have studied this page, but I have instructions that I should do it by calculating the evidences. To be honest, I don't completely understand the way that model selection is done in this page, and as a result I wasn't able to persuade that a similar way is the right way to go. – Low Yield Bond Aug 21 '14 at 18:41
• Can you calculate the evidence by hand? – Dan Aug 21 '14 at 18:42
• Also, you say you have a uniform prior but what parameter is that prior for? The Normal distribution has two parameters so which parameter are you referring to (the mean or the standard deviation)? – Dan Aug 21 '14 at 18:45
• Based now on your edits, do you think you could explain how the parameter from your partial differential equations system enters into the likelihood? We usually only place priors on parameters of the likelihood (or functions of it) so without seeing your models more explicitly I'm not sure an answer exists without just saying calculate the integral by hand for Bayes' Theorem denominator. – Dan Aug 21 '14 at 18:54
• I have set the precision (1/variance) as a Uniform, but I am interested in a (physical) parameter in the original PDEs system (which I also set as Uniform), which is linked deterministically to the mean of the likelihood. Therefore I am little confused and can't think a way to approach this problem by hand. – Low Yield Bond Aug 21 '14 at 18:55

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.