####################
#  (I) Compile all these functions which set up the full conditionals 
given 
#  above. 
############################


I = diag(9)    # Create an identity matrix

sampleb0 = function()
{    # b0 is alpha, the constant parameter in the regression model
   z = y - X%*%b
   m = mean(z) 
   S = s2/length(y)
   b0= m+sqrt(S)*rnorm(1) 
   b0  # returns new value for the constant in the regression model
}

sampleb=function()
{ # samples b ~ p(b|...) = N(m,S)   b is slope parameters S is the V in 
the text
#y1= y -b0
z=y-b0  # we remove the constant regression parameter
lf = lsfit(X,z,intercept=F)  # fit regression model without intercept
Bhat=lf$coef     # slope regression coefficients
R=qr.R(lf$qr)   # QR decomposition gives X'X matrix
VB =t(R)%*% R   # V_Beta^-1
S=solve(I/sb2 + VB/s2)
m=S %*% ( (mu/sb2)+(1/s2)*VB %*% Bhat )  
V=chol(0.5*S+t(S)) # i.e. VV'=S
z=rnorm(9)
b=m+V %*% z
b  # returns new value for beta vector
}

samplemu=function()
{ # samples mu ~ p(mu | .....) = N(m,s)  # sample conditional for mu
   m=mean(b)
   S=sb2/9
   mu=m+sqrt(S)*rnorm(1)
   mu   #returns new value for mu
 }
 
 samples2 = function()   # sample conditional for sigma square
 {
 z=y-b0-X%*%b
 sig2=rchisq(1,n) 
 A=mean(z*z)
 s2=(1/sig2)*(n*A)
 s2
 }
 
 samplesb2 = function()  # sample conditional for b's sigma square
 {
 sb2=rchisq(1,k)
 sb2=(1/sb2)*(sum((b-mu)*(b-mu)))
 sb2
 }

###########################
#   (II) Enter the data, scale it  and standardize it 
###########################
stdz= function(x) 
{
return(  (x-mean(x))/sqrt(var(x)) )
}
# read in data y, x
temp= (1/100)*c(rep(1300,6),rep(1200,6),rep(1100,4))
temp=stdz(temp)
ratio=c(7.5,9.0,11.0,13.5,17.0,23.0,5.3,7.5,11.0,13.5,17.0,23.0,5.3,7.5,11.0,17.0)
ratio=stdz(ratio)
zeit= 
100*c(0.0120,0.0120,0.0115,0.0130,0.0135,0.0120,0.0400,0.0380,0.0320,0.0260,0.0340,0.0410,0.0840,0.0980,0.0920,0.0860)
zeit=stdz(zeit)
X=cbind(temp,ratio,zeit,temp*ratio,temp*zeit,ratio*zeit,temp*temp,ratio*ratio,zeit*zeit)
y=c(49,50.2,50.5,48.5,47.5,44.5,28.0,31.5,34.5,35.0,38.0,38.5,15.0,17.0,20.5,19.5)
n=length(y)  # sample size 


############################
#(III)  Initialize parameters  using the standard regression analysis
############################

lf = lsfit(X,y)
b=lf$coef[2:10]
b
b0 =lf$coef[1] # check whether intercept is last or first coeff
b0
s2= (sum(lf$resid^2))/(n-10)
s2
mu=mean(b)
mu
sb2=var(b)
sb2
th0=c(b0,b,mu,s2,sb2)  # th vector
th0
k = length(b)
k
n = length(y) 
n
###########################
#  (IV) Now do the Gibbs with those initial values to start with
############################

TH= th0  # initialize a matrix of simulated par values 
T0 = 1000    # try 1000 iterations to start with and see what happens

for(i in 1:T0) { # try different simulation lengths 
b = sampleb()
b0= sampleb0()
s2=samples2()
mu=samplemu()
sb2=samplesb2()
tht=c(b0,b,mu,s2,sb2) #th^t vector
TH=rbind(TH,tht)    # add current iteration to matrix of simulations
}

#######################
# (5) Plot the trajectories 
####################

par(mfrow=c(5,3))  # put all traces in the same frame

plot(1:(nrow(TH)-1000), TH[1000:10000,1],xlab="ITERATION", ylab="alpha", 
type="l")
for(i in 2:10) { 
plot(1:nrow(TH),TH[,i],xlab="ITERATION",ylab="beta[i]",type="l") 
}

plot(1:nrow(TH), TH[,11],xlab="ITERATION",ylab="MU",type="l")
plot(1:nrow(TH), TH[,12],xlab="ITERATION",ylab="sigma of y",type="l")
plot(1:nrow(TH), TH[,13],xlab="ITERATION",ylab="sigma of theta", type="l")

# look at your traces and make a decision
# if things don't yet look conversed, save your plot in your file,  and 
# type in the console the following: 

 dev.off()   # to be able to restart another plot later 
   
#and rep the last for loop 
# for another 1000,  to see if things get better. 



##########################################
# (6) Second Chain. Now with different starting values 
###########################################
TH2 =NULL  # clear the matrix of simulated values
b=rep(0,9)  #initial values for b
b0 =mean(y)   # initial value for the intercept
s2=5     # initial value for s2
mu=20    # initial value for mu
sb2=50    # initial value for sb2
th2=c(b0,b,mu,s2,sb2)   # row of initial values

#################
# (7) Simulate another 1000 iterations now
#################

TH2 = th2 # clear th again
                   # it's important though that TH "exists", (else the 
                   # rbind(TH, tht) below would not work. 
 
 T0 = 1000
 for(i in 1:T0) { #try different simulation lengths 
   b = sampleb()
b0= sampleb0()
s2=samples2()
mu=samplemu()
sb2=samplesb2()
tht=c(b0,b,mu,s2,sb2) #th^t vector
TH2=rbind(TH2,tht)    # add current iteration to matrix of simulations
}


###############################
# (8)  Write commands to plot the new chain. Repeat the commands 
# I gave you in (5) above for the plot of TH, but chain what is needed.
#  Remember to close your par(mfrow) window with dev.off()
###############################################







 
#############################
#  8.-Plot both chains together. Run these plots one by one and 
# copy paste. 
#############################

TT = min(nrow(TH), nrow(TH2)) 

matplot(1:TT, 
cbind(TH[1:TT,1],TH2[1:TT,1],xlab="Iteration",ylab="alpha",type="l")
for(i in 2:10){
  matplot(1:TT, cbind(TH[1:TT,i],TH2[1:TT,i],xlab="Iteration", 
ylab="BETA",type="l") 
      # note beta1 is the second parameter in the th vector 
  }

matplot(1:TT, cbind(TH[1:TT, 11], TH2[1:TT,11], 
xlab="ITERATION",ylab="MU",type="l")
matplot(1:TT, cbind(TH[1:TT, 12], TH2[1:TT,12], 
xlab="ITERATION",ylab="sigma of Y",type="l")
matplot(1:TT, cbind(TH[1:TT, 13], TH2[1:TT,13], 
xlab="ITERATION",ylab="sigma of beta",type="l")

# if things don't yet look converged, repeat the last for loop for another 
100 or so iterations. 



############## NOTE: FROM NOW ON, YOU MUST CHANGE THE RANGES 
####### EACH OF YOU MAY HAVE RUN THE CHAIN FOR LONGER OR SHORTER 
######## NUMBER OF ITERATIONS. CHANGE CODE. 

###################################################
#  Box plots of coefficients. This will summarize the marginal 
# posterior distributions p(betaj | y). Plot only when you think things 
converged
# NOTE: YOU MUST CHANGE THE RANGE 500:1000 TO WHAT YOU THINK IS 
# NEEDED
###################################################

boxplot(TH2[500:1000,1],TH2[500:1000,2],TH2[500:1000,3],TH2[500:1000,4],TH2[500:1000,5],TH2[500:1000,6],TH2[500:1000,7],TH2[500:1000,8],TH2[500:1000,9])


##############################################
#  Make a scatter plot of the simulated values for (beta1, beta2). This 
will summarize 
# the joint posterior distribution p(beta1, beta2 | y) 
###############################################

plot(TH2[500:1000,2], 
TH2[500:1000,3],xlab='beta1',ylab='beta2',main='simulated beta1 and 
beta2')

###### You can do summaries

summary(TH2[500:1000,1])   # etc...


######################################
#  95\% intervals of the last 500 simulations for coefficients
######################################

blow=rep(0,10)
bhigh=rep(0,10)

for(i in 1:10){
  blow[i] =quantile(TH2[500:1000], 0.25) 
  bhigh[i]= quantile(TH2[500:1000],0.925)
  }
  
  ####################################
#  Predictive for new observation at x1= 1150, x2=10, x3 = 0.09
# THERE COULD BE SOME CODE TYPOS HERE. FIX THEM. 
######################################

x1=((1/100)*1150 -12.125)/0.80622
x2=(10.0-12.44)/5.66204
x3=(9.0-4.03125)/3.163852
xtilde =c(1,x1,x2,x3,x1*x2,x1*x3,x2*x3,x1*x1,x2*x2,x3*x3)
ytilde=rep(0,500)  #initialize vector of y tilde values
for( i in 1:500){
beta[i]=TH2[500+i,1:10] 
s2 =TH2[500+i,12]
ytilde[i] =rnorm(1, m=t(xtilde)%*%betat, sd=sqrt(s2))
}
hist(ytilde,col=0,xlab="Y[n+1]") # plot the histogram

quantile(ytilde,c(0.025,0.975))






COPY PAST ALL YOUR CODE.