## Calculate Sizes for Balanced Normal Case

K<-20
ALPHA<-.05
NSIM<-10^5
results<-matrix(NA,8,1)

prob.reject.s<-function(k=K,alpha=ALPHA){
1- pf( ((k-1)/k)* qchisq(1-alpha,1),1,k-1 )
}

### test 1
results[1,1]<- prob.reject.s(K)

## test 2 and 4 by simulation
simulation<-function(nsim,k=K,alpha=ALPHA){
   reject.d2<-rep(0,nsim)
   reject.d4<-rep(0,nsim)
    sqrtk<-sqrt(k)
	for (i in 1:nsim){
	z<-rnorm(k)
	x<-(z-mean(z))^2
   df<-   (sum(x)^2) / sum(x^2)
   Ta2<- k*(mean(z))^2/ ( (1/(k-1))* sum(x) ) 
   Ts2<- (k/(k-1))*Ta2
#print(paste("i=",i))
#print(paste("Ta2=",Ta2))
#print(paste("Ts2=",Ts2))
#print(paste("df=",df))
#print(paste("qf(.95,1,df)=",qf(.95,1,df)))
   if (Ts2 > qf(1-alpha,1,df))   reject.d2[i]<-1
   if (Ta2 > qf(1-alpha,1,df))   reject.d4[i]<-1
}

out<-list(size.d2=mean(reject.d2),size.d4=mean(reject.d4))
out
}

.Random.seed<- c(29,3,18,18,15,1,51,20,35,12,9,3) 

simout<-simulation(nsim=NSIM)
results[2,1]<-simout$size.d2
results[4,1]<-simout$size.d4

### test 3
prob.reject.delta3<-function(k=K,alpha=ALPHA){
1- pf( ((k-1)/k)* qf(1-alpha,1,k-1),1,k-1 )
}
results[3,1]<-prob.reject.delta3()

### test 5
prob.reject.delta5<-function(k=K,alpha=ALPHA){
1- pf( qf(1-alpha,1,k-1),1,k-1 )
}
results[5,1]<-prob.reject.delta5()

## test 6
prob.reject.md<-function(k=K,alpha=ALPHA){
1- pf( (k/(k-1))* qchisq(1-alpha,1),1,k-1)
}
results[6,1]<-prob.reject.md()

## test 7
prob.reject.md2<-function(k=K,alpha=ALPHA,p=P){
1- pf( (k/(k-1))* qf(.95,1,k-p),1,k-1)
}
results[7,1]<-prob.reject.md2(p=1)
results[8,1]<-prob.reject.md2(p=2)


round(results,3)