## -------------------------------------------
## Purpose: Casella and Berger Boostrapping/multinomial example simulation (pg 14 in book)
## Author: Amber J. Hackstadt
## -------------------------------------------

## ^^^^^^^^^^^^^^
## simulation function
## ^^^^^^^^^^^^^^
sim_funct = function(set,B){
  avg = rep(-1,B)
  n = length(set)
  for(b in 1:B){
    x_temp = sample(set,n,replace = TRUE)
    avg[b] = mean(x_temp)
  }
  pmf_est = table(avg)/length(avg)
  out = list(avg = avg,pmf_est = pmf_est)
}

## ^^^^^^^^^^^^^^
## run simulation when B = 500
## ^^^^^^^^^^^^^^
sim1 = sim_funct(set = c(2,4,9,12),B =500)
hist(sim1$avg)
sim1$pmf_est

## from book P(avg = 4.75) = 3/64 = 0.046875
## from book P(avg = 8) = 9/128 = 0.0703125

## ^^^^^^^^^^^^^^
## run simulation when B = 5000
## ^^^^^^^^^^^^^^
sim1 = sim_funct(set = c(2,4,9,12),B =5000)
hist(sim1$avg)
sim1$pmf


## ^^^^^^^^^^^^^^
## run simulation when B = 50000
## ^^^^^^^^^^^^^^
sim1 = sim_funct(set = c(2,4,9,12),B =50000)
hist(sim1$avg)
sim1$pmf

## ^^^^^^^^^^^^^^
## run simulation when B = 500000
## ^^^^^^^^^^^^^^
sim1 = sim_funct(set = c(2,4,9,12),B =500000)
hist(sim1$avg)
sim1$pmf