####
set.seed(54321)
n      = 100
sigma2 = 2
beta   = 3
true   = c(beta,sigma2)
x      = rnorm(n)
y      = beta*x+rnorm(n,0,sqrt(sigma2))

par(mfrow=c(1,1),mai=0.9*rep(1,4))
plot(x,y)

# Sufficient statistics
sumx2 = sum(x^2)
sumxy = sum(x*y)

# Prior hyperparameters
tau02  = 9
b0     = 1

# MCMC set up - for half-cauchy prior for sigma2
# ----------------------------------------------
burnin   = 1000
M        = 2000
niter    = burnin+M
draws    = matrix(0,niter,2)
sigma2.d = 1
set.seed(31416)
for (iter in 1:niter){
 
  # Sampling beta|sigma2 - Gibbs step
  var = 1/(1/tau02+sumx2/sigma2.d)
  mean = var*(b0/tau02+sumxy/sigma2.d)
  beta.d = rnorm(1,mean,sqrt(var))

  # Sampling sigma2|beta - MH step
  par1 = (n-2)/2
  par2 = sum((y-x*beta.d)^2)/2
  sig2 = 1/rgamma(1,par1,par2)
  prob = min(1,(1+sigma2.d)/(1+sig2))
  if (runif(1)<prob){
     sigma2.d  = sig2
  }
  draws[iter,]  = c(beta.d,sigma2.d)
}

names = c("beta","sigma2")
par(mfrow=c(2,3),mai=c(0.35,0.35,0.5,0.5))
for (i in 1:2){
  ts.plot(draws[,i],xlab="Iterations",main=names[i])
  abline(h=true[i],col=2,lwd=3)
  acf(draws[,i],main="")
  hist(draws[(burnin+1):niter,i],prob=TRUE,main="")
  abline(v=true[i],col=2,lwd=3)
}