# Edps 590BAY
# Spring 2018
# C.J.Anderson
#
# simple function for Gibbs Sampling or y ~ N(mu,sigma)
#
libaray(coda)

# The function
gibbs <- function(S,y,mu0,t02,s02,nu0,seed){ 
  # Set up 
    ybar <- mean(y)
    s2   <- var(y)
    set.seed(seed)              # to be able to repeat results exactly
    X<-matrix(nrow=S,ncol=2)    # for saving samples at each iteration
    X[1,] <- c(ybar,1/s2)

  for(s in 2:S) {
    
  # new theta value from its full conditional
    mu.n<-  ( mu0/t02 + n*ybar*x[2] ) / ( 1/t02 + n*x[2] )
    t2n<- 1/( 1/t02 + n*x[2] )
    x[1] <- rnorm(1,mu.n,sqrt(t2n))
 
  # new sigma^2 value from its full conditional
    nun<- nu0+n
    s2n<- (nu0*s02 + (n-1)*s2 + n*(ybar-x[1])^2 ) /nun  
    x[2]<- rgamma(1, nun/2, nun*s2n/2)
    X[s,]<- x  
    }
 return(X)
}


# Create some data:
mu = 4 
sigma = 2 
n = 50  
y <- rnorm(n,mu,sigma)

# Flat prior
mu0 <- 0
t02 <- 100
s02 <- 1
seed <- 2374
S <- 1000

sim1 <- gibbs(S,y,mu0,t02,s02,nu0,seed)
sim1 <- mcmc(sim1)

par(mfrow=c(2,2))
traceplot(sim1,smooth=F,type='l',xlab='Iterations', ylab="",main="Marginal Trace, Little Gibbs Function")

# Joint trace:
iter.n<-10
plot( sim1[1:iter.n,],type="l",xlim=range(sim1[1:S,1]), ylim=range(sim1[1:S,2]),
       lty=1,col="grey",xlab=expression(theta),ylab=expression(tilde(sigma)^2),
       main="First 10 iterations: precision x theta")
text(  sim1[1:iter.n,1], sim1[1:iter.n,2], c(1:iter.n) )

iter.n<-S
plot( sim1[1:iter.n,],type="l",xlim=range(sim1[1:S,1]), ylim=range(sim1[1:S,2]),
       lty=1,col="gray",xlab=expression(theta),ylab=expression(tilde(sigma)^2),
       main='All iterations: precision x theta')


mu.post <- mean(sim1[226:S,1])
tau.post <- sd(sim1[226:S,1])
var.post <- 1/sim1[,2]
var.post <- mean(var.post)
sd.post <- sqrt(var.post)
post.stat <- c(mu.post,tau.post,var.post, sd.post)
names(post.stat) <- c("Posterior: mean","tau","Sigma","sqrt(Sigma)")
post.stat