パワー 非心カイ二乗分布 モンテカルロ

• 昨日、検出力と非心カイ二乗分布のことを書いた
• 中で使用している"RandomSphere()"はこちら

Ntrial<-100
Ma<-30
as<-seq(from=0,to=Ma,length.out=Ntrial) # 検定の有意水準はa^2できることにする
N<-10000 # 乱数点の数指定
k<-sample(2:10,1) # 次元を指定
mr<-runif(1)*5 # その原点からの距離
xs<-RandomSphere(k,n=N) # 球面上の均一ランダム点座標を作る
C<-xs[,1] # そのコサイン
S<-sqrt(1-C^2) # そのサイン
# 非心カイ二乗分布で出すパワーと球面乱数で出すパワーを格納
p1<-p2<-rep(0,Ntrial)

for(xx in 1:Ntrial){
	a<-as[xx]
	tmp<-a^2-mr^2*S^2
	tmp[which(tmp<0)]<-0
	tmp2<-mr*C
	tmp3<-sqrt(tmp)
	tmp4<-tmp2-tmp3
	tmp5<-tmp2+tmp3
	tmp4<-abs(tmp4)
	tmp5<-abs(tmp5)
	R12<-matrix(c(tmp4,tmp5),ncol=2)
	R1<-apply(R12,1,min)
	R2<-apply(R12,1,max)
	plot(R1,R2)
# 検定閾値と真値との大小で処理が異なる
	if(a<mr){
		Q2<-0.5 + (pchisq(R1^2,k)+pchisq(R2^2,k,lower.tail=FALSE))/2
	}else{
		Q2<-(pchisq(R1^2,k,lower.tail=FALSE)+pchisq(R2^2,k,lower.tail=FALSE))/2
	}

p1[xx]<-mean(Q2)
p2[xx]<-pchisq(a^2,k,mr^2,lower.tail=FALSE)

}
# 一致する
plot(p1,p2)