#
#
####### Simulating SETE.
#

################################################################################

N <- 4500

library(gb)
gld.params <- c(8,-.043,-.05,-.25)
pop.n <- rrsgld(n = N + 50, gld.params)
pop.n <- pop.n[pop.n >= 1]
pop.n <- round(pop.n)

hist(pop.n, breaks = 80)
summary(pop.n)

detach("package:gb")

N <- length(pop.n)
pop.means <- rnorm(N, 800, 10)
summary(pop.means)
hist(pop.means)

pop.sd <- rnorm(N, 12, 2)
summary(pop.sd)
hist(pop.sd)

pop <- data.frame(pop.means, pop.sd, pop.n)
rm(pop.means, pop.sd, pop.n, gld.params)
ls()

summary(pop)
nrow(pop)

################################################################################

library(lme4)
intervals.range <- as.vector(0)
course.means.prediction.errors <- as.list(0)

iterations <- 12  ######## THIS IS THE NUMBER OF OUTER LOOPS RUN.

for (i in 1:iterations){  #__# <<< Beginning of the i loop. >>> #__#

data.list <- as.list(0)
course.means <- as.vector(0)
cl <- seq(1:N)

for(j in 1:N){
 xbar <- sample(pop[,1], 1)
 stdev <- sample(pop[,2], 1)
 n <- sample(pop[,3], 1)
 course <- rep(cl[j], n)
 scores <- rnorm(n, xbar, stdev)
 data.list[[j]] <- data.frame(course, scores)
 course.means[j] <- mean(scores)
 rm(xbar, stdev, n, course, scores)
}
rm(cl,j)
ls()

####

temp <- rbind(NULL)

for (k in 1:N){
 d <- matrix(unlist(data.list[[k]]), byrow = FALSE, ncol = 2)
 temp <- rbind(temp, d)
 }
rm(k,d)
data.df <- data.frame(temp)
names(data.df) <- c("course", "scores")
rm(temp)

head(data.df)
tail(data.df)
nrow(data.df)
summary(data.df)

ls()

####

model.lmer <- lmer(scores ~ 1 + (1|course), data.df)
summary(model.lmer)

null.posterior <- mcmcsamp(model.lmer, n = 5000, saveb = TRUE)
predicted.means <- null.posterior@ranef[,1] + null.posterior@fixef[,1] # Predicted mean of each course.
course.means.prediction.errors[[i]] <- predicted.means - course.means

HPD.intervals <- HPDinterval(null.posterior, prob = .90)

intervals.range[i] <- mean(abs(HPD.intervals$ranef[,1] - HPD.intervals$ranef[,2]))
intervals.range[i]

rm(model.lmer, null.posterior, course.means, predicted.means, HPD.intervals, data.df)

} #__# <<< End of the i loop. >>> #__#

################################################################################

rm(i, N, iterations)
ls()

mean(course.means.prediction.errors[[1]])

intervals.range
summary(intervals.range)


detach("package:lme4")
detach("package:Matrix")
detach("package:lattice")