# ============================================================================
#  04-07  Sampling & CLT Simulator
#  Chapters 08-01, 08-02
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

set.seed(2026)
pop <- read.csv("data/population.csv", stringsAsFactors = TRUE)
N   <- nrow(pop)
v   <- pop$value

# ── 1. TRUE POPULATION PARAMETERS ───────────────────────────────────────────
cat("\n===== 1. POPULATION PARAMETERS (known exactly) =====\n")

mu    <- mean(v)
sigma <- sqrt(sum((v - mu)^2) / N)        # POPULATION sd: divide by N, not n-1

cat(sprintf("N            %d\n", N))
cat(sprintf("mu           %.4f\n", mu))
cat(sprintf("sigma        %.4f    (population formula)\n", sigma))
cat(sprintf("sample sd    %.4f    (n-1 version, for contrast)\n", sd(v)))
cat(sprintf("skewness     %+.4f   -> the population is right-skewed\n",
            sum((v - mu)^3) / N / sigma^3))

cat("\nBy region:\n")
print(round(cbind(N_i  = table(pop$region),
                  mean = tapply(v, pop$region, mean),
                  sd   = tapply(v, pop$region, sd)), 3))

# ── 2. FOUR DESIGNS, ONE SAMPLE EACH (n = 40) ───────────────────────────────
cat("\n===== 2. FOUR SAMPLING DESIGNS, n = 40 =====\n")

n <- 40

srs <- pop[sample(N, n), ]

k     <- floor(N / n)
start <- sample(k, 1)
sys   <- pop[seq(start, N, by = k)[1:n], ]

alloc <- round(n * table(pop$region) / N)
strat <- do.call(rbind, lapply(names(alloc), function(r) {
  rows <- which(pop$region == r)
  pop[sample(rows, alloc[[r]]), ]
}))

pop$cluster <- rep(1:50, each = 20)                 # 50 clusters of 20
chosen <- sample(50, 2)
clus   <- pop[pop$cluster %in% chosen, ]

cat(sprintf("proportional allocation: %s  (sums to %d)\n",
            paste(names(alloc), alloc, sep = "=", collapse = ", "), sum(alloc)))

one_shot <- data.frame(
  design = c("SRS", "Systematic", "Stratified", "Cluster"),
  n      = c(nrow(srs), nrow(sys), nrow(strat), nrow(clus)),
  mean   = round(c(mean(srs$value), mean(sys$value),
                   mean(strat$value), mean(clus$value)), 3),
  error  = round(abs(c(mean(srs$value), mean(sys$value),
                       mean(strat$value), mean(clus$value)) - mu), 3)
)
print(one_shot, row.names = FALSE)
cat("\nOne sample proves nothing -- repeat each design 2,000 times.\n")

# ── 3. PRECISION COMPARISON (2,000 replications) ────────────────────────────
cat("\n===== 3. PRECISION COMPARISON (2,000 replications each) =====\n")

R <- 2000

srs_means <- replicate(R, mean(sample(v, n)))

sys_means <- replicate(R, {
  st <- sample(k, 1)
  mean(v[seq(st, N, by = k)[1:n]])
})

idx_by_region <- split(seq_len(N), pop$region)
strat_means <- replicate(R, {
  picked <- unlist(lapply(names(alloc), function(r)
    sample(idx_by_region[[r]], alloc[[r]])))
  mean(v[picked])
})

clus_means <- replicate(R, {
  ch <- sample(50, 2)
  mean(v[pop$cluster %in% ch])
})

prec <- data.frame(
  design    = c("SRS", "Systematic", "Stratified", "Cluster"),
  mean_of_means = round(c(mean(srs_means), mean(sys_means),
                          mean(strat_means), mean(clus_means)), 3),
  bias      = round(c(mean(srs_means), mean(sys_means),
                      mean(strat_means), mean(clus_means)) - mu, 4),
  std_error = round(c(sd(srs_means), sd(sys_means),
                      sd(strat_means), sd(clus_means)), 4)
)
prec$relative_to_SRS <- round(prec$std_error / prec$std_error[1], 3)
print(prec, row.names = FALSE)

cat(sprintf("\ntheoretical SRS standard error  sigma/sqrt(n) = %.4f\n", sigma / sqrt(n)))
cat("Every design is UNBIASED (bias ~ 0) but their PRECISION differs.\n")
cat("Stratified wins because the regional means genuinely differ; cluster\n")
cat("loses because each draw commits to whole regions at a time.\n")

# ── 4. CENTRAL LIMIT THEOREM ────────────────────────────────────────────────
cat("\n===== 4. CENTRAL LIMIT THEOREM =====\n")

draw_means <- function(m, reps = 3000) replicate(reps, mean(sample(v, m)))
sizes <- c(1, 5, 30, 100)
clt <- lapply(sizes, draw_means)
names(clt) <- paste0("n=", sizes)

clt_tab <- data.frame(
  n          = sizes,
  sim_mean   = round(sapply(clt, mean), 4),
  sim_se     = round(sapply(clt, sd), 4),
  theory_se  = round(sigma / sqrt(sizes), 4),
  sim_skew   = round(sapply(clt, function(z)
                 sum((z - mean(z))^3) / length(z) / sd(z)^3), 4)
)
print(clt_tab, row.names = FALSE)

cat("\nThree things to notice:\n")
cat(" 1. sim_mean is mu at EVERY n -- the sample mean is unbiased.\n")
cat(" 2. sim_se tracks sigma/sqrt(n) almost exactly.\n")
cat(" 3. sim_skew falls toward 0 as n grows -- that IS the CLT.\n")

# ── 5. BIAS: A CONVENIENCE SAMPLE ───────────────────────────────────────────
cat("\n===== 5. WHY A BIGGER SAMPLE CANNOT FIX BIAS =====\n")

conv40  <- mean(head(sort(v, decreasing = TRUE), 40))
conv400 <- mean(head(sort(v, decreasing = TRUE), 400))
cat(sprintf("true mu                        %.3f\n", mu))
cat(sprintf("convenience sample of 40       %.3f   (bias %+.3f)\n", conv40, conv40 - mu))
cat(sprintf("convenience sample of 400      %.3f   (bias %+.3f)\n", conv400, conv400 - mu))
cat("The bigger convenience sample is still wrong -- just more precisely wrong.\n")

# ── 6. PLOTS ────────────────────────────────────────────────────────────────
png("sampling_plots.png", width = 1200, height = 900)
par(mfrow = c(2, 3), mar = c(4.5, 4.5, 3, 1))

hist(v, breaks = 40, col = "#8A5FBF", border = NA,
     main = "Population (right-skewed)", xlab = "value")
abline(v = mu, col = "#0FA3A3", lwd = 2)

for (i in seq_along(sizes)) {
  hist(clt[[i]], breaks = 40, col = "#0FA3A3", border = NA,
       main = sprintf("Sampling distribution, n = %d", sizes[i]),
       xlab = "sample mean", xlim = range(clt[[1]]))
  abline(v = mu, col = "#5B2A86", lwd = 2)
}

boxplot(list(SRS = srs_means, Systematic = sys_means,
             Stratified = strat_means, Cluster = clus_means),
        col = c("#5B2A86", "#8A5FBF", "#0FA3A3", "#0B7A7A"), las = 2,
        main = "Precision by design (2,000 reps)", ylab = "sample mean")
abline(h = mu, lty = 2, col = "grey40")

par(mfrow = c(1, 1))
dev.off()
cat("\nWrote sampling_plots.png\n")