library(magrittr)


# Function Declaration ----------------------------------------------------

type_1_errors <- function(n_coeff, alpha, n_samples = 100000) {
  # Resampling approach to Type I Error estimation
  #
  # Args:
  #   n_coeff (int): Number of coefficients in model.
  #   alpha (dbl): Confidence level to accept alternative hypothesis.
  #   n_samples (int): Number of iterations to build sampling distribution.
  #     Defaults to 100,000.

  p <- 1- alpha
  samples <- rbinom(n_samples, n_coeff, p)
  outcome_probs <-samples %>%
    table() %>%
    prop.table()

  expected_value <- sum(as.double(names(outcome_probs)) * outcome_probs)  # weighted average of outcomes

  list(
    n_coeff = n_coeff,
    alpha = alpha,
    outcomes = outcome_probs,
    percent_outcome = outcome_probs * 100,
    expected_errors = expected_value,
    raw = samples
    )
}


# Example -----------------------------------------------------------------

n_coeff <- 63
alpha <- 0.95
error <- type_1_errors(n_coeff, alpha)

print(paste0('model with ', n_coeff, ' coefficients ...'))
print(paste0('expected type I errors with alpha=', error$alpha, ': ', error$expected_errors))
print('error distribution (%):')
print(error$percent_outcome)