人文2012共通問題②参照

適合度検定量を、 \(T=\displaystyle \sum_{i=1}^{k}\dfrac{(X_i-Np_i)^2}{Np_i}\)　とする。
\(N\rightarrow \infty\) のとき
　　　\(T=\displaystyle \sum_{i=1}^{k}\dfrac{(X_i-Np_i)^2}{Np_i}\rightarrow \chi^2(k-1)\)　　（ピアソンの定理）
（自由度は級の数 k - 1 だが、未知母数を含む分布の検定する場合、推定した未知母数の数だけ自由度が減る）

基本については当サイト「2_0_適合度検定」を参照

0.1 カイ二乗値の計算例

多項分布を与える乱数 rmultinom を用いて 50回試行における４水準の出現回数を出力し、カイ二乗値を計算する。

d <- rmultinom(n=1, size = 50, prob = c(1,5,3,1))  # 乱数を発生
t(d)#転置
##      [,1] [,2] [,3] [,4]
## [1,]    5   29   12    4
m <-50*c(1,5,3,1)/10;m
## [1]  5 25 15  5
(d-m)
##      [,1]
## [1,]    0
## [2,]    4
## [3,]   -3
## [4,]   -1
(d-m)^2
##      [,1]
## [1,]    0
## [2,]   16
## [3,]    9
## [4,]    1
(d-m)^2/m
##      [,1]
## [1,] 0.00
## [2,] 0.64
## [3,] 0.60
## [4,] 0.20
sum((d-m)^2/m)
## [1] 1.44

0.2 ピアソンの定理

ピアソンの定理は、一般に \(Np_i\geq 5\) であれば近似が良いと言われている。
このことをシミュレーションを通じて確かめよう。

K=c()
m <-50*c(1,5,3,1)/10#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 50, prob = c(1,5,3,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

不思議なことに、試行数を変えても分布の形は影響を受けない。

K=c()
m <-20*c(1,5,3,1)/10#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 20, prob = c(1,5,3,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

確率を変えてみる.

K=c()
m <-50*c(2,4,3,1)/10#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 50, prob = c(2,4,3,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

やはりヒストグラムに大きな変化は見られない。
多項分布の確率分布に左右されないことが特徴である。

0.3 サイコロ

K=c()
m <-60*c(1,1,1,1,1,1)/6#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 60, prob = c(1,1,1,1,1,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,5),add=T)

カテゴリーのまとめ方に左右されるのだろうか？

K=c()
m <-60*c(2,1,1,2)/6#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 60, prob = c(2,1,1,2))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

K=c()
m <-60*c(3,1,1,1)/6#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 60, prob = c(3,1,1,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

驚くべきことに、カテゴリーのまとめ方に左右されない。
自由度のみに左右される。

0.4 イカサマのサイコロ

K=c()
m <-60*c(1,1,1,1,1,1)/6#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 60, prob = c(1,1,1,1,1,1))
  k <- sum((d-m)^2/m)
  K <- c(K,k)
}
L=c()
m <-60*c(1,1,1,1,1,1)/6#期待度数
for (i in 1:3000){
  d <- rmultinom(n=1, size = 60, prob = c(2,0.5,0.5,1,1,1))
  l <- sum((d-m)^2/m)
  L <- c(L,l)
}
hist(K,prob=T)
hist(L,prob=T,add=T,col="lightblue")
curve(dchisq(x,5),add=T)

sum(L>=11.07)/3000#検出力
## [1] 0.8276667

検出力は８０％を超える。
しかし、検出されないことがあることも念頭においておく必要がある。