#----------------------------------------------------
#永田（2003）第9章 1元配置分散分析―誤差分散が既知の場合
#----------------------------------------------------


sig.level <- 0.05; k <- 3; n <- 5; u <- c(10,11.6,12); sigma <- sqrt(2)^2

#####################################
#1元配置分散分析-誤差分散が既知の場合（検出力）
#水準の平均値を定める場合
#####################################
power.f1 <- function(sig.level,k,n,u,sigma){
delta <- sum((u-sum(u)/length(u))^2)/sigma
lamda <- n * delta
pchisq(q=qchisq(df=k-1,p=sig.level,lower.tail=F),df=k-1,ncp=lamda,lower.tail=F)
}

power.f1(sig.level=0.05,k=3,n=5,u=c(10,11.6,12),sigma=2)


#####################################
#1元配置分散分析-誤差分散が既知の場合（検出力）
#母平均の範囲(d)を定める場合
#####################################
power.f2 <- function(sig.level,k,n,d,sigma){
delta <- (d^2)/(2*sigma)
lamda <- n * delta
pchisq(q=qchisq(df=k-1,p=sig.level,lower.tail=F),df=k-1,ncp=lamda,lower.tail=F)
}


#####################################
#1元配置分散分析-誤差分散が既知の場合（検出力）
#効果量(Δ)を定める場合を定める場合
#####################################
power.f3 <- function(sig.level,k,n,delta){
lamda <- n * delta
pchisq(q=qchisq(df=k-1,p=sig.level,lower.tail=F),df=k-1,ncp=lamda,lower.tail=F)
}





##ex 9.2
power.f1(sig.level=0.05,k=3,n=5,u=c(10,11.6,12),sigma=2)
power.f2(sig.level=0.05,k=3,n=5,d=2,sigma=2)

##ex 9.3
delta <- c(.1,.2,.3,.5,.7,1,1.5,2,2.5,3,4,5,6,7)
power.f3(sig.level=0.05,k=3,n=3,delta=delta)

## k =3の検出力曲線
delta <- seq(0,7,by=.01)
plot(power.f3(sig.level=0.05,k=3,n=3,delta=delta)~delta,type="l")
lines(power.f3(sig.level=0.05,k=3,n=5,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=3,n=7,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=3,n=10,delta=delta)~delta)

## k = 5 の検出力曲線
delta <- seq(0,7,by=.01)
plot(power.f3(sig.level=0.05,k=5,n=3,delta=delta)~delta,type="l")
lines(power.f3(sig.level=0.05,k=5,n=5,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=5,n=7,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=5,n=10,delta=delta)~delta)


## k = 7 の検出力曲線
delta <- seq(0,7,by=.01)
plot(power.f3(sig.level=0.05,k=7,n=3,delta=delta)~delta,type="l",ylim=c(0,1))
lines(power.f3(sig.level=0.05,k=7,n=5,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=7,n=7,delta=delta)~delta)
lines(power.f3(sig.level=0.05,k=7,n=10,delta=delta)~delta)