#----------------------------------------------------
#永田（2003）第6章 2つの母分散の比の検定
#----------------------------------------------------

################################
#2つの母分散の比の検定（検出力）
#(1) sigma_1  =/ sigma_2
################################
power.sigma1 <- function(delta,sig.level,n1,n2){
	pf(1/(delta^2*qf(sig.level/2,df1=n2-1,df2=n1-1,lower.tail=F)),df1=n1-1,df2=n2-1,lower.tail=T)+
	pf(qf(sig.level/2,df2=n2-1,df1=n1-1,lower.tail=F)/delta^2,df1=n1-1,df2=n2-1,lower.tail=F)
}

power.sigma1(n1=10,n2=8,delta=2,sig.level=.05)

################################
#2つの母分散の比の検定（検出力）
#(2) sigma_1  > sigma_2
################################
power.sigma2 <- function(delta,sig.level,n1,n2){
	pf(qf(sig.level,df2=n2-1,df1=n1-1,lower.tail=F)/delta^2,df1=n1-1,df2=n2-1,lower.tail=F)
}

power.sigma2(n1=10,n2=8,delta=2,sig.level=.05)


################################
#2つの母分散の比の検定（検出力）
#(3) sigma_1  < sigma_2
################################
power.sigma3 <- function(delta,sig.level,n1,n2){
	pf(1/(delta^2*qf(sig.level,df1=n2-1,df2=n1-1,lower.tail=F)),df1=n1-1,df2=n2-1,lower.tail=T)
}

power.sigma3(n1=10,n2=8,delta=.5,sig.level=.05)



################################
#2つの母分散の比の検定（サンプルサイズ）
#(1) sigma_1  =/ sigma_2
################################
sample.sigma <- function(delta,sig.level,power){
	1+((qnorm(sig.level/2,lower.tail=F)-qnorm(power,lower.tail=F))/
	log(delta))^2
	}

sample.sigma(delta=5/2,sig.level=0.05,power=0.9)


################################
#2つの母分散の比の検定（サンプルサイズ）
#(2) sigma_1  > sigma_2
################################
sample.sigma2 <- function(delta,sig.level,power){
	1+((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	log(delta))^2
	}

sample.sigma2(delta=2,sig.level=0.05,power=0.8)
power.sigma2(delta=2,n1=14,n2=14,sig.level=.05)
power.sigma2(delta=2,n1=15,n2=15,sig.level=.05)


################################
#2つの母分散の比の検定（サンプルサイズ）
#(3) sigma_1  < sigma_2
################################
sample.sigma3 <- function(delta,sig.level,power){
	1+((qnorm(sig.level,lower.tail=F)-qnorm(power,lower.tail=F))/
	log(delta))^2
	}

sample.sigma3(delta=1/3,sig.level=0.05,power=0.95)
power.sigma3(delta=1/3,n1=10,n2=10,sig.level=.05)
power.sigma3(delta=1/3,n1=11,n2=11,sig.level=.05)

##ex 6.2
delta <- seq(0,3.5,by=.01)
plot(power.sigma1(n1=10,n2=8,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma1(n1=21,n2=16,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma1(n1=31,n2=25,delta=delta,sig.level=.05)~delta,type="l")

##ex 6.4
delta <- seq(0,3.5,by=.01)
plot(power.sigma2(n1=10,n2=8,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma2(n1=21,n2=16,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma2(n1=31,n2=25,delta=delta,sig.level=.05)~delta,type="l")


##ex 6.6
delta <- seq(0,3.5,by=.01)
plot(power.sigma3(n1=10,n2=8,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma3(n1=21,n2=16,delta=delta,sig.level=.05)~delta,type="l")
lines(power.sigma3(n1=31,n2=25,delta=delta,sig.level=.05)~delta,type="l")


##ex 6.7
sample.sigma(delta=1/2,sig.level=0.05,power=0.9)