## Simple linear regression using UMacs package ##load libraries library(Umacs);library(rv);library(MASS);library(mvtnorm) ##simulate data n=500 b1 = 10 b2 = 6 beta = matrix(c(b1,b2),2,1) sig = 4 x=cbind(1,rnorm(n,seq(1,20,length=n),sd=1)) y=matrix(rnorm(n,x%*%beta,sqrt(sig)),n,1) #priors on beta and beta variance b.prior = as.matrix(ncol=1,nrow=2,c(0,0)) v.prior = solve(diag(1000,2)) s1.prior = .1 #for gamma prior on sigma s2.prior = .1 #set up the initial values and define the dimension of each variable b.init <- function () as.matrix(rnorm(2,0,1)) # b[1] is intercept, b[2] is slope v.init <- function () rnorm(1,0,1)^2 # error #define the conditional probabilities for the gibbs sampler for the betas and sigma b.update <- function (){ #From Clark text (page 192) and Clark workbook (pg 81) if(length(v) == 1) { sx = crossprod(x) * 1/v sy = crossprod(x,y) * 1/v } if(length(v) > 1) { ss = t(x) %*% 1/v sx = ss %*% x sy = ss %*% y } bigv = solve(sx + v.prior) smallv = sy + v.prior %*% b.prior b = t(rmvnorm(1,bigv%*%smallv,bigv)) return(b) } v.update <- function () { sx = crossprod((y - x%*%b)) u1 = s1.prior + 0.5*n u2 = s2.prior + 0.5*sx return(1/rgamma(1,u1,u2)) } s <- Sampler( .title = "Linear Regression", x = x, y = y, b = Gibbs (b.update, b.init), v = Gibbs (v.update, v.init), Trace("b[1]"), Trace("b[2]"), Trace("v") ) s(n.iter=2000) r=s() r=as.rv(r) r paste("beta 1 = ", b0,", beta 2 = ",b1," sigma = ",sig) #compare with original values