################################################################################ #Monte Carlo simulation to evaluate the true confidence levels for sigma^2 # # C.I.s # ################################################################################ options(width = 60) #Settings mu<-2.713333 sigma<-2.195581 sigma^2 n<-9 R<-500 alpha<-0.05 #Simulate data all at once set.seed(9811) y.sim<-matrix(data = rnorm(n = n*R, mean = mu, sd = sigma), nrow = R, ncol = n) y.sim[1,] #Sample # var(y.sim[1,]) #t_1 ################################################################################ # PP with the foreach package library(boot) 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) } 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]) c(normal.based, asym, basic, percentile, bca, student) } .libPaths() #Verifes the location of my installed packages is correct #I included export R_LIBS=/lustre/work/stattools/bilder/packages in my # .bash.profile file. Now, I no longer need to use the .libpaths() code below #.libPaths(new = "/work/stattools/bilder/packages") library(foreach) library(doParallel) library(foreach) library(doParallel) cl<-makeCluster(spec = 10) #Don't need to load parallel package because doParallel does it registerDoParallel(cl = cl) start.time<-proc.time() clusterSetRNGStream(cl = cl, iseed = 9182) #Multiple streams of seeds save.all2<-foreach(i = 1:500, .combine = rbind, .packages = "boot") %dopar% { sim.func(y = y.sim[i,], alpha = 0.05) } stopCluster(cl) end.time<-proc.time() save.time<-end.time-start.time cat("\n Number of minutes running:", save.time[3]/60, "\n \n") head(save.all2) ################################################################################ # 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.all2[,1:2], sigma.sq = sigma^2) asym <- summarize(low.up = save.all2[,3:4], sigma.sq = sigma^2) basic <- summarize(low.up = save.all2[,5:6], sigma.sq = sigma^2) percentile <- summarize(low.up = save.all2[,7:8], sigma.sq = sigma^2) bca <- summarize(low.up = save.all2[,9:10], sigma.sq = sigma^2) student <- summarize(low.up = save.all2[,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)) #