# Real calculation of Bayesian model selection

I've been reviewing Bayesian model selection, where probability of data given the model is defined with the following equation: $$p(y|M_k) = \int p(y|\theta)p(\theta|M_k)d\theta$$ As I understand, $p(y|\theta)$ is simply the distribution of y under the parameters we've already defined. I'm a bit of confused by the term $p(\theta|M_k)$. I know that we can rewrite it with the Bayes rule: $$p(\theta|M_k) = p(\theta)p(M_k|\theta)$$ I suppose that we can find the value of $p(\theta)$ using the Beta probability distribution. But what is happening with $p(M_k|\theta)$? How do we estimate the probability of the model given the parameters? It feels like we should estimate the probability of the probability distribution given its parameters, so maybe it is a Dirichlet distribution problem?

I don't know what to do with this $p(M_k|\theta)$ term, will be very grateful if you can explain it to me.

It seems to me that each of the models you are comparing also come with a prior that is specific to that model. In a simple regression setting, for example, you are trying to choose between explaining $y$ with:

Model 1 ($M_1$), which has a single regressor $x_1$ and priors $f_1(\beta)$ and $g_1(\sigma)$. These priors are your $p(\theta|M_1$), with $\theta=(\beta,\sigma)$.

Model 2 ($M_2$), which has a single regressor $x_2$ and priors $f_2(\beta)$ and $g_2(\sigma)$. These priors are your $p(\theta|M_2$), with $\theta=(\beta,\sigma)$.

You can find a nice example in slides 42-46 here: https://drive.google.com/file/d/0B081GdveJIEkVmwzTUJzdk15WE0/view?usp=sharing

• Thank you for the answer. Do you mean that we can not decompose probability of theta given model with Bayes, but just directly refer to the model we have to compute this value? Like not to calculate probability of theta and probability of model given theta separately? – olejnik_ Dec 14 '16 at 0:46
• That is correct. $p( \theta|M_1)$ and $p(\theta|M_2)$ are the priors that you determined for each of the models you are comparing. They are something you already have. – suckrates Dec 14 '16 at 1:19
• No need to decompose them or try to derive them some other way. – suckrates Dec 14 '16 at 1:25
• ok, I got it. Thank you very much for explanation. – olejnik_ Dec 14 '16 at 19:25

On the issue of model priors:

It is usually impossible to assign a specific prior to each of the $2^p=K$ models, where $p$ is the number of predictors. A conventional way to go about this problem is to elicit a non-informative prior to the models as it is hard to know anything specific about the model space. However, if you are certain about the probability inclusion of some variables then it is possible to convey this information in the set up.

The idea is to decentralise prior mass such that the MCMC can discover the majority of all good models. An example of a conventional set up for the model prior is as follows:

$\textrm{for}\quad M_k\in\mathcal{M} \quad \textrm{let} \quad M_k\sim{}Bin(p,\phi)\quad \textrm{where} \quad \phi\sim{}Beta(a,b)$

By letting $a=1$ you get $b=\frac{p-\mathbb{E}(M_k)}{\mathbb{E}(M_k)}$, thus one has only to specify the expected prior models size, although one wants this expectation to carry as little influence as possible.

You can use a tessellation model prior and a ridge prior for the coefficients if your design matrix is singular, that way you'll avoid centring prior mass around bad models, making it easier for the MCMC to discover good models.

• Thank you very much for your answer, but as long as I'm not very familiar with MCMC it makes the answer a little hard for me to understand. First of all, what is the expected value of the model? As I understand, if our model is probability distribution, than it is the mean, but if our model is from ML, what value does expectation take? Average assigned value by ML algorithm? And did I understand it right: you're assuming that we have 2^p models because each of the variables/predictors can be one of the 2 models we're working with? – olejnik_ Dec 15 '16 at 17:30
• Im not exactly sure what you mean by ML? All models are to be assigned a probability distribution unless you are performing a non-parametric analysis, which complicates things in the Bayesian world. And yes, you are right that it is the expectation taken from the specified distribution, although it can be user specified where, as I mentioned, it shouldn't carry much influence in the end. The $2^p$ models represent all possible model combinations, where in the case of two predictors there are 4 different combinations that should be tested. – Herzriesig Dec 15 '16 at 22:58
• ok, thank you very much, I got what you mean. You made it easier for me to test models without referring to particular specified probability distributions parameters assigned for each model. – olejnik_ Dec 16 '16 at 12:46
• No problem! If you are up for frequentist model selection then a LASSO/Elastic Net approach is an excellent way to go and to some extent easier to grasp. If you want to know about MCMC then start with the ergodic theorem and get a feel of the Metroplis Hastings Algorithm – Herzriesig Dec 16 '16 at 13:53