bmi.sim <- read.csv("http://www.stat.ufl.edu/~winner/data/nhl_nba_ebl_bmi.csv", header=TRUE) attach(bmi.sim); names(bmi.sim) N <- rep(0,3) N[1] <- 717 N[2] <- 505 N[3] <- 526 sigma <- rep(0,3) (sigma[1] <- sd(NHL_BMI,na.rm=T)) sigma[2] <- sd(NBA_BMI,na.rm=T) (sigma[3] <- sd(EPL_BMI,na.rm=T)) mu <- rep(0,3) (mu[1] <- mean(NHL_BMI,na.rm=T)) mu[2] <- mean(NBA_BMI,na.rm=T) (mu[3] <- mean(EPL_BMI,na.rm=T)) sampsz <- rep(3,0) sampsz[1] <- 25 sampsz[2] <- 0 sampsz[3] <- 25 sumssz <- sum(sampsz) set.seed(4721) ybar1 <- rep(0,100000) sd1 <- rep(0,100000) ybar2 <- rep(0,100000) sd2 <- rep(0,100000) for (i in 1:100000) { samp1 <- sample(1:N[1],sampsz[1],replace=F) # samp2 <- sample(1:N[2],sampsz[2],replace=F) samp3 <- sample(1:N[3],sampsz[3],replace=F) y1 <- NHL_BMI[samp1] y3 <- EPL_BMI[samp3] ybar1[i] <- mean(y1); sd1[i] <- sd(y1) ybar2[i] <- mean(y3); sd2[i] <- sd(y3) } meandiff <- ybar1-ybar2 mean(meandiff) sd(meandiff) hist(meandiff, breaks=50) sp <- sqrt((sd1^2+sd2^2)/2) sediff <- sp*sqrt(1/sampsz[1] + 1/sampsz[3]) t.stat <- meandiff/sediff sum(abs(t.stat) > qt(.975,(sampsz[1]+sampsz[3]-2)))/100000 varratio <- (sd1^2/sd2^2) sum(varratio >= qf(.975,sampsz[1]-1,sampsz[3]-1)) sum(varratio <= qf(.025,sampsz[1]-1,sampsz[3]-1)) varratio_LB <- varratio*qf(.025,sampsz[3]-1,sampsz[1]-1) varratio_UB <- varratio*qf(.975,sampsz[3]-1,sampsz[1]-1) varratio.pop <- sigma[1]^2/sigma[3]^2 sum(varratio_LB >= varratio.pop) sum(varratio_UB <= varratio.pop)