############################## ## Simple linear regression using UMacs package ##load libraries library(Umacs) library(rv) library(MASS) library(mvtnorm) ##simulate data n=256 b1 = 30 b2 = 20 beta = matrix(c(b1,b2),2,1) sig = 10 data=read.table('DadRegLiPonto.txt',header=T) X=data$PH # Y=data$SB x=cbind(1,X) 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 = .001 # gamma prior s2.prior = .001 #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 (){ 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) #note original values are (approximately) returned