# Edps 590BAY:  Getting what you paid for 
# Sept 23, 2021
# 
# Gibbs sampling mean and variance
#
# 



mydata <- list( y = df$ave_tot,
                N = nrow(df),
				meanY=mean(df$ave_tot),
				sdY = sd(df$ave_tot)
			)


jagsModel <- "model {
    for (i in 1:N) {
	   y[i] ~ dnorm(mu, precision)
	}
     mu ~ dnorm(meanY, 1/(100*sdY^2))
	 precision ~ dgamma(.01,.01)
	 sigma <- 1/sqrt(precision)
}"

writeLines(jagsModel, con="jagsModel.txt")

################### rjags #########################
meanY <- mean(df$ave_tot)
sdY <- sd(df$ave_tot)
prec <- 1/var(df$ave_tot)
initValues <- list(mu=meanY, precision=prec)

jModel <- jags.model(file="jagsModel.txt",
              data=mydata,
			  inits=initValues,
			  n.chains=4,
			  n.adapt=500)
			

samples1 <- coda.samples(jModel,
                variable.names = c("mu", "precision", "sigma"),
				n.iter=5000
				)
 
plot(samples1)

gelman.diag(samples1)

geweke.diag(samples1)