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 NHL_BMI <- NHL_BMI[1:N[1]] NBA_BMI <- NBA_BMI[1:N[2]] EPL_BMI <- EPL_BMI[1:N[3]] par(mfrow=c(3,1)) hist(NHL_BMI,breaks=25) hist(NBA_BMI,breaks=25) hist(EPL_BMI,breaks=25) par(mfrow=c(1,1)) #NHL_RAN <- runif(N[1],-0.5,0.5) #NBA_RAN <- runif(N[2],-0.5,0.5) #EPL_RAN <- runif(N[3],-0.5,0.5) sigma <- rep(0,3) (sigma[1] <- sd(NHL_BMI)) (sigma[2] <- sd(NBA_BMI)) (sigma[3] <- sd(EPL_BMI)) mu <- rep(0,3) (mu[1] <- mean(NHL_BMI)) (mu[2] <- mean(NBA_BMI)) (mu[3] <- mean(EPL_BMI)) sampsz <- rep(3,0) sampsz[1] <- 7 sampsz[2] <- 7 sampsz[3] <- 7 sumssz <- sum(sampsz) N.sim <- 100000 set.seed(4721) ybar1 <- numeric(N.sim) sd1 <- numeric(N.sim) ybar2 <- numeric(N.sim) sd2 <- numeric(N.sim) ybar3 <- numeric(N.sim) sd3 <- numeric(N.sim) ybar <- numeric(N.sim) for (i in 1:N.sim) { 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] y2 <- NBA_BMI[samp2] y3 <- EPL_BMI[samp3] ybar1[i] <- mean(y1); sd1[i] <- sd(y1) ybar2[i] <- mean(y2); sd2[i] <- sd(y2) ybar3[i] <- mean(y3); sd3[i] <- sd(y3) ybar[i] <- (mean(y1)+mean(y2)+mean(y3))/3 } SST <- 7*((ybar1-ybar)^2 + (ybar2-ybar)^2 + (ybar3-ybar)^2) SSE <- (7-1)*(sd1^2 + sd2^2 + sd3^2) MST <- SST/(3-1) MSE <- SSE/(sumssz-3) F <- MST/MSE (f.alpha <- qf(.95,3-1,sumssz-3)) sum(F >= f.alpha)/N.sim