SingleDiceEquivSim <- function(m.sims, n.rolls, pv, t1e){
  #Simulate dice rolls
  dr <- rmultinom(n = m.sims, size = n.rolls, prob =pv)
  
  #Calculate mean
  sim.means <- c((t(dr) %*% matrix(1:6, ncol = 1))/n.rolls)
  
  #Calculate lower and upper test statistic values
  sim.t.lower <- sqrt(n.rolls)*(sim.means - 3.15)/1.69
  sim.t.upper <- sqrt(n.rolls)*(sim.means - 3.85)/1.69
  sim.t.lower.reject <- sim.t.lower > qnorm(p = t1e, lower.tail = FALSE)
  sim.t.upper.reject <- sim.t.upper < qnorm(p = t1e, lower.tail = TRUE)
  
  # Only reject if the test statistic is 
  both.reject <- sim.t.lower.reject & sim.t.upper.reject
  reject.percent <- sum(both.reject)/length(both.reject)
  return(reject.percent)
}

pi0 <- c(1/6 - .07, rep(1/6, times = 4), 1/6 + .07)
pi1 <- rep(1/6, 6)
#Type-I Error 
SingleDiceEquivSim(1e5, 200, pi0, .10)
#Type-II Error 
SingleDiceEquivSim(1e5, 200, pi1, .10)



# Probability of having at least one face show up more than 84 or fewer than 69 rolls in 458 rolls
dr <- rmultinom(n = 1e5, size = 458, prob = rep(1/6, times = 6))
min.vals <- apply(dr, 2, min)
max.vals <- apply(dr, 2, max)
sum(max.vals > 84 | min.vals < 69)/length(max.vals)