対立仮説からnoncentral parameterを推定

対立仮説があると、それに合致する分割表が作れる。その分割表から算出されるカイ自乗統計量は、非心カイ自乗分布の最も高確率で観測される値である。
他方、非心カイ自乗分布は、自由度と平均-自由度、をパラメタとする分布である。最も高確率で観測される値が分かっても、そのような対立仮説のときに観測される統計量の平均は分からない。最大確率値において、確率密度関数が極大になることを利用して、以下のように、最大確率カイ自乗値xと自由度dfから、一定の精度accur,accur2とでを推定する関数を作ってみた。

estimateLambda<-function(df,x,accur,accur2){
 
ret<-0
ncp<-c(0.0000001,x+df,100*(x+df))
trio<-c(x-accur2,x,x+accur2)
ans<-matrix(c(0,0,0,0,0,0,0,0,0),ncol=3)
loop<-1
while(loop==1){
 for(i in 1:3){
  for(j in 1:3){
   ans[i,j]=dchisq(trio[i],df,ncp[j])
  }
 }
 delta<-c(ans[1,1]-ans[3,1],ans[1,2]-ans[3,2],ans[1,3]-ans[3,3])

 if(delta[1]*delta[2]<0){
  ncp<-c(ncp[1],(ncp[1]+ncp[2])/2,ncp[2])


 }else{
  ncp<-c(ncp[2],(ncp[2]+ncp[3])/2,ncp[3])
 }
 if(ncp[3]-ncp[1]<accur){loop=0}
}
ret=ncp[2]
return(ret)

}

実行例は

accur<-0.0001
accur2<-0.0001
x<-14
df<-3
est<-estimateLambda(df=df,x=x,accur=accur,accur2=accur2)
k=seq(0,x*2,length=700)

plot(k,dchisq(x=k,df=df,ncp=est),type="l",xlab="",ylab="",ylim=c(0,0.3),cex.axis=1.5)

結果が掲載図。分布の極大がx=14になっていることが分かる。
非心カイ自乗分布に関連する群大青木先生の問答はこちら