How to use Bayesian Evidence to Compare Models I'm exploring how to use Bayesian Estimation to compare models, and need to compare with a way to implement it. So far I just found this article is easier to understand, but I'm still not sure how to implement. If you go to page 20 "Model comparison using the Bayesian evidence".
We just need to calculate this for each model and compare:
$$P(M_1|D)=\frac{P(D|M_1)P(M_1)}{\sum_iP(D|M_i)P(M_i)}$$
Also as the author mentioned, "$P(M)$ is the model prior probability. If we have no reason to favour one model
over another, then we just set all of these to be equal". $P(D|M_i)$ is the evidence and "The important point is that the evidence is a normalized PDF over the data".
These raised me 2 questions:

*

*$P(M_i)$ can be equal for all the models, it means I can set them all as 1 or the sum of them should be 1?


*If evidence $P(D|M_i)$ is normalized pdf, it should be a distribution, right? Then how can I convert a normalized pdf to a value? What confuses me most is, imagine I have features X, ground truth y, and predicted classes y_pred. How can I calculate $P(D|M_i)$ for each model?
Is that possible to give me some suggestions with a simple example with features X, ground truth y, and predicted classes y_pred? If you can do that in python or R, it's even better.
 A: Let me give a shot at a basic answer (albeit without code)...
Question 1: Yes, it means the sum of them should be 1.
Question 2: $P(D|M_i)$ is indeed a distribution. It provides the probability of observing the data $D$ given the generating underlying process is Model $M_i$. To get a specific value, you need to plugin the model you want to evaluate and the data.
Let's make an example: Assume your data $D$ consists of points $(x, y)$ distributed (modeled) according to a straight line with zero-mean i.i.d. Gaussian residuals:


*

*Data $D = \{(1, 2), (2, 4.1), (3, 6)\}$

*Statistical Model $M = (S, \mathcal{P}) = (\mathbb{R}, \text{straight line with i.i.d. Gaussian residuals})$

*i.e. $Y = \beta X + \epsilon$ with $X \sim \mathcal{N}(0, 1)$ and $\epsilon \sim \mathcal{N}(0, \sigma^2)$

*i.e. $\mathcal{P} = \left\{ P(x,y) = P(y|x) P(x) = \mathcal{N}(\beta X, \sigma^2) \cdot \mathcal{N}(0, 1) : \beta \in \mathbb{R}, \sigma > 0\right\}$

*explicitely written $$P(x,y) = \mathcal{N}(\beta X, \sigma^2+1) = \frac{1}{2\pi\sigma} e^{-\frac{(y-\beta x)^2}{2\sigma^2}} e^{-\frac{x^2}{2}}$$


Now let's assume you used some method to estimate the model parameters and found a candidate you think is good. Let's call that model $M_1$ and the parameters you've found are $\beta = 1$ and $\sigma = 1$, i.e. $P(x,y) = \frac{1}{2\pi} e^{-\frac{(y-x)^2}{2}} e^{-\frac{x^2}{2}}$
Based on this calculate $P(D|M_1)$:
\begin{align}
P(D|M_1) 
&= P(\{(1, 2), (2, 4.1), (3, 6)\}| \{P(x,y) = \mathcal{N}(X, 1) \cdot \mathcal{N}(0, 1)\})\\
&= P((1, 2)|P(x,y)) \cdot P((2, 4.1)|P(x,y)) \cdot P((3, 6.1)|P(x,y)) \\
&=
\frac{1}{2\pi} e^{-\frac{(2-1)^2}{2}} e^{-\frac{1^2}{2}} \cdot
\frac{1}{2\pi} e^{-\frac{(4.1-2)^2}{2}} e^{-\frac{2^2}{2}} \cdot
\frac{1}{2\pi} e^{-\frac{(6-3)^2}{2}} e^{-\frac{3^2}{2}} \\
&=
\frac{1}{2\pi} e^{-\frac{1}{2}} e^{-\frac{1}{2}} \cdot
\frac{1}{2\pi} e^{-\frac{4.41}{2}} e^{-\frac{4}{2}} \cdot
\frac{1}{2\pi} e^{-\frac{9}{2}} e^{-\frac{9}{2}} \\
&=
\frac{1}{2\pi} e^{-1} \cdot
\frac{1}{2\pi} e^{-4.205} \cdot
\frac{1}{2\pi} e^{-9} \\
&= \frac{1}{8 \pi^3} e^{-14.205}
\end{align}
As you can see the "converting it to a value" depends on the specific model (including parameters) and the observed data.
