############################## ## 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 # gamma prior 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 "Models for Ecological Data" ISBN 0-691-12178-8 (page 192) and Clark workbook (pg 81) # the If() statements allows this function to be used for a full covariance matrix or just a single variance parameter 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 ###################################################### Here is one of my simple MH samplers using a simple linear regression with a Cauchy error term: ####################################################### x <- c(1.808,1.799,1.179,0.574,3.727,0.946,3.719,1.566,3.596,3.253) y <- c(1.816,1.281,-1.382,0.573,3.793,0.935,1.775,1.474,3.679,3.889) fn = function(x,a=0,b=1){ a+b*x } sample.ab <-function(x,y,a,b,s,da,db){ bstar = runif(1,b-db,b+db) astar = runif(1,a-da,a+da) logalpha = sum(dcauchy(y,location=fn(x,astar,bstar),scale=s,log=T) - dcauchy(y,location=fn(x,a,b),scale=s,log=T)) logu = log(runif(1,0,1)) acc = (logu < logalpha) b = acc*bstar + (1-acc)*b a = acc*astar + (1-acc)*a list(b=b,a=a,acc=acc) } samples = function(x,y,a,b,s,ds){ sstar = runif(1,s-ds,s+ds) while(sstar <= 0){ sstar = runif(1,s-ds,s+ds) } logalpha = sum( dcauchy(y,location=fn(x,a,b),scale=sstar,log=T) - dcauchy(y,location=fn(x,a,b),scale=s,log=T)) - log(sstar/s) logu = log(runif(1,0,1)) acc = (logu < logalpha) s = acc*sstar + (1-acc)*s list(s=s,accs=acc) } sample.abs<-function(n=10000,x,y,a=0,b=1,s=2,da=.2,db =.2,ds=1) { accab <- 0 accs <- 0 A = B = S = rep(NaN,n) for(i in 1:n){ z = sampleab(x,y,a,b,s,da,db) q <- samples(x,y,a,b,s,ds) A[i] = a = z$a B[i] = b = z$b S[i]=s=q$s accab = accab + z$acc accs <- accs +q$accs } invisible(list(a=A, b=B, s=S, accab=accab/n,accs=accs/n)) }