### 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))
ts.plot(alpha.Gibbs,main="alpha")
ts.plot(beta.Gibbs,main="beta")
ts.plot(tau.Gibbs,main="tau")

sigma.Gibbs = 1/sqrt(tau.Gibbs)

ts.plot(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)