#### Computing power for the specific sample sizes in Section 16.10


r <- 4
mu <- c(12.5,13,18,21)
n <- c(5,5,4,5)
nT <- sum(n)

(mu_all <- sum(mu*n)/nT)

sigma <- 3.5 

(lambda <- (1/sigma^2)*sum(n*(mu-mu_all)^2))

(power <- 1-pf(qf(1-0.05,r-1,nT-r),r-1,nT-r,lambda))


#### Computing power for increasing (balanced) sample sizes in Section 16.10 

r <- 4
mu <- c(12.5,13,18,21)
sigma <- 3.5
mu_all <- sum(mu)/r

n_vec <- 2:50
power_vec <- rep(0,49)

for (n in 2:50) {
nT <- r*n

lambda <- (1/sigma^2)*sum(n*(mu-mu_all)^2)

power_vec[n-1] <- 1-pf(qf(1-0.05,r-1,nT-r),r-1,nT-r,lambda)

}

print(cbind(n_vec,power_vec))