######################################################################
# STAT 7630 Bayesian Statistics 
# Peng Zeng @ Auburn University
# 09-02-2026 
######################################################################

library(rstan)

midge = c(1.64, 1.70, 1.72, 1.74, 1.82, 1.82, 1.82, 1.90, 2.08)
midge_data = list(n = length(midge), y = midge)

midge_fit = stan(
    file = "midge.stan",    # Stan program
    data = midge_data,      # named list of data
    chains = 4,             # # of Markov chains
    warmup = 1000,          # # of warmup iterations per chain
    iter = 2000,            # # of iterations per chain
    cores = 1,              # number of cores per chain
    refresh = 0             # no progress shown
)

# print summary of samples 
print(midge_fit)
print(midge_fit, pars = c("mu", "sigma"), probs = c(0.025, 0.5, 0.975));

# plot posterior mean with credible intervals 
plot(midge_fit)

# trace plots 
traceplot(midge_fit, pars = c("mu", "sigma"))

# extract all samples for parameters and log-density 
draws = extract(midge_fit) 
hist(draws$mu)
summary(draws$mu)

# sample from posterior predictive distribution 
n_sims = length(draws$lp__);
y_pred = rnorm(n_sims, draws$mu, draws$sigma); 
hist(y_pred)

######################################################################
# THE END
######################################################################