> pdf("E:\\blue_drive\\Rmisc\\graphs\\meniscus.pdf")
Error in pdf("E:\\blue_drive\\Rmisc\\graphs\\meniscus.pdf") : 
  cannot open file 'E:\blue_drive\Rmisc\graphs\meniscus.pdf'
> 
> meniscus <- read.table("http://www.stat.ufl.edu/~winner/data/meniscus.dat",header=F,
+                        col.names=c("method","loadfail","displace","stiff"))
> 
> attach(meniscus)
> 
> set.seed(13579)
> 
> r <- 3
> n_T <- 18
> 
> (ybarall <- mean(loadfail))
[1] 84.95
> (ybar1 <- mean(loadfail[1:6]));  (var1 <- var(loadfail[1:6]))
[1] 102.05
[1] 129.251
> (ybar2 <- mean(loadfail[7:12]));  (var2 <- var(loadfail[7:12]))
[1] 48.88333
[1] 75.44167
> (ybar3 <- mean(loadfail[13:18]));  (var3 <- var(loadfail[13:18]))
[1] 103.9167
[1] 937.6937
> 
> (sstr <- 6*((ybar1-ybarall)^2 + (ybar2-ybarall)^2 + (ybar3-ybarall)^2))
[1] 11717.69
> (sse <- (6-1)*(var1 + var2 + var3))
[1] 5711.932
> 
> (f_obs <- (sstr/(r-1))/(sse/(n_T-r)))
[1] 15.38581
> 
> num_samp <- 10^5-1
> f_rand <- rep(0,num_samp)
> 
> for (i1 in 1:num_samp) {
+ randgrp <- sample(n_T,size=n_T,replace=FALSE)
+ ybar1 <- mean(loadfail[randgrp[1:6]]); var1 <- var(loadfail[randgrp[1:6]])
+ ybar2 <- mean(loadfail[randgrp[7:12]]); var2 <- var(loadfail[randgrp[7:12]])
+ ybar3 <- mean(loadfail[randgrp[13:18]]); var3 <- var(loadfail[randgrp[13:18]])
+ 
+ sstr <- 6*((ybar1-ybarall)^2 + (ybar2-ybarall)^2 + (ybar3-ybarall)^2)
+ sse <- (6-1)*(var1 + var2 + var3)
+ 
+ f_rand[i1] <- (sstr/(r-1))/(sse/(n_T-r))
+ }
> 
> summary(f_rand)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 0.00003  0.28910  0.71170  1.16200  1.51200 30.46000 
> 
> (p_rand <- (sum(f_rand >= f_obs)+1)/(num_samp+1))
[1] 0.00037
> 
> hist(f_rand,xlab="F*",main="Randomization Distribution for F-test")
> abline(v=f_obs)