#MC simulation section main program #Note to Chris: These simulations are the same as for the STAT 950 project #5 except # that I use a lower value of R here due to using v.jack instead of v.L ################################################################################ # Simulate data options(width = 60) #Change width of output in R to help format my notes #Settings mu <- 2.713333 sigma <- 2.195581 sigma^2 n <- 9 num.sim <- 500 #Number of simulations - do not use R here because R is number of resamples in boot section #Simulate data all at once set.seed(9811) y.sim <- matrix(data = rnorm(n = n*num.sim, mean = mu, sd = sigma), nrow = num.sim, ncol = n) y.sim[1,] #Sample #1 var(y.sim[1,]) #t_1 ################################################################################ # Examine true confidence level library(boot) alpha <- 0.05 calc.t2 <- function(data, i) { d2 <- data[i] var(d2) } calc.t <- function(data, i) { d <- data[i] n <- length(d) l.jack <- empinf(data = d, statistic = calc.t2, stype = "i", type = "jack") v.jack <- var.linear(L = l.jack) c(var(d), v.jack) } calc.t(data = y.sim[1,], i = 1:n) #Calculate the confidence intervals for a data set sim.func <- function(y, alpha = 0.05, R = 1999) { n <- length(y) t <- var(y) normal.based <- (n - 1)*t / qchisq(p = c(1-alpha/2, alpha/2), df = n - 1) mu.hat4 <- 1/n*sum((y - mean(y))^4) asym <- t + qnorm(p = c(alpha/2, 1-alpha/2))*sqrt((mu.hat4 - t^2)/n) boot.res <- boot(data = y, statistic = calc.t, R = R, sim = "ordinary") save.int <- boot.ci(boot.out = boot.res, conf = 1-alpha, type = "all") basic <- c(save.int$basic[4], save.int$basic[5]) percentile <- c(save.int$perc[4], save.int$perc[5]) bca <- c(save.int$bca[4], save.int$bca[5]) student <- c(save.int$student[4], save.int$student[5]) #bca.quant <- save.int$bca[2:3] #This was used to determine the cause of warning messages #c(normal.based, asym, basic, percentile, bca, student, bca.quant) c(normal.based, asym, basic, percentile, bca, student) } #Test function set.seed(7127) #need to set seed due to bootstrap sim.func(y = y.sim[1,]) ################################### #Estimate time to complete using 10 simulations start.time<-proc.time() set.seed(7127) save.intervals.temp<-t(apply(X = y.sim[1:10,], MARGIN = 1, FUN = sim.func)) end.time<-proc.time() save.time<-end.time-start.time cat("\n Number of minutes running:", save.time[3]/60, "\n \n") #Estimated number of minutes for 500 simulated data sets save.time[3]/60 * 500/10 ################################### ################################### #All 500 start.time<-proc.time() set.seed(7127) save.intervals<-t(apply(X = y.sim, MARGIN = 1, FUN = sim.func)) end.time<-proc.time() save.time<-end.time-start.time cat("\n Number of minutes running:", save.time[3]/60, "\n \n") ################################### ################################################################################ # Summarize summarize <- function(low.up, sigma.sq) { true.conf <- mean(ifelse(test = sigma.sq > low.up[,1], yes = ifelse(test = sigma.sq < low.up[,2], yes = 1, no = 0), no = 0), na.rm = TRUE) exp.length <- mean(low.up[,2] - low.up[,1], na.rm = TRUE) exclude <- sum(is.na(low.up[,1])|is.na(low.up[,2])) data.frame(true.conf, exp.length, exclude) } normal.based <- summarize(low.up = save.intervals[,1:2], sigma.sq = sigma^2) asym <- summarize(low.up = save.intervals[,3:4], sigma.sq = sigma^2) basic <- summarize(low.up = save.intervals[,5:6], sigma.sq = sigma^2) percentile <- summarize(low.up = save.intervals[,7:8], sigma.sq = sigma^2) bca <- summarize(low.up = save.intervals[,9:10], sigma.sq = sigma^2) student <- summarize(low.up = save.intervals[,11:12], sigma.sq = sigma^2) interval.names <- c("Normal", "Asymptotic", "Basic", "Percentile", "BCa", "Student") data.frame(interval = interval.names, rbind(normal.based, asym, basic, percentile, bca, student)) #