library(tidyverse)
library(infer)   # Contains rep_slice_sample() function

dragon_population <- read_csv("http://www.cknuckles.com/statsprogramming/data/dragon_data.csv")

# Example - Sampling Distribution of Means of 10,000 Dragon Samples of Size 75
########################################################################################

dragon_population |>
  rep_slice_sample(n = 75, reps = 10000)  # 10,000 independent samples of size 75

means_10000_samples <- dragon_population |>
  rep_slice_sample(n = 75, reps = 10000) |>   # 10,000 independent samples each of size 75
  group_by(replicate) |>                      # subset the data into individual samples
  summarize(
    sample_mean = mean(wingspan)              # separate mean for each sample
  )
means_10000_samples

# Frequency Histogram
# ggplot warns that it removed 2 rows: the empty bars at the two ends fall outside the x-axis limits
ggplot(means_10000_samples, mapping = aes(x = sample_mean)) +
  geom_histogram(bins=50, color="white") +
  scale_x_continuous(limits = c(6, 9), breaks = c(6,6.6,7,7.34,8,9))

# Relative Frequency Histogram
# ggplot internally populates count, which after_stat(count / sum(count)) turns into a proportion.
ggplot(means_10000_samples, mapping = aes(x = sample_mean, y = after_stat(count / sum(count)))) +
  geom_histogram(bins=50, color="white") +
  scale_x_continuous(limits = c(6, 9), breaks = c(6,6.6,7,7.2,7.34,7.5,8,9)) +
  labs(y = "Probability")

# Empirical Sampling Distribution mean and standard deviation
means_10000_samples |>
  summarize(
    samp_dist_mean = mean(sample_mean),
    SE = sd(sample_mean)
  )


# Example - Effect of Sample Size on the Sampling Distribution
########################################################################################

sample_size <- 25   # Compare for 25, 75, 150

means_10000_samples <- dragon_population |>
  rep_slice_sample(n = sample_size, reps = 10000) |>
  group_by(replicate) |>
  summarize(
    sample_mean = mean(wingspan)
  )

means_10000_samples |>
  summarize(
    samp_dist_mean = mean(sample_mean),
    SE = sd(sample_mean)
  )

# The same x-axis for every sample size, so the spreads can be compared
ggplot(means_10000_samples, mapping = aes(x = sample_mean, y = after_stat(count / sum(count)))) +
  geom_histogram(bins=50, color="white") +
  scale_x_continuous(limits = c(5, 9.5), breaks = seq(5,9.5,.5)) +
  labs(title = paste("n =", sample_size), y = "Probability")