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