# This file is posted as power.ssc

n _ 2:50  
power _ 1 - pf(qf(.95,2,3*n-3),2,3*n-3,n/2) 

# Plot the power curve

par(mfrow=c(1,1),fin=c(5,5))
plot(c(0,50), c(0,1), type="n", 
      xlab="Sample Size",ylab="Power",
	   main="Power versus Sample Size\n F-test for equal means")
lines(n, power, type="l",lty=1)


#  Create a loop to serch for the 
#  smallest sample size to achieve
#  the specified power level.  Here
#  we create the model matrix for
#  the one way anova needed to 
#  evaluate the noncentrality 
#  parameter

nmax=50
alpha= .05
power= .90

#  Enter the C matrix for the null
#  hypothesis of CB=0

 C <- matrix(c(0, 1, 0, -1, 0, 0, 1, -1),2,4,byrow=T)

#  Enter deviations from the null hypothesis as 
#  multiples of the stadard deviation for the
#  random errors.

d <- matrix(c(0.5, 1), 2,1,byrow=T)

#  Enter the portion of the model matrix
#  that will be replicated as the sample
#  size is increased

X <- cbind(matrix(1,3,1),diag(rep(1,3)))


#  Set the initial sample size
   
n<-2
X<- rbind(X,X)

#  Create a function for computing power

powerfun<-function(n,X,C,d,alpha){
	nc<-t(d)%*%solve(C%*%ginverse(t(X)%*%X)%*%t(C))%*%d
	df1 <- nrow(C)
	df2 <- nrow(X)-qr(X)$rank
   1-pf(qf(.95,df1,df2),df1,df2,nc)
}

#  Enter the loop to evaluate the
#  required sample size

npower<-0
while (n<nmax){
        k <- nrow(X)/n
        n <-n+1
        X <- rbind(X,X[1:k, ])
        if(powerfun(n,X,C,d,alpha)<power) npower<-n+1
}

#  List the sample size.  

npower

# If the value of npower>nmax you have
# not found the desired sample size and you 
# must increase  nmax