# ============================================================================
#  04-06  Normal Distribution Explorer
#  Chapters 04-02, 07-02
#
#  Run from the project folder:   Rscript r/analysis.R
# ============================================================================

x  <- read.csv("data/heights_cm.csv")$height_cm
n  <- length(x)
mu <- mean(x)
s  <- sd(x)

# ── 1. TWO-WAY CALCULATOR ───────────────────────────────────────────────────
cat("\n===== 1. NORMAL CALCULATOR =====\n")
cat(sprintf("Using the sample estimates: mean = %.3f, sd = %.3f, n = %d\n\n", mu, s, n))

cat("VALUE -> PROBABILITY\n")
for (v in c(160, 168, 175, 182, 190)) {
  cat(sprintf("  x = %3d   z = %+6.3f   P(X < x) = %.4f   P(X > x) = %.4f\n",
              v, (v - mu) / s, pnorm(v, mu, s), pnorm(v, mu, s, lower.tail = FALSE)))
}
cat(sprintf("\n  P(168 < X < 182) = %.4f\n",
            pnorm(182, mu, s) - pnorm(168, mu, s)))

cat("\nPROBABILITY -> VALUE  (the inverse direction)\n")
for (p in c(0.01, 0.05, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99)) {
  cat(sprintf("  P%-3.0f = %7.2f cm   (critical z = %+6.4f)\n",
              100 * p, qnorm(p, mu, s), qnorm(p)))
}
cat(sprintf("\n  middle 95%%: %.2f to %.2f cm\n",
            qnorm(0.025, mu, s), qnorm(0.975, mu, s)))

# ── 2. Z-SCORES ─────────────────────────────────────────────────────────────
cat("\n===== 2. Z-SCORES =====\n")

z <- as.numeric(scale(x))
cat(sprintf("CHECK  mean(z) = %.10f   sd(z) = %.10f   (must be 0 and 1)\n",
            mean(z), sd(z)))
cat(sprintf("|z| > 2 (unusual):      %d of %d  (%.1f%%)\n",
            sum(abs(z) > 2), n, 100 * mean(abs(z) > 2)))
cat(sprintf("|z| > 3 (very unusual): %d of %d  (%.1f%%)\n",
            sum(abs(z) > 3), n, 100 * mean(abs(z) > 3)))

top <- order(-abs(z))[1:5]
cat("\nFive most unusual observations:\n")
print(data.frame(height_cm = x[top], z = round(z[top], 3)), row.names = FALSE)

# ── 3. EMPIRICAL RULE CHECK ─────────────────────────────────────────────────
cat("\n===== 3. EMPIRICAL RULE CHECK =====\n")

emp <- data.frame(
  k          = 1:3,
  lower      = round(mu - 1:3 * s, 2),
  upper      = round(mu + 1:3 * s, 2),
  observed   = round(sapply(1:3, function(k) mean(abs(z) <= k)), 4),
  empirical  = c(0.6827, 0.9545, 0.9973),
  chebyshev  = round(c(NA, 1 - 1/4, 1 - 1/9), 4)
)
print(emp, row.names = FALSE)
cat("\nChebyshev is always satisfied but always weaker; the empirical rule is\n",
    "tight here because the data really is bell-shaped.\n")

# ── 4. NORMALITY ASSESSMENT ─────────────────────────────────────────────────
cat("\n===== 4. NORMALITY ASSESSMENT =====\n")

skew <- sum((x - mu)^3) / n / s^3
kurt <- sum((x - mu)^4) / n / s^4 - 3

cat(sprintf("mean %.3f   median %.3f   difference %.3f\n", mu, median(x), mu - median(x)))
cat(sprintf("skewness         %+.4f   (0 if normal)\n", skew))
cat(sprintf("excess kurtosis  %+.4f   (0 if normal)\n", kurt))
cat(sprintf("Pearson SK       %+.4f\n", 3 * (mu - median(x)) / s))

sw <- shapiro.test(x)
cat(sprintf("\nShapiro-Wilk  W = %.5f, p = %.4f\n", sw$statistic, sw$p.value))
cat(if (sw$p.value > 0.05)
      "=> FAIL TO REJECT normality. The normal model is reasonable.\n"
    else
      "=> REJECT normality. Inspect the Q-Q plot before deciding what to do.\n")
cat("\nRemember: with large n these tests reject trivial departures.\n",
    "The Q-Q PLOT is the primary evidence.\n")

# ── 5. NORMAL APPROXIMATION TO THE BINOMIAL ─────────────────────────────────
cat("\n===== 5. NORMAL APPROXIMATION TO THE BINOMIAL =====\n")
cat("n = 200, p = 0.55, find P(X >= 120)\n\n")

nb <- 200; pb <- 0.55
cat(sprintf("np = %.0f >= 5  and  nq = %.0f >= 5   -> approximation is valid\n",
            nb * pb, nb * (1 - pb)))
mu_b <- nb * pb; sd_b <- sqrt(nb * pb * (1 - pb))
cat(sprintf("mu = %.1f   sigma = %.4f\n\n", mu_b, sd_b))

exact  <- pbinom(119, nb, pb, lower.tail = FALSE)
withcc <- pnorm(119.5, mu_b, sd_b, lower.tail = FALSE)
nocc   <- pnorm(120,   mu_b, sd_b, lower.tail = FALSE)

cat(sprintf("EXACT binomial              %.6f\n", exact))
cat(sprintf("normal WITH correction      %.6f   (error %.6f)\n",
            withcc, abs(withcc - exact)))
cat(sprintf("normal WITHOUT correction   %.6f   (error %.6f)\n",
            nocc, abs(nocc - exact)))
cat(sprintf("\nThe correction reduces the error by a factor of %.0f.\n",
            abs(nocc - exact) / abs(withcc - exact)))

# ── 6. PLOTS ────────────────────────────────────────────────────────────────
png("normal_plots.png", width = 1100, height = 850)
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1))

hist(x, breaks = 20, freq = FALSE, col = "#8A5FBF", border = "white",
     main = "Heights with fitted normal curve", xlab = "Height (cm)")
curve(dnorm(xx, mu, s), xname = "xx", add = TRUE, col = "#0FA3A3", lwd = 2)
abline(v = mu, col = "#0B7A7A", lwd = 2, lty = 2)

qqnorm(x, pch = 19, col = "#5B2A86", main = "Normal Q-Q plot")
qqline(x, col = "#0FA3A3", lwd = 2)

boxplot(x, horizontal = TRUE, col = "#0FA3A3",
        main = "Boxplot", xlab = "Height (cm)")

curve(dnorm(xx, mu, s), from = mu - 4*s, to = mu + 4*s, xname = "xx",
      col = "#5B2A86", lwd = 2, ylab = "density",
      main = "Shaded tail: P(X > 185)")
xs <- seq(185, mu + 4*s, length.out = 200)
polygon(c(185, xs, mu + 4*s), c(0, dnorm(xs, mu, s), 0),
        col = "#0FA3A380", border = NA)
text(mu + 2.6*s, dnorm(mu, mu, s) * 0.5,
     sprintf("%.4f", pnorm(185, mu, s, lower.tail = FALSE)), col = "#0B7A7A")

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