## 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