rm(list=ls())
zeroes <- rep.int(0,2000)
ones   <- rep(1,2000)

p <- sum(ones) / length(c(ones,zeroes))
p
## [1] 0.5
pop <- c(zeroes, ones)


## Take a random sample phat with sample size 50 and calculate a test statistic ##

samp_size <- 75

phat <- sum(sample(pop,samp_size)) / samp_size


## Calculate a test statistic.

#z <- (phat - .5 - .01) / sqrt(.5*.5/samp_size)

#z


### Now do that procedure a bunch of times ###


results <- vector(mode = "numeric", length = 1000)


for (i in 1:1000) {
  
  phat <- sum(sample(pop,samp_size)) / samp_size
  z <- (phat - .5)/ sqrt(.5*.5/samp_size)
  results[i] <- z
  rm(phat)
  rm(z)
  }


#results

mistakes <-results[results < -1.96 | results > 1.96]

mistakes
##  [1] -2.193931 -2.655811  2.193931  2.655811  3.117691  1.962991  2.886751
##  [8]  2.655811  2.424871 -1.962991  1.962991 -2.193931  2.424871  1.962991
## [15] -2.193931  2.193931 -2.193931  1.962991  2.193931  1.962991 -2.193931
## [22] -2.193931 -3.348632 -1.962991 -2.193931 -2.655811 -1.962991 -1.962991
## [29]  2.193931  1.962991 -1.962991 -2.193931  2.655811 -2.193931  2.193931
## [36]  1.962991 -1.962991 -1.962991 -1.962991 -1.962991  1.962991 -2.655811
## [43]  2.655811 -2.424871  2.424871  2.193931  1.962991 -2.424871  1.962991
## [50] -2.886751 -1.962991
length(mistakes) / 1000
## [1] 0.051