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