## Read data from website and attach data frame and obain variable names rr.mar <- read.csv( "http://www.stat.ufl.edu/~winner/data/rocknroll_marathon_mf2015a.csv") attach(rr.mar); names(rr.mar) mph.f <- mph[Gender=="F"] mph.m <- mph[Gender=="M"] N.f <- length(mph.f) N.m <- length(mph.m) set.seed(23456) num.sim <- 100000 n.sample.f <- 8 n.sample.m <- 12 samp.var.f <- rep(0, num.sim) samp.var.m <- rep(0, num.sim) for (i in 1:num.sim) { sample.f <- sample(1:N.f, n.sample.f, replace=FALSE) samp.var.f[i] <- var(mph.f[sample.f]) sample.m <- sample(1:N.m, n.sample.m, replace=FALSE) samp.var.m[i] <- var(mph.m[sample.m]) } sigma2.f <- (N.f-1)*var(mph.f)/N.f sigma2.m <- (N.m-1)*var(mph.m)/N.m X2.s2.f <- (n.sample.f-1)*samp.var.f/sigma2.f X2.s2.m <- (n.sample.m-1)*samp.var.m/sigma2.m X2.L.f <- qchisq(.025,n.sample.f-1) X2.U.f <- qchisq(.975,n.sample.f-1) X2.L.m <- qchisq(.025,n.sample.m-1) X2.U.m <- qchisq(.975,n.sample.m-1) sum(X2.s2.f >= X2.L.f & X2.s2.f >= X2.U.f)/num.sim sum(X2.s2.m >= X2.L.m & X2.s2.m >= X2.U.m)/num.sim hist(X2.s2.f,breaks=seq(0,ceiling(max(X2.s2.f)),.2),freq=F,main="Females") x <- seq(0,ceiling(max(X2.s2.f)),.01) lines(x,dchisq(x,n.sample.f-1),lwd=3,col="red") hist(X2.s2.m,breaks=seq(0,ceiling(max(X2.s2.m)),.2),freq=F,main="Males") x <- seq(0,ceiling(max(X2.s2.m)),.01) lines(x,dchisq(x,n.sample.m-1),lwd=3,col="blue") F.fm <- (samp.var.f/sigma2.f) / (samp.var.m/sigma2.m) F.L.fm <- qf(.025,n.sample.f-1,n.sample.m-1) F.U.fm <- qf(.975,n.sample.f-1,n.sample.m-1) sum(F.fm >= F.L.fm & F.fm <= F.U.fm)/num.sim hist(F.fm,breaks=seq(0,ceiling(max(F.fm)),.1),freq=F,main="Females/Males") x <- seq(0,ceiling(max(F.fm)),.01) lines(x,df(x,n.sample.f-1,n.sample.m-1),lwd=3,col="green") qf(c(.025,.5,.975),n.sample.f-1,n.sample.m-1) quantile(F.fm,c(.025,.5,.975))