### Linear Regression using Gibbs Sampling

###

### The model: Y(i) = alpha + beta*(x(i)-x.bar) + epsilon(i)

### i = 1,...,n

### n = 5

### Y = [1,3,3,3,5]

### x = [1,2,3,4,5]

###

###################################################################################################

# Assign the data and fit the standard regression model

n <- 5

n

Y <- c(1,3,3,3,5)

Y

x <- c(1,2,3,4,5)

x

Y.bar <- mean(Y)

Y.bar

YmYbar <- Y-Y.bar

YmYbar

SSY <- sum(YmYbar^2)

x.bar <- mean(x)

x.bar

xmxbar <- x-x.bar

xmxbar

SSx <- sum(xmxbar^2)

SSx

SSYx <- sum(YmYbar*xmxbar)

SSYx

par(mfrow=c(2,2))

plot(x,Y)

plot(xmxbar,Y)

Y.reg <- lm(Y ~ xmxbar)

Y.reg

Y.fit <- Y.reg$fit

plot(xmxbar,Y)

lines(xmxbar,Y.fit)

###################################################################################################

### Bayesian Regression

# plot priors on alpha, beta, and tau

par(mfrow=c(2,2))

mu.alpha <- 0

sigma.alpha <- 10

alpha <- seq(mu.alpha-3*sigma.alpha,mu.alpha+4*sigma.alpha,0.1)

sigmasq.alpha <- sigma.alpha^2

tau.alpha <- 1/sigmasq.alpha

palpha <- dnorm(alpha,mu.alpha,sigma.alpha)

plot(alpha,palpha,type="l",main="alpha")

mu.beta <- 0

sigma.beta <- 10

beta <- seq(mu.beta-3*sigma.beta,mu.beta+4*sigma.beta,0.1)

sigmasq.beta <- sigma.beta^2

tau.beta <- 1/sigmasq.beta

pbeta <- dnorm(beta,mu.beta,sigma.beta)

plot(beta,pbeta,type="l",main="beta")

a.tau <- 0.1

b.tau <- 0.1

mean.tau <- a.tau/b.tau

sd.tau <- sqrt(a.tau/(b.tau^2))

tau <- seq(0.01,mean.tau+3*sd.tau,0.01)

ptau <- dgamma(tau,a.tau,b.tau)

plot(tau,ptau,type="l",main="tau")

###################################################################################################

### Run Gibbs Sampler

Burn <- 1000

Reps <- 2000

# storage vectors

alpha.Gibbs <- numeric(Reps)

beta.Gibbs <- numeric(Reps)

tau.Gibbs <- numeric(Reps)

# starting values for chains

alpha.in <- 0

beta.in <- 0

tau.in <- 1

alpha.Gibbs[1] <- alpha.in

beta.Gibbs[1] <- beta.in

tau.Gibbs[1] <- tau.in

# the loop

for(h in 2:Reps){

 

# generate an alpha Gibbssample

pm.alpha <- ((tau.Gibbs[(h-1)]*n)/(tau.Gibbs[(h-1)]*n+a.tau))*Y.bar

+(a.tau/(tau.Gibbs[(h-1)]*SSx + a.tau))*mu.alpha # marginal posterior mean

pp.alpha <- tau.Gibbs[(h-1)]*n + a.tau # marginal posterior precision

alpha.Gibbs[h] <- rnorm(1,pm.alpha,sqrt(1/pp.alpha))

 

# generate a beta Gibbs sample

pm.beta <- ((tau.Gibbs[(h-1)]*SSx)/(tau.Gibbs[(h-1)]*SSx + b.tau))*(SSYx/SSx)

+(b.tau/(tau.Gibbs[(h-1)]*SSx+b.tau))*mu.beta

pp.beta <- tau.Gibbs[(h-1)]*SSx + b.tau

beta.Gibbs[h] <- rnorm(1,pm.beta,sqrt(1/pp.beta))

 

# generate a tau Gibbs sample

pa.tau <- (n/2) + a.tau

SSE <- sum((Y - (alpha.Gibbs[h] + beta.Gibbs[h]*xmxbar))^2)

pb.tau <- (SSE/2) + b.tau

tau.Gibbs[h] <- rgamma(1,pa.tau,pb.tau)

}

alpha.Gibbs[1:10]

beta.Gibbs[1:10]

tau.Gibbs[1:10]

 

###################################################################################################

### Posterior Analysis

# Plot chains

par(mfrow=c(4,1))

tsplot(alpha.Gibbs,main="alpha")

tsplot(beta.Gibbs,main="beta")

tsplot(tau.Gibbs,main="tau")

sigma.Gibbs <- 1/sqrt(tau.Gibbs)

tsplot(sigma.Gibbs,main="sigma")

 

# Posterior Estimates

mean(alpha.Gibbs[(Burn+1):Reps])

median(alpha.Gibbs[(Burn+1):Reps])

sqrt(var(alpha.Gibbs[(Burn+1):Reps]))

mean(beta.Gibbs[(Burn+1):Reps])

median(beta.Gibbs[(Burn+1):Reps])

sqrt(var(beta.Gibbs[(Burn+1):Reps]))

mean(tau.Gibbs[(Burn+1):Reps])

median(tau.Gibbs[(Burn+1):Reps])

sqrt(var(tau.Gibbs[(Burn+1):Reps]))

mean(sigma.Gibbs[(Burn+1):Reps])

median(sigma.Gibbs[(Burn+1):Reps])

sqrt(var(sigma.Gibbs[(Burn+1):Reps]))

# Prior/Posterior Plots

par(mfrow=c(2,2))

# alpha

y.density <- density(alpha.Gibbs[(Burn+1):Reps], width = 1)

y.max <- max(y.density$y)

plot(alpha,palpha,type="l",main="alpha",ylim=c(0,y.max))

lines(y.density,type="h")

# beta

y.density <- density(beta.Gibbs[(Burn+1):Reps], width = 1)

y.max <- max(y.density$y)

plot(beta,pbeta,type="l",main="beta",ylim=c(0,y.max))

lines(y.density,type="h")

# tau

y.density <- density(tau.Gibbs[(Burn+1):Reps], width = 1)

y.max <- max(y.density$y)

plot(tau,ptau,type="l",main="tau",ylim=c(0,y.max))

lines(y.density,type="b")

# sigma

y.density <- density(sigma.Gibbs[(Burn+1):Reps], width = 1)

y.max <- max(y.density$y)

plot(y.density,type="b", main="sigma",ylim=c(0,y.max))

# Fit the model, using posterior means

alpha.est <- mean(alpha.Gibbs[(Burn+1):Reps])

beta.est <- mean(beta.Gibbs[(Burn+1):Reps])

Y.est <- alpha.est + beta.est*(x-x.bar)

par(mfrow=c(1,1))

plot(xmxbar,Y)

lines(xmxbar,Y.est)