# ============================================================================
#  04-02  Frequency Distribution Builder
#  Chapters 02-01, 02-02
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

x <- read.csv("data/service_times.csv")$minutes
n <- length(x)

# ── 1. CLASS DESIGN ─────────────────────────────────────────────────────────
cat("\n===== 1. CLASS DESIGN =====\n")

rng <- diff(range(x))
k_sturges <- ceiling(1 + 3.322 * log10(n))
k_sqrt    <- ceiling(sqrt(n))
w_raw     <- ceiling(rng / k_sturges)
w         <- ceiling(w_raw / 10) * 10          # round up to a clean multiple of 10
start     <- floor(min(x) / w) * w

cat(sprintf("n                 %d\n", n))
cat(sprintf("min / max         %.1f / %.1f\n", min(x), max(x)))
cat(sprintf("range R           %.1f\n", rng))
cat(sprintf("k (Sturges)       %d      k (sqrt n) = %d\n", k_sturges, k_sqrt))
cat(sprintf("raw width         %d  ->  chosen width w = %d\n", w_raw, w))
cat(sprintf("first lower limit %d\n", start))

# boundaries sit half a unit outside the limits; the data is recorded to 0.1,
# so use 0.05 -- but 0.5 is the textbook convention for whole-unit data.
lower <- seq(start, by = w, length.out = ceiling((max(x) - start) / w) + 1)
upper <- lower + w - 1
breaks <- c(lower - 0.5, tail(upper, 1) + 0.5)

# ── 2. THE FREQUENCY TABLE ──────────────────────────────────────────────────
cat("\n===== 2. FREQUENCY DISTRIBUTION =====\n")

f <- as.integer(table(cut(x, breaks = breaks, right = TRUE)))

dist <- data.frame(
  class      = paste0(lower, "-", upper),
  lower_bnd  = head(breaks, -1),
  upper_bnd  = tail(breaks, -1),
  midpoint   = (lower + upper) / 2,
  f          = f,
  rel_f      = round(f / n, 4),
  cum_f      = cumsum(f),
  cum_pct    = round(100 * cumsum(f) / n, 1)
)
print(dist, row.names = FALSE)

cat(sprintf("\nCHECK  sum(f) = %d  (must equal n = %d)\n", sum(f), n))
cat(sprintf("CHECK  sum(rel_f) = %.4f  (must equal 1)\n", sum(dist$rel_f)))

# ── 3. GROUPED vs EXACT STATISTICS ──────────────────────────────────────────
cat("\n===== 3. GROUPED vs EXACT =====\n")

xbar_g <- sum(f * dist$midpoint) / n
s_g    <- sqrt(sum(f * (dist$midpoint - xbar_g)^2) / (n - 1))

cat(sprintf("grouped mean %.3f    exact mean %.3f    error %.2f%%\n",
            xbar_g, mean(x), 100 * abs(xbar_g - mean(x)) / mean(x)))
cat(sprintf("grouped sd   %.3f    exact sd   %.3f    error %.2f%%\n",
            s_g, sd(x), 100 * abs(s_g - sd(x)) / sd(x)))

# ── 4. SHAPE ────────────────────────────────────────────────────────────────
cat("\n===== 4. SHAPE =====\n")

cat(sprintf("mean   %.2f\n", mean(x)))
cat(sprintf("median %.2f\n", median(x)))
cat(sprintf("Pearson SK = 3(mean - median)/s = %.3f\n",
            3 * (mean(x) - median(x)) / sd(x)))
cat("five-number summary: ", paste(round(fivenum(x), 2), collapse = "  "), "\n")
cat("IQR:", round(IQR(x), 2), "\n")

shape <- if (mean(x) > median(x) * 1.02) "RIGHT-SKEWED"
         else if (mean(x) < median(x) * 0.98) "LEFT-SKEWED"
         else "roughly SYMMETRIC"
cat("\nShape:", shape, "\n")
cat("=> report the",
    if (shape == "roughly SYMMETRIC") "MEAN and STANDARD DEVIATION"
    else "MEDIAN and IQR", "\n")

# ── 5. THE FOUR PLOTS ───────────────────────────────────────────────────────
png("frequency_plots.png", width = 1100, height = 850)
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1))

# histogram — bars touch by construction in hist()
hist(x, breaks = breaks, col = "#8A5FBF", border = "white",
     main = "Histogram", xlab = "Service time (minutes)", ylab = "Frequency")

# frequency polygon — MIDPOINTS, closed on the axis
mids <- c(dist$midpoint[1] - w, dist$midpoint, tail(dist$midpoint, 1) + w)
freq <- c(0, f, 0)
plot(mids, freq, type = "b", pch = 19, col = "#5B2A86", lwd = 2,
     main = "Frequency polygon", xlab = "Class midpoint", ylab = "Frequency")
abline(h = 0, col = "grey80")

# ogive — upper BOUNDARIES, starting at 0
plot(breaks, c(0, cumsum(f)), type = "b", pch = 19, col = "#0FA3A3", lwd = 2,
     main = "Ogive (cumulative frequency)",
     xlab = "Upper class boundary", ylab = "Cumulative frequency")
abline(h = n / 2, lty = 2, col = "grey50")
abline(v = median(x), lty = 2, col = "grey50")

# boxplot
boxplot(x, horizontal = TRUE, col = "#0B7A7A",
        main = "Boxplot", xlab = "Service time (minutes)")

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