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

# Example - Sampling with Replacement
########################################################################################

set.seed(5)   # Same random results as in the book

box_of_balls <- tibble(
  balls = 1:10
)
box_of_balls

sample_with_replace1 <- slice_sample(box_of_balls, n=5, replace=TRUE)
sample_with_replace1
sample_with_replace2 <- slice_sample(box_of_balls, n=5, replace=TRUE)
sample_with_replace2


# Example - Bootstrap Distribution of Means for the n=75 Dragon Sample
########################################################################################

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

set.seed(6)   # Same dragon sample as in the Sampling Theory chapter, and same bootstrap results as in the book

dragon_sample <- dragon_population |>
  slice_sample(n = 75)       # sample of size 75

dragon_sample |>
  summarize(
    sample_mean = mean(wingspan),
    sample_sd =  sd(wingspan),
    n = n()
  )

# 10,000 re-samples WITH replacement
bootstrap_resamples <- dragon_sample |>
  rep_slice_sample(n=75 , reps=10000 ,  replace=TRUE)
bootstrap_resamples

# Calculate a separate mean for each of the 10,000 re-samples
bootstrap_distribution_of_means <- bootstrap_resamples |>
  group_by(replicate) |>
  summarize(resample_mean = mean(wingspan))
bootstrap_distribution_of_means

# Histogram of Bootstrap Distribution for the 10,000 bootstrap_distribution_of_means
bootstrap_distribution_of_means |>
  ggplot(aes(x = resample_mean)) +
  geom_histogram(bins = 30, color = "white") +
  labs(x = "Bootstrap Sample Means")

# Calculate the Bootstrap Statistics
bootstrap_distribution_of_means |>
  summarize(
      bootstrap_distribution_mean = mean(resample_mean),
      bootstrap_distribution_SE = sd(resample_mean)
  )


# Example - Bootstrap Distributions for Different Numbers of Repetitions
########################################################################################

reps <- 500  # Compare for 500, 1000, 10,000, more

dragon_sample |>
  rep_slice_sample(n=75 , reps=reps ,  replace=TRUE) |>
  group_by(replicate) |>
  summarize(resample_mean = mean(wingspan)) |>
  ggplot(aes(x = resample_mean)) +
    geom_histogram(bins = 40, color = "white") +
    labs(title = paste(reps, "Means of Bootstrap Re-Samples"), x = "Bootstrap Sample Means")