カイ二乗適合度検定

カイ二乗適合度検定は、観測された度数分布が、理論分布と等しいかどうかを調べるときに利用される。これは、多項分布の出力値から算出された統計量がカイ二乗分布に従うことを根拠としている。
ここでは多項分布に従う乱数を出力し、統計量がカイ二乗分布に従うことをシミュレーションによって確かめてみよう。

1 多項分布とカイ二乗値

nは試行回数、sizeは１回の試行における頻度の総数、probは各クラスの実現確率 （確率は自動的に正規化される） 関数chisq.testの出力からカイ二乗値をrsとして取り出す。
これを1000回繰り返してヒストグラムを描く。
理論上の分布は自由度１のカイ二乗分布

X = c()
s <- rmultinom(n=1000, size = 10, prob = c(0.5,0.5)) 
for (i in 1:1000) {
  res <- chisq.test(s[,i], p = c(0.5,0.5))
  rs <- res$ statistic
  X <-c(X,rs)
}
hist(X,prob=T)
curve(dchisq(x,1),add=T)

今度は関数を用いずカイ二乗値の定義に従って計算する。

s <- rmultinom(n=1000, size = 10, prob = c(0.5,0.5)) 
#ds <- data.frame(X=s[1,],Y=s[2,])
Z=c()
for (i in 1:1000){
  z <- sum((s[,i]-10*c(0.5,0.5))^2/(10*c(0.5,0.5)))
  Z <- c(Z,z)
}
#Z
hist(Z,prob=T)

市原p123 例4-4

X = c()
s <- rmultinom(n=1000, size = 60, prob = c(1,1,1,1,1,1)) 
for (i in 1:1000) {
  res <- chisq.test(s[,i], p = c(1,1,1,1,1,1)/6)
  rs <- res$ statistic
  X <-c(X,rs)
}
hist(X,prob=T)
curve(dchisq(x,5),add=T)

検出力
どの程度比率が変われば、有意となるか？

X = c()
s <- rmultinom(n=1000, size = 60, prob = c(2,0.5,0.5,1,1,1)/6) 
for (i in 1:1000) {
  res <- chisq.test(s[,i], p = c(3,0.5,0.5,0.5,0.5,1)/6)
  rs <- res$ statistic
  X <-c(X,rs)
}
hist(X,prob=T)
curve(dchisq(x,5),add=T)

セルの数が５以下のものがある場合は、Rが警告を出す。
このような場合はフィッシャーの正確検定の考えで直接確率を計算する必要がある。

chisq.testを使わずにやってみるが、分布の当てはまりはそれほど悪くない。

s <- rmultinom(n=1000, size = 6, prob = c(1,1,1,1,1,1)/6) 
#ds <- data.frame(X=s[1,],Y=s[2,])
Z=c()
for (i in 1:1000){
  z <- sum((s[,i]-6*c(1,1,1,1,1,1)/6)^2/(6*c(1,1,1,1,1,1)/6))
  Z <- c(Z,z)
}
#Z
hist(Z,prob=T)
curve(dchisq(x,5),add=T)