library(tidyverse)

marginError <- function(x, confidenceLevel = 0.95) {
  n <- length(x)
  sem <- sd(x) / sqrt(n)
  p <- (1 + confidenceLevel) / 2
  critT <- qt(p, df = n - 1)
  marginError <- critT * sem
  marginError
}

simulChol <- tibble(
  Cholesterol = rnorm(600, mean=101.2, sd=7.63),
  Simulation = rep(1:100, each=6) 
)

simulChol

simulCholIntervals <- simulChol %>% 
  group_by(Simulation) %>% 
  summarise(
    MeanChol = mean(Cholesterol), 
    MarginError = marginError(Cholesterol), 
    CILower = MeanChol-MarginError, 
    CIUpper = MeanChol+MarginError, 
    ContainsTrueMean = (CILower<=101.2 & CIUpper>=101.2)
  )

ggplot(simulCholIntervals, aes(x=Simulation, ymin=CILower, ymax=CIUpper, colour=ContainsTrueMean)) + 
    geom_errorbar() + geom_hline(yintercept=101.2)


simulChol <- tibble(
  Cholesterol = rnorm(600, mean=101.2, sd=7.63),
  Simulation = rep(1:100, each=6) 
)

simulCholIntervals <- simulChol %>% 
  group_by(Simulation) %>% 
  summarise(
    MeanChol = mean(Cholesterol), 
    MarginError = marginError(Cholesterol, confidenceLevel = 0.9), 
    CILower = MeanChol-MarginError, 
    CIUpper = MeanChol+MarginError, 
    ContainsTrueMean = (CILower<=101.2 & CIUpper>=101.2)
  )

ggplot(simulCholIntervals, aes(x=Simulation, ymin=CILower, ymax=CIUpper, colour=ContainsTrueMean)) + 
  geom_errorbar() + geom_hline(yintercept=101.2)

runSimulation <- function(
  confidenceLevel = 0.95, 
  sampleSize = 6, 
  numberOfSimulations = 100,
  mean = 0,
  sd = 1) {
  
  simul <- tibble(
    Data = rnorm(sampleSize * numberOfSimulations, mean=mean, sd=sd),
    Simulation = rep(1:numberOfSimulations, each=sampleSize) 
  )
  simulIntervals <- simul %>% 
    group_by(Simulation) %>% 
    summarise(
      Mean = mean(Data), 
      MarginError = marginError(Data, confidenceLevel = confidenceLevel), 
      CILower = Mean-MarginError, 
      CIUpper = Mean+MarginError, 
      ContainsTrueMean = (CILower<=mean & CIUpper>=mean)
    )
  print(ggplot(simulIntervals, aes(x=Simulation, ymin=CILower, ymax=CIUpper, colour=ContainsTrueMean)) + 
    geom_errorbar() + geom_hline(yintercept=mean))
  simulIntervals
}

runSimulation(mean = 101.2, sd=7.63)