setwd("C:/Users/Lisa/Dropbox/CSUMB stats class")

library("BRugs") 

## some usual steps (like clicking in WinBUGS):
modelCheck("Normalmodel.txt") # check model file
modelData("Norm_tst300data.txt") # read data file
modelCompile(numChains=2) # compile model with 2 chains
modelInits(rep("Norm_tstinits.txt", 2)) # read init data file
modelUpdate(1500) # burn in 
samplesSet(c("mu", "s2")) # parameters to be monitored
modelUpdate(15000) # 15000 more iterations, no thinning ....
samplesStats("*") # the summarized results

## some plots
samplesHistory("*") # plot the chain,
samplesDensity("*") # plot the densities,
samplesBgr("*") # plot the bgr statistics, and
samplesAutoC("*", 1) # plot autocorrelations of 1st chain
samplesAutoC("*", 2) # plot autocorrelations of 2nd chain

## Pulling out posterior samples of the parameters.
# mu is in the first list, variance (s2) is in the second list.
Para_Out1 <- samplesHistory("*", plot = FALSE)
Para_Out10 <- samplesHistory("*", plot = FALSE)
Para_Out20 <- samplesHistory("*", plot = FALSE)
Para_Out50 <- samplesHistory("*", plot = FALSE)
Para_Out100 <- samplesHistory("*", plot = FALSE)
Para_Out300 <- samplesHistory("*", plot = FALSE)

## Sampling from the posterior parameter distribution.
# What do the results look like as we increase sample size?

par(mfcol = c(3,2))
hist(rnorm(n = 15000, mean = rnorm(n = 15000, mean = 0, sd = sqrt(1/0.0001)),
           sd = sqrt(rgamma(n = 15000, shape = 0.001, scale = 0.001))),
     main = "N = 0", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out1$mu), mean = Para_Out1$mu, sd = sqrt(Para_Out1$s2)), main = "N = 1", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out10$mu), mean = Para_Out10$mu, sd = sqrt(Para_Out10$s2)), main = "N = 10", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out20$mu), mean = Para_Out20$mu, sd = sqrt(Para_Out20$s2)), main = "N = 20", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out50$mu), mean = Para_Out50$mu, sd = sqrt(Para_Out50$s2)), main = "N = 50", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out100$mu), mean = Para_Out100$mu, sd = sqrt(Para_Out100$s2)), main = "N = 100", xlab = "Estimated value")
hist(rnorm(n = length(Para_Out300$mu), mean = Para_Out300$mu, sd = sqrt(Para_Out300$s2)), main = "N = 300", xlab = "Estimated value")

tmp <- rnorm(n = length(Para_Out300$mu), mean = Para_Out300$mu, sd = sqrt(Para_Out300$s2))
mean(tmp)
sd(tmp)