library(tidyverse)
share <- .7
population_size <- 10000
sample_size <- 30

binary_var <- c(rep(T, population_size * share),
            rep(F, population_size * (1-share)))

many_sample_means <- function(x, sample_size, n_samples = 10000) {
  map(1:n_samples, ~sample(x, sample_size)) %>%
    map(mean) %>%
    unlist %>%
    table
}


sample_sizes <- c(10, 30, 50, 100)
xlabs <- str_c("sample size: ", sample_sizes)

l <- c("n = 10" = 10, n_30 = 30, n_50 = 50, n_100 = 100) %>%
  map(~many_sample_means(binary_var, .))

par(mfrow = c(2,2))

for(i in seq_along(l)) {
  plot(l[[i]], main = xlabs[i], ylab = "count of realizations", xlab = "mean", xlim = c(0,1))
}

par(mfrow = c(1,1))