# 2009 年 12 月 29 日作成

# 有意水準 alpha の尤度比検定統計量（n が自由度）のしきい値の条件
g<-function(x,y,n){
x**(n/2)*exp(-x/2)-y**(n/2)*exp(-y/2)
}

n<-5
c1<-0.1
f<-function(x){g(x,c1,n)}
c2<-uniroot(f,interval=c(n,100))$root

# 有意水準 alpha の尤度比検定統計量（n が自由度）のしきい値の計算
myquantile2<-function(alpha,n){
repp<-20000
res<-0
resx<-0
for (i in 1:repp){
c1<-n-0.001*i
f<-function(x){g(x,c1,n)}
c2<-uniroot(f,interval=c(n,100))$root
p<-1-pchisq(c2,n)+pchisq(c1,n)
res<-append(res,p)
resx<-append(resx,n-0.001*i)
i<-ifelse(p-alpha>= 0,i,break) 
}
res<-res[2:length(res)]
resx<-res[2:length(resx)]
return(c(c1,c2,p))
}

alpha<-0.05
n<-10
# 通常の検定のしきい値
u1<-qchisq(alpha/2,n)
u2<-qchisq(1-alpha/2,n)
# power1 は通常の検定の検出力
# power2 は尤度比検定の検出力
power1<-0
power2<-0
sigma2<-0
for (i in 20:200){
sigma<-0.01*i
sigma2<-append(sigma2,sigma)
power<-pchisq(u1/sigma,n)+1-pchisq(u2/sigma,n)
power1<-append(power1,power)
power<-pchisq(c1/sigma,n)+1-pchisq(c2/sigma,n)
power2<-append(power2,power)
}
power1<-power1[2:length(power1)]
power2<-power2[2:length(power2)]
sigma2<-sigma2[2:length(sigma2)]

# power の比較のグラフの作図
plot(sigma2,power1,type="l",col=1,ylim=c(0,1),ylab="power")
par(new=T)
plot(sigma2,power2,type="l",col=2,ylim=c(0,1),ylab="")
abline(h=alpha)