

# Clear workspace
rm(list = ls())

# Install package if not already installed
if (!requireNamespace("gapminder", quietly = TRUE)) {
  install.packages("gapminder")
}

# Load package
library(gapminder)

# Verify dataset is available
data("gapminder")

# Extract life expectancy variable
life_data <- gapminder$lifeExp

# Check data
cat("Number of observations:", length(life_data), "\n")
cat("Mean life expectancy:", mean(life_data), "\n")

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

sample_means <- replicate(
  1000,
  mean(sample(life_data, size = 30, replace = TRUE))
)

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

##################################################
# Compare n = 10 and n = 100
##################################################

means_10 <- replicate(
  1000,
  mean(sample(life_data, size = 10, replace = TRUE))
)

means_100 <- replicate(
  1000,
  mean(sample(life_data, size = 100, replace = TRUE))
)

par(mfrow = c(1, 2))

hist(
  means_10,
  main = "Sample Size = 10",
  xlab = "Mean",
  col = "lightgreen"
)

hist(
  means_100,
  main = "Sample Size = 100",
  xlab = "Mean",
  col = "lightpink"
)

##################################################
# Function to Simulate Means
##################################################

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

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

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

par(mfrow = c(2, 2))

for (n in sample_sizes) {
  
  simulated_values <- simulate_means(n)
  
  hist(
    simulated_values,
    main = paste("Sample Mean Distribution n =", n),
    xlab = "Sample Mean",
    col = "lightblue"
  )
}

# Reset plotting window
par(mfrow = c(1, 1))

##################################################
# Summary Statistics
##################################################

results <- data.frame(
  Sample_Size = sample_sizes,
  Mean = sapply(sample_sizes, function(n) mean(simulate_means(n))),
  SD = sapply(sample_sizes, function(n) sd(simulate_means(n)))
)

print(results)



