I didn't rewrite the whole thing. But here's the first part. It uses the `histogram` function from the lattice package.

    population_data <- data.frame(
        uniform = runif(10000, min = -20, max = 20),
        normal = rnorm(10000, mean = 0, sd = 4),
        binomial = rbinom(10000, size = 1, prob = .5),
        beta = rbeta(10000, shape1 = .9, shape2 = .5),
        exponential = rexp(10000, .4),
        chisquare = rchisq(10000, df = 2)
    )
    
    histogram(~ values|ind, stack(population_data),
              layout = c(6, 1),
              scales = list(x = list(relation="free")),
              breaks = NULL)
    
    take_random_sample_mean <- function(data, sample_size) {
        x <- sample(data, sample_size)
        c(mean = mean(x), sd = sqrt(var(x)))
    }
    
    sample_statistics <- replicate(20000, sapply(population_data, take_random_sample_mean, 60))
    
    sample_mean <- as.data.frame(t(sample_statistics["mean", , ]))
    sample_sd <- as.data.frame(t(sample_statistics["sd", , ]))
    
    histogram(sample_mean[["uniform"]])
    histogram(sample_mean[["binomial"]])
    
    histogram(~values|ind, stack(sample_mean), layout = c(6, 1),
              scales = list(x = list(relation="free")),
              breaks = NULL)