[カイ二乗検定モンテカルロでカイ二乗検定の検出力

sphere MoteCarloPowerChisq 多次元球 分割表
  • 予定パッケージ名
  • 関数名
    • MoteCarloPowerChisq
  • タイトル
MoteCarloPowerChisq
  • 説明
指定次元、真値としてカイ二乗統計量を与え、棄却水準(p値もしくは、カイ二乗値)に対応するパワーを算出
非心カイ二乗分布での算出と、モンテカルロ算出との２方法を返す
  • 使用例
MoteCarloPowerChisq(K=15,df=4)
  • ソース
MoteCarloPowerChisq<-function(as=c(0.1,0.01,0.001),K,df,N=10000,p=TRUE){
	if(p){
		as<-qchisq(as,df,lower.tail=FALSE)
	}
	mr<-K
	xs<-RandomSphere(df,n=N) # 球面上の均一ランダム点座標を作る
	C<-xs[,1] # そのコサイン
	S<-sqrt(1-C^2) # そのサイン
	# 非心カイ二乗分布で出すパワーと球面乱数で出すパワーを格納
	p1<-p2<-rep(0,length(a))

for(xx in 1:length(as)){
	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,df)+pchisq(R2^2,df,lower.tail=FALSE))/2
	}else{
		Q2<-(pchisq(R1^2,df,lower.tail=FALSE)+pchisq(R2^2,df,lower.tail=FALSE))/2
	}

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

}
# 一致する
list(p.ncchisq=p1,p.sim=p2)
}
  • Rdファイル
\name{MoteCarloPowerChisq}
\alias{MoteCarloPowerChisq}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
MoteCarloPowerChisq
}
\description{
指定次元、真値としてカイ二乗統計量を与え、棄却水準(p値もしくは、カイ二乗値)に対応するパワーを算出
非心カイ二乗分布での算出と、モンテカルロ算出との２方法を返す
}
\usage{
MoteCarloPowerChisq(K=15,df=4)
}
\arguments{
\item{as}{検定有意水準のp値リスト、もしくは、カイ二乗値リスト}
\item{K}{カイ二乗値の真値}
\item{df}{自由度}
\item{N}{モンテカルロ乱数点の個数}
\item{p}{TRUEはasがp値、FALSEはasがカイ二乗値}
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
\item{p.ncchisq}{非心カイ二乗分布によるパワー}
\item{p.sim}{モンテカルロシミュレーションによるパワー}
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}
\examples{
MoteCarloPowerChisq(K=3,df=4)
#####
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)