# Edps 590BAY
# Fall 2019
# c.j.anderson
#
# Getting what paid for example
#
#  I commented out plotting by function to speed things up

# Metropolois:  just the mean
#                 
#
# Metropolois:  Mean and variance
#
#
library("coda")

#
# Read in data
#
sat <- read.table("D:\\Dropbox\\edps 590BAY\\Lectures\\5 MCMC\\getting_what_paid_for_data.txt",header=TRUE)

head(sat)

##########################################
# sample statistics
##########################################

ybar <- mean(sat$ave_tot)
s2   <- var(sat$ave_tot)
s    <- sqrt(s2)
n    <- length(sat$ave_tot)
min.y <- min(sat$ave_tot)
max.y <- max(sat$ave_tot)
sample.stat <- cbind(n,ybar,s2,s,min.y,max.y)
sample.stat <- as.data.frame(sample.stat)
names(sample.stat) <- c("n","ybar","s2","s","min.y","max.y")
sample.stat
hist(sat$ave_total)

###########################################
# Read in metroLogNorm
###########################################

iter.set <- 1000

###########################################
# Some code to work on tuning
###########################################
iter.set <- 500
jump.sd <- 30

# tried jumps
# very bad: .2, 3, 10
# getting better: 20
# Reasonable:  30

chain1 <- metroLogNorm(iter.set,sat$ave_tot,ybar,s,jump.sd)
t.accept <-table(chain1[,2])
t.accept/iter.set * 100
#
# Check trace
par(mfrow=c(2,1))
iter.x <- seq(1:iter.set)
plot(iter.x,chain1[,1],type='l',
         main='Working on Tuning',
         xlab="Iteration",
         ylab="Simulated Value", 
         lwd=2,
         col="blue")
# density		 
hist(chain1[,1])		 
		 
################################################
# Run 4 chains
################################################		 
iter.set <- 10000
jump.sd <- 30
chain1 <- metroLogNorm(iter.set,sat$ave_tot,ybar,s,jump.sd)
chain2 <- metroLogNorm(iter.set,sat$ave_tot,(ybar+s),s,jump.sd)
chain3 <- metroLogNorm(iter.set,sat$ave_tot,(ybar-s),s,jump.sd)
chain4 <- metroLogNorm(iter.set,sat$ave_tot,(ybar+2*s),s,jump.sd)
		 
par(mfrow=c(1,1))
iter.x <- seq(1:iter.set)
plot(iter.x,chain1[,1],type='l',
         main='5000 iterations, jump=30',
         xlab="Iteration",
         ylab="Simulated Mean",
		 ylim=c(910,1000),
         lwd=2,
         col="blue")		 
lines(iter.x,chain2[iter.x,1], type="l", lwd=2, col="red")
lines(iter.x,chain3[iter.x,1], type='l', lwd=2, col="green")
lines(iter.x,chain4[iter.x,1], type='l', lwd=2, col="orange")
legend("bottomright",  title="Staring value",
  legend = c("ybar", "ybar+s", "ybar-s","ybar+2*s"),  
  col = c("blue","red","green","orange"),
  text.col = c("blue","red","green","orange"),  
  lty=c(1,1),  bty = "y",  cex = 1.0)

  
# Histograms
par(mfrow=c(2,2))
hist(chain1[,1],breaks=25,main="Chain 1: ybar")
hist(chain2[,1],breaks=25,main="Chain 1: ybar")
hist(chain3[,1],breaks=25,main="Chain 1: ybar")
hist(chain4[,1],breaks=25,main="Chain 1: ybar")

# To use in lecture:   coda package
#
#  make chains mcmc objects
#
c1 <- mcmc(chain1[,1])
c2 <- mcmc(chain2[,1])
c3 <- mcmc(chain3[,1])
c4 <- mcmc(chain4[,1])

par(mfrow=c(2,2))
traceplot(c1,col="blue",main="chain1",ylim=c(900,1000))
traceplot(c2,col="blue",main="chain2",ylim=c(900,1000))
traceplot(c3,col="blue",main="chain3",ylim=c(900,1000))
traceplot(c4,col="blue",main="chain4",ylim=c(900,1000))

#
# Combine chains
#
post.dist <- mcmc.list(c1,c2,c3,c4)

par(mfrow=c(1,1))
traceplot(post.dist, main="Chains Mixing?",ylim=c(920,1000))

densplot(post.dist)

autocorr.plot(post.dist)
	
gelman.plot(post.dist)

gelman.diag(post.dist)

effectiveSize(post.dist)

geweke.diag(post.dist,frac1=0.1,frac2=0.5)

summary(post.dist,quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975))

###########################################################33
#  Now estimate mean and variance
#
# Read in function metroNorm2
#
#
#########################################################
y <- sat$ave_tot
jump.sd <- 15
iter.set <- 10000
start <- c(ybar,s)

#
# Intial run to make sure things look good
#
start2a <- c(ybar,s)
chain2a <- metroNorm2(start2a,jump.sd,iter.set)

traceplot(chain2a)
densplot(chain2a)

#
# Now do some more chains
#
start2b <- c(ybar+20,s+10)
chain2b <- metroNorm2(start2b,jump.sd,iter.set)
traceplot(chain2b)
densplot(chain2b)

start2c <- c(ybar+10,s-10)
chain2c <- metroNorm2(start2c,jump.sd,iter.set)
traceplot(chain2c)
densplot(chain2c)

start2d <- c(ybar-20,s-20)
chain2d <- metroNorm2(start2d,jump.sd,iter.set)
traceplot(chain2d)
densplot(chain2d)

post.dist2 <- mcmc.list(chain2a,chain2b,chain2c,chain2d)

# do some diagnositics here and then check results