Gibbs sampling for a simple linear model — need help with the likelihood function

Question

So in order to better acquaint myself with Gibbs sampling, I've been working on a fairly simple linear model, written in Python/R. Basically, I have 2-dimensional input data (the x_i) and a scalar output vector (the y_i). I am looking to fit a beta vector, i.e. β^T * x_i = y_i + ε_i (ε_i is the noise).

So I decided I'd use a Gaussian prior for the betas, plus an inverse gamma prior for their variance, giving me a posterior distribution that is Gaussian (thanks conjugacy!). And so I used the formula from the Wikipedia page on conjugate distributions to get the likelihood function I need to generate random samples of my beta coefficients and their variances, given the data (by their notation, x_i's):

So, it would seem in order to do Gibbs sampling I simply iterate through β₁, γ₁, β₂, and γ₂, generating samples from the distributions listed on that wikipedia page. My confusion is, what exactly are my data points (the "samples" I'm using to calculate the distribution parameters, which I need to sample my own betas, i.e. the x_i's on the Wikipedia page)??? The way I see it, it could be one of two things. For discussion purposes, let's discuss β₁ for now.

1. First thought is, I'll just iterate over each i, generating my data points using (y_i - β₂*x_i,2) / x_i,1 — basically, I subtract out the influence of the second factor and divide the "leftovers" by the first factor, to obtain the impact of the first factor on the response variable. Then I can just use that population to get my distribution parameters, and finally can sample my β₁ and γ₁.

2. Second thought is, I'll just go through each i, and generate my data points as (y_i / x_i), meaning I do NOT subtract out the influence of the other factors on the response variable; it treats everything independently.

So far I've been using #1, but even with artificial data sets, I'm finding that it is implying enormous beta variances that don't seem to make any sense — BUT, they are darn close to the sample variances in the populations I'm generating (my "data points" I use to calc those distribution parameters). What am I doing wrong? Is #1 or #2 the right way to do a linear model? Or am I missing something entirely?

Thanks in advance for your help! Please let me know if any info would help.

Answer 1

This is a statistics question, not a programming question, and would better be asked on CrossValidated. At least, the LaTeX code is getting parsed there automatically :). Also, this is more complicated than what is readily available on that webpage. I'll give some guidance, but as long as you want to learn how to do things, this won't be the complete answer. (If you don't want to do that, we can locate the cooked answers on the web, too.)

Each sampling of betas relies on the complete data set. If you are doing this with individual y_i's and x_i's, you are not doing this right. Before you start working with the code, you need to sit down with a piece of paper (letter size or A4, depending on your geography) and derive the posterior distributions of betas:

1. This is given: y|beta is normal with mean x'beta and precision tau
2. This is given: prior for beta is normal with mean mu and precision gamma
3. Obtain this: the marginal distribution of y, by integrating the betas out (which is easy to do, since the joint distribution of y and beta is multivariate normal, and you can do this by kernel matching: the part that depends on beta is going to be exp[ a quadratic form in beta], so you recognize this to be a relevant part of a normal distribution distribution to integrate over; whatever's left after integration should be a normal density in y and the prior parameters)
4. Obtain this: the posterior distribution of beta given y, by Bayes theorem (the likelihood times the prior divided by the posterior; again this should be a moderately complicated combination of exp[ jointly quadratic in y and beta ])
5. Obtain this: the conditional distribution of beta_1 given beta_2 and y, one of the margins of the multivariate normal distribution obtained at the previous step.

You need to know how to manipulate the multivariate normal distribution and get conditional and marginal distributions out of it. Again, if this is over your head, we can find the ready solutions.

Note that you also need a sampler for the variance of regression errors, unless you treat it as known (which is hardly a practical situation). This will be slightly more complicated, as you would need to incorporate another dimension into your integration procedures.