##
##Example R Code - Running a Simple Gibbs Sampler for mu and phi for normal
##data using informative non-conjugate priors
##

##Set mu and phi (precision = 1/variance)
mu <- 5
sigma2 <- 2
phi <- 1/sigma2

#Draw a Sample of size 53 From a N(mu,sigma2)
x <- rnorm(53,mu,sqrt(sigma2))

#Write the data to a data file
write(x,file="./sampledata.dat",sep=",")

#Read the data file and save it as an object
mydata <- scan("./sampledata.dat",sep=",")

#Prior Values
m <- 0
s2 <- 10
a <- 2
b <- 3
n <- length(mydata)

#Allocate Space to save Gibbs sampling draws
totdraws <- 5000
mudraws <- rep(0,5000) #Initial Value of mu is 0
sigma2draws <- rep(1,5000) #Initial Value of s2 is 1

#Begin the Gibbs sampler
for(i in 2:totdraws){

  #Generate a Draw for mu
  mstar <- ((1/s2)*m + (n/sigma2draws[i-1])*mean(mydata))/((1/s2)+(n/sigma2draws[i-1]))
  s2star <- 1/((1/s2)+(n/sigma2draws[i-1]))
  mudraws[i] <- rnorm(1,mstar,sqrt(s2star))

  #Generate a Draw for Sigma2
  astar <- a + (n/2)
  bstar <- (sum((mydata-mudraws[i])^2)/2) + b
  phitemp <- rgamma(1,astar,bstar)
  sigma2draws[i] <- 1/phitemp

}

#Discard the first 500 Draws as Burn-in
mudraws <- mudraws[501:length(mudraws)]
sigma2draws <- sigma2draws[501:length(mudraws)]

#Partition plot into 4 subplots
par(mfrow=c(2,2))

#Plot Trace Plots
plot(mudraws,type="l")
plot(sigma2draws,type="l")

#Plot Marginal Posterior Densities
plot(density(mudraws),col="red",main=expression(paste("Marginal Posterior Density of ",mu)))
plot(density(sigma2draws),col="blue",main=expression(paste("Marginal Posterior Density of ",sigma^2)))