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

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 xi) and a scalar output vector (the yi). I am looking to fit a beta vector, i.e. βT * xi = yi + ε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, xi's):

So, it would seem in order to do Gibbs sampling I simply iterate through β1, γ1, β2, and γ2, 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 xi's on the Wikipedia page)??? The way I see it, it could be one of two things. For discussion purposes, let's discuss β1 for now.

1. First thought is, I'll just iterate over each i, generating my data points using (yi - β2*xi,2) / xi,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 β1 and γ1.

2. Second thought is, I'll just go through each i, and generate my data points as (yi / xi), 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?

## migrated from stackoverflow.comAug 15 '11 at 9:14

This question came from our site for professional and enthusiast programmers.

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.

• Thanks so much for answering. I thought if I did derive it by hand (like I used to do for prob sets in college) I would just end up with those formulas for the Wikipedia page, but with the data transformed in some manner. I was trying to circumvent the derivation step by transforming the data directly and using those formulas. Guess I will have to suck it up and get down to business on that. Will report back after I do - but one question first. In part #1, you mention a variable tau. I assume that is what I am supposed to sample from the gamma distribution? How do I derive the form of that?Thx – Jason Aug 12 '11 at 20:53
• Precision is the inverse of variance. Just a different name. As you've seen on Wikipedia page, it is easier to work with precisions, as they are additive as far as computing the posteriors are concerned. You could check your calculations by having a regression with a single explanatory variable, an intercept. Then your formulae should reduce to the ones on Wikipedia page. – StasK Aug 12 '11 at 21:54