##################################################
# Central Limit Theorem using Gapminder Data
##################################################

# Clear workspace
rm(list = ls())

# Install and load package
if (!require("gapminder")) {
  install.packages("gapminder", repos = "https://cloud.r-project.org")
  library(gapminder)
}

# Load data
data(gapminder)

# Extract life expectancy
life_data <- gapminder$lifeExp

# Basic statistics
cat("Number of observations:", length(life_data), "\n")
cat("Population Mean Life Expectancy:", mean(life_data), "\n\n")

# Set seed for reproducibility
set.seed(123)

##################################################
# Function to generate sample means
##################################################

simulate_means <- function(n, reps = 1000) {
  replicate(
    reps,
    mean(sample(life_data, n, replace = TRUE))
  )
}

##################################################
# Sampling Distribution for n = 30
##################################################

sample_means_30 <- simulate_means(30)

hist(
  sample_means_30,
  main = "Sampling Distribution (n = 30)",
  xlab = "Sample Mean",
  col = "lightblue",
  border = "black"
)

##################################################
# Compare Sample Sizes
##################################################

sample_sizes <- c(10, 30, 50, 100)

par(mfrow = c(2, 2))

for (n in sample_sizes) {
  
  sim_data <- simulate_means(n)
  
  hist(
    sim_data,
    main = paste("Sample Size =", n),
    xlab = "Sample Mean",
    col = "lightgreen",
    border = "black"
  )
}

par(mfrow = c(1, 1))

##################################################
# Create Results Table
##################################################

results <- data.frame(
  Sample_Size = sample_sizes,
  Mean = NA,
  SD = NA
)

for (i in seq_along(sample_sizes)) {
  
  sim_data <- simulate_means(sample_sizes[i])
  
  results$Mean[i] <- mean(sim_data)
  results$SD[i] <- sd(sim_data)
}

results$Mean <- round(results$Mean, 3)
results$SD <- round(results$SD, 3)

##################################################
# Display Results
##################################################

cat("\nSummary Statistics\n")
print(results)

##################################################
# End of Program
##################################################
