lm.bayes returns N draws from the posterior distribution of the coefficients and the variance using an MCMC algorithm (with Gibbs steps).  QR factorization is used for stability, instead of matrix inversion. 

### input: vector y of responses
### input: X the design matrix  (n by p, this can easily be constructed with lm function, see the example).  X should contain a column of ones if an intercept term is desired.
### input: N - a value of how many posteriors draws you want
### output: posterior draws for b0,b1,...,bp, sigma^2

lm.bayes=function(y,X,N=10000)
{    n=length(y)
     k=ncol(X)-1                        #X matrix has the row of ones
     lf <- lsfit(X,y, intercept=FALSE)  #intercept is in the design matrix
     bhat= lf$coef
     R=qr.R(lf$qr)                     #QR decomp for stable matrix inversion
     L=solve(R)
     s2=var(y)
     sig2= (n-k)*s2/rchisq(N,n-k)      #transforms a chi-sq into an Inverse chi-sq
     z = matrix(rnorm((k+1)*N,sd=sqrt(sig2)),byrow=T,ncol=N) 
     b=bhat+L%*%z
     summ=cbind(t(b),sig2)
     summ
}

##### Example
x=seq(0,5, length=100)            #simulate data
y=5+8*x+rnorm(100,0,1)
plot(x,y)
reg1=lm(y~x, x=TRUE)              #use OLS regression to intialize the parameters
X = reg1$x                              #build the X matrix
reg2=lm.bayes(y,X, 100000)