# STAT 540 Lec 08 Monte Carlo experiments


# Demonstrate the law of large numbers
n <- 10000
mu <- 10
sigma <- 2
X <- rnorm(n,mu,sigma)

# demonstrate cumsum() function
x <- c(1,2,3,4)
cumsum(x)

Xbar <- cumsum(X)/c(1:n)

plot(Xbar,type="l")
abline(h = mu , lty=3)


## From a Bernoulli distribution

n <- 10000
p <- 0.4
X <- rbinom(n,size=1, prob = p) # size is the number of Bernoulli trials
Xbar <- cumsum(X)/c(1:n)

plot(Xbar, type = "l")
abline(h = p, lty = 3)


# Monte Carlo simulation/experiment to check the coverage probability
# of the t-based CI for mu when the data come from an 
# exponential (instead of a normal) distribution.

N <- 10000
n <- 30
lambda <- 1/2
alpha <- 0.05

t_val <- qt(1-alpha/2,n-1)
mu <- 1/lambda # expected of X
cvr <- logical(N)

for(i in 1:N){  

  X <- rexp(n,lambda)
  xbar <- mean(X)
  sn <- sd(X)
  lo <- xbar - t_val * sn / sqrt(n)
  up <- xbar + t_val * sn / sqrt(n)
  
  # check if my CI covered the true mean mu
  cvr[i] <- (lo <= mu) & (up >= mu)

}

mean(cvr) # will coerce logical values to numeric


# Monte Carlo experiment to check of the Type I error rate
# of the F-test for equal means across treatments is 
# sensitive to non-normality in the error terms.

K <- 3 # number of treatment groups
n <- 4 # number of subjects in each treatment group
alpha <- 0.05 # set the significance level
trt <- c(1:K) %x% rep(1,n) # make treatment levels
crit_val <- qf(1-alpha,K-1,K*(n-1)) # this is the critical value
N <- 1000

rej <- logical(N) # makes a vector of length N of logical values
for(i in 1:N){
    
  Y <- rchisq(n*K,1) # all response values come from the same distribution
  out <- lm(Y ~ as.factor(trt)) # this is for fitting linear models
  summ_out <- summary(out)
  Fstat <- summ_out$fstatistic[1] # this is the F statistic value
  rej[i] <- Fstat > crit_val # do I reject H_0?

}

typeIrate <- mean(rej)
typeIrate