## R-code for power calculation with the F-distribution
## AHW3
## March 5, 2008


## power 
t <- 4
sigma.sq <- 5
alpha <- 0.05
avg.sq.treat.eff <- 1

r <- 5
N <- r*t
lambda <- N*(avg.sq.treat.eff)*(1/sigma.sq) 
f.crit <- qf(1-alpha, t-1, N-t)
1-pf(f.crit, t-1, N-t, ncp=lambda)

## plot

## Under H0
x <- seq(0,8, by=0.01)
plot(x, df(x, t-1, N-t), type="l", lwd=2, col="blue")
abline(v=f.crit, lwd=2, col="red")

## Under H1
lines(x, df(x, t-1, N-t, ncp=lambda), lwd=2, col="dark green")

bet.div.with <- avg.sq.treat.eff/sigma.sq
bet.div.with