Kolmogorov-Smirnov Test

# One-Sample Kolmogorov-Smirnov Test
# Test if exam scores follow N(75, 15)

# Load packages
library(ggplot2)
library(gridExtra)
library(boot)

# Set seed for reproducibility
set.seed(42)

# Generate exam score data (n=50)
# True distribution: N(73, 16) - slightly off from N(75, 15)
exam_scores <- rnorm(50, mean = 73, sd = 16)

# ====================
# 1. CHECK ASSUMPTIONS
# ====================

cat("ASSUMPTION CHECKS\n")
cat("==================\n\n")

# A1: Data is continuous
cat("1. Continuous data:\n")
cat("   Unique values:", length(unique(exam_scores)), "out of", length(exam_scores), "\n")
cat("   Data type:", class(exam_scores), "\n\n")

# A2: Independence
cat("2. Independence: Verified by design(different students)\n\n")

# A3: Fully specified distribution
cat("3. Fully specified distribution:\n")
cat("   Testing against N(75, 15) - parameters from test design\n")
cat("   NOT estimated from data(would require Lilliefors test)\n\n")

# A4: Ties check
cat("4. Ties check:\n")
ties_prop <- 1 - length(unique(exam_scores)) / length(exam_scores)
cat("   Proportion of ties:", round(ties_prop * 100, 2), "%\n")
if (ties_prop < 0.05) {
  cat("   ✓ Few ties - KS test appropriate\n\n")
} else {
  cat("   ⚠ Many ties - consider Anderson-Darling\n\n")
}

# ====================
# 2. RUN KS TEST
# ====================

cat("KOLMOGOROV-SMIRNOV TEST\n")
cat("========================\n\n")

# One-sample KS test
ks_result <- ks.test(exam_scores, "pnorm", mean = 75, sd = 15)

cat("Hypothesis Test Results:\n")
cat("  D statistic:", round(ks_result$statistic, 4), "\n")
cat("  p-value:", round(ks_result$p.value, 4), "\n")
cat("  Sample size: n =", length(exam_scores), "\n\n")

# Decision
alpha <- 0.05
if (ks_result$p.value < alpha) {
  cat("Decision: REJECT H₀ (α = 0.05)\n")
  cat("Interpretation: Scores do NOT follow N(75, 15)\n\n")
} else {
  cat("Decision: FAIL TO REJECT H₀ (α = 0.05)\n")
  cat("Interpretation: Insufficient evidence that scores differ from N(75, 15)\n\n")
}

# ====================
# 3. BOOTSTRAP CI FOR D
# ====================

cat("Bootstrap 95% CI for D statistic\n")
cat("=================================\n\n")

boot_ks <- function(data, indices) {
  d <- data[indices]
  ks_test <- ks.test(d, "pnorm", mean = 75, sd = 15)
  return(ks_test$statistic)
}

boot_results <- boot(exam_scores, boot_ks, R = 2000)
boot_ci <- boot.ci(boot_results, type = "perc")

cat("  95% CI for D: [", round(boot_ci$percent[4], 4), ",", 
    round(boot_ci$percent[5], 4), "]\n\n")

# ====================
# 4. EFFECT SIZE
# ====================

cat("EFFECT SIZE\n")
cat("============\n\n")

cat("  D statistic:", round(ks_result$statistic, 4), "\n")
cat("  Interpretation:\n")
if (ks_result$statistic < 0.10) {
  cat("    Small deviation from N(75, 15)\n")
} else if (ks_result$statistic < 0.25) {
  cat("    Moderate deviation from N(75, 15)\n")
} else {
  cat("    Large deviation from N(75, 15)\n")
}

# Cohen's d for location shift
mean_diff <- mean(exam_scores) - 75
cohens_d <- mean_diff / 15
cat("\n  Location shift(Cohen's d):", round(cohens_d, 3), "\n")
cat("    (Difference in means / theoretical SD)\n\n")

# ====================
# 5. VISUALIZATIONS (6 plots)
# ====================

# Plot 1: ECDF vs Theoretical CDF
p1 <- ggplot(data.frame(x = exam_scores), aes(x)) +
  stat_ecdf(geom = "step", color = "blue", size = 1) +
  stat_function(fun = pnorm, args = list(mean = 75, sd = 15), 
                color = "red", linetype = "dashed", size = 1) +
  labs(title = "ECDF vs Theoretical CDF",
       subtitle = "Blue = Empirical, Red = N(75, 15)",
       x = "Exam Score", y = "Cumulative Probability") +
  theme_minimal()

# Plot 2: Histogram with theoretical density
p2 <- ggplot(data.frame(x = exam_scores), aes(x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 15, 
                 fill = "lightblue", color = "black", alpha = 0.7) +
  stat_function(fun = dnorm, args = list(mean = 75, sd = 15), 
                color = "red", size = 1) +
  labs(title = "Histogram with Theoretical Density",
       subtitle = "Red curve = N(75, 15)",
       x = "Exam Score", y = "Density") +
  theme_minimal()

# Plot 3: Q-Q plot
p3 <- ggplot(data.frame(sample = exam_scores), aes(sample = sample)) +
  stat_qq(distribution = qnorm, dparams = list(mean = 75, sd = 15)) +
  stat_qq_line(distribution = qnorm, dparams = list(mean = 75, sd = 15), 
               color = "red") +
  labs(title = "Q-Q Plot vs N(75, 15)",
       x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()

# Plot 4: Deviation plot (ECDF - Theoretical CDF)
ecdf_fun <- ecdf(exam_scores)
x_seq <- seq(min(exam_scores), max(exam_scores), length.out = 200)
deviations <- ecdf_fun(x_seq) - pnorm(x_seq, mean = 75, sd = 15)

p4 <- ggplot(data.frame(x = x_seq, dev = deviations), aes(x, dev)) +
  geom_line(color = "darkgreen", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_hline(yintercept = c(-ks_result$statistic, ks_result$statistic), 
             linetype = "dotted", color = "red") +
  labs(title = "CDF Deviation Plot",
       subtitle = "ECDF - Theoretical CDF(red lines = ±D)",
       x = "Exam Score", y = "Deviation") +
  theme_minimal()

# Plot 5: Bootstrap distribution of D
p5 <- ggplot(data.frame(D = boot_results$t), aes(D)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "black", alpha = 0.7) +
  geom_vline(xintercept = ks_result$statistic, color = "red", 
             linetype = "dashed", size = 1) +
  labs(title = "Bootstrap Distribution of D Statistic",
       subtitle = paste0("Red line = observed D = ", round(ks_result$statistic, 3)),
       x = "D Statistic", y = "Frequency") +
  theme_minimal()

# Plot 6: Descriptive statistics comparison
stats_df <- data.frame(
  Statistic = c("Mean", "SD", "Median", "IQR"),
  Sample = c(mean(exam_scores), sd(exam_scores), 
             median(exam_scores), IQR(exam_scores)),
  Theoretical = c(75, 15, 75, 15 * 1.349)
)

p6 <- ggplot(stats_df, aes(x = Statistic)) +
  geom_point(aes(y = Sample, color = "Sample"), size = 4) +
  geom_point(aes(y = Theoretical, color = "Theoretical"), size = 4) +
  geom_segment(aes(xend = Statistic, y = Sample, yend = Theoretical), 
               linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("Sample" = "blue", "Theoretical" = "red")) +
  labs(title = "Sample vs Theoretical Statistics",
       y = "Value", color = "") +
  theme_minimal() +
  theme(legend.position = "top")

# Display all plots
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)

# ====================
# 6. INTERPRETATION
# ====================

cat("\nINTERPRETATION\n")
cat("===============\n\n")

cat("The one-sample Kolmogorov-Smirnov test compared exam scores\n")
cat("to the theoretical distribution N(75, 15).\n\n")

cat("Key findings:\n")
cat("  • D =", round(ks_result$statistic, 3), "(maximum CDF deviation)\n")
cat("  • p =", round(ks_result$p.value, 3), "\n")
cat("  • Sample mean:", round(mean(exam_scores), 1), "vs theoretical: 75\n")
cat("  • Sample SD:", round(sd(exam_scores), 1), "vs theoretical: 15\n\n")

if (ks_result$p.value < 0.05) {
  cat("Conclusion: Scores significantly differ from N(75, 15).\n")
  cat("The ECDF shows notable deviation from the theoretical CDF.\n")
} else {
  cat("Conclusion: Scores are consistent with N(75, 15).\n")
  cat("No significant evidence of distributional mismatch.\n")
}