ピアソンの定理

1 はじめに

1.1 数値シミュレーション.

当サイトでは、乱数シミュレーションによる体験的学習を「ひょっこり方式」と呼んでいます。
あとで知ったのですが、数値シミュレーションを活用した統計の学習方法は、
専修大学教授 小杉考司先生の.

　 「数値シミュレーションで読み解く統計のしくみ～Rでためしてわかる心理統計」.

において詳しく紹介されています。 統計ソフトRについても丁寧に説明されており、Rを使って統計を学びたい方にお勧めです。 また、当サイトについてお知らせしたところ、小杉先生から励ましの言葉をいただき、統計学習のモチベーションが高まりました。小杉先生の教材は「分散分析のアルゴリズム」など、とても勉強になりました。この場を借りて御礼申し上げます。

1.2 ピアソンの定理.

ピアソンの定理は、
分割表の適合度検定量　　\(T=\displaystyle \sum_{i=1}^{k}\dfrac{(X_i-Np_i)^2}{Np_i}\)　がカイ二乗分布に分布収束するというものです。
高校の分割表関連では「数学Ⅰでクロス集計」「情報Ⅰでカイ二乗検定」が扱われています。
数学的証明を高校で扱うことはありませんが、実験も証明もなしで使うことには抵抗があります。
そこで、探究学習の一環として、
　　ピアソンの定理を、乱数シミュレーションで体験的に理解する　　という試みを行います。

2 二項分布.

二項分布 \(B(10,0.5)\)に従う二項乱数を100個出力する。
str(hist)をみると、

• breakesは「より大きく、〜以下」であることが分かる.
• 最小区間より小さい値は最小区間に含められている.
  
• plot = FALSE とするとグラフは出力されない.

2.1 乱数によるデータ.

a <- rbinom(100,10,0.5)#二項乱数
#sort(a)
h <- hist(a,breaks=seq(0, 10, by = 1), border=1,xaxp=c(0, 10,1))

str(h)

## List of 6
##  $ breaks  : num [1:11] 0 1 2 3 4 5 6 7 8 9 ...
##  $ counts  : int [1:10] 3 5 7 29 17 20 11 7 1 0
##  $ density : num [1:10] 0.03 0.05 0.07 0.29 0.17 0.2 0.11 0.07 0.01 0
##  $ mids    : num [1:10] 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5
##  $ xname   : chr "a"
##  $ equidist: logi TRUE
##  - attr(*, "class")= chr "histogram"

#sum(h$ counts)

2.2 理論確率の分布.

hist関数の区切り h$breaks に合わせて理論確率の分布を計算する.

• 累積確率の階差として確率分布を求める.
• 累積確率の階差として分布確率を計算しているので、合計は１にならない.
• 累積確率の初項pnom[1]を加えると１になる.

二項分布の累積確率関数は pbinom(q, size, prob) である.

ヒストグラムの区切りで累積確率を計算.

#累積確率
pnom <- pbinom(h$breaks, 10, 0.5)#累積確率
str(pnom)

##  num [1:11] 0.000977 0.010742 0.054687 0.171875 0.376953 ...

差分関数 diff により階差数列を計算.

#density
qnom <-diff(pnom, lag = 1)#差分
str(qnom)

##  num [1:10] 0.00977 0.04395 0.11719 0.20508 0.24609 ...

sum(qnom)

## [1] 0.9990234

c(pnom[1],qnom)

##  [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250
##  [6] 0.2460937500 0.2050781250 0.1171875000 0.0439453125 0.0097656250
## [11] 0.0009765625

sum(pnom[1],qnom)

## [1] 1

初項を付け加えた (pnom[1],qnom) が確率分布を表す。

2.2.1 ヒストグラムの比較.

par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(1,2))
plot(h,xlab="乱数出力")
barplot(100*c(pnom[1],qnom),names.arg=c(seq(0,10,by=1)),xlab="理論分布")#初項を付け加える

#barplot(100*pnom,names.arg=c(seq(0,10,by=1)))#累積度数曲線

2.3 乖離度の計算.

乱数による出現分布と理論確率による分布との乖離度をカイ二乗値によって計算する.

a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(0, 10, by = 1),plot=FALSE)
#累積確率
pnom <- pbinom(h$breaks, 10, 0.5)#累積確率
qnom <-diff(pnom, lag = 1)#差分
str(qnom)

##  num [1:10] 0.00977 0.04395 0.11719 0.20508 0.24609 ...

y <- 100*c(1-sum(qnom[2:10]),qnom[2:10])
str(y)

##  num [1:10] 1.07 4.39 11.72 20.51 24.61 ...

sum(y<5)/10

## [1] 0.5

度数５未満のセルが多くあり、近似の精度に影響していると考えられる。
コクランの規則に従ってセルを合併すると近似の精度は上がる。

2.4 シミュレーション.

1000回繰り返して、自由度９のカイ二乗分布に近づくかをみる。

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(0, 10, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1]),h$ counts[2:10])
  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott")
curve(dchisq(x,9),xlim=c(0,20),add=T)

#fit <- fitdistr(a, densfun = "gamma")
#fit

理論分布と少しずれがあるように感じられる。
そこで、コクランの規則に従って修正を行う」。

2.4.1 コクランの規則.

ピアソンの定理は、\(Np_i\geq 5\) （𝑖 = 1,2, ⋯ , 𝑘） であれば、カイ二乗分布による近似がよいと言われている。

しかし、histによる分割では \(Np_i< 5\) となる区間も出てくる。 そこで、区間を結合するなどして、期待値が 5 未満のセルの割合を20%以下にするというのがコクランの規則である。

［コクランの規則］Cochran’s rule.

• 期待値が 5 未満のセルが全体の20％以上にならないようにする.

具体的には、
期待度数５以下の区分が２０％を超えないよう左端の区間と右端の区間を合併する.

2.4.2 左端を合併.

a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(0, 10, by = 1),plot=FALSE)

y <- 100*c(1-sum(qnom[3:10]),qnom[3:10])
str(y)

##  num [1:9] 5.47 11.72 20.51 24.61 20.51 ...

sum(y<5)/9

## [1] 0.3333333

まだ 5 未満のセルが多いようである。
シミュレーションすると、

K <- c()

for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(0, 10, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1:2]),h$ counts[3:10])
  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott")
curve(dchisq(x,8),xlim=c(0,20),add=T)

まだずれがあるように感じられる。両端を合併してみよう。

2.4.3 両端を合併.

a <- rbinom(100,10,0.5)#二項乱数
h <- hist(a,breaks=seq(0, 10, by = 1),plot=FALSE)
#累積確率
pnom <- pbinom(h$breaks, 10, 0.5)#累積確率
qnom <-diff(pnom, lag = 1)#差分
y <- 100*c(1-sum(qnom[3:10]),qnom[3:7],sum(qnom[8:10]))
str(y)

##  num [1:7] 5.47 11.72 20.51 24.61 20.51 ...

sum(y<5)/7

## [1] 0

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(0, 10, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1:2]),h$ counts[3:7],sum(h$ counts[8:10]))
  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T)
curve(dchisq(x,6),xlim=c(0,20),add=T)

自由度６のカイ二乗分布の理論曲線はヒストグラムによく当てはまっており、近似の精度は良くなったように思う。

参考サイト　生物科学研究所 井口研究室 様.

2.4.4 推定母数があるとき.

期待度数を求めるために観測値からr個の母数を推定していれば、自由度は r 減る。

　　　\(T=\displaystyle \sum_{i=1}^{k}\dfrac{(X_i-Np_i)^2}{Np_i}\rightarrow \chi^2(k-r-1)\).

p=0.5 の代わりに、推定値 p = mean(a)/10 を用いる

• 理論累積確率 pbinom(h$breaks, 10, 0.5)　のところを.

• 推定値を用いて　pbinom(h$breaks, 10, mean(a)/10).

としてシミュレーションすると自由度は１減るはずである.

K <- c()
for(i in 1:1000){
  a <- rbinom(100,10,0.5)#二項乱数
  h <- hist(a,breaks=seq(0, 10, by = 1), plot = FALSE)
  x <- c(sum(h$ counts[1:2]),h$ counts[3:7],sum(h$ counts[8:10]))
  
  pnom <- pbinom(h$breaks, 10, mean(a)/10)#累積確率
  qnom <-diff(pnom, lag = 1)#差分
  y <- 100*c(1-sum(qnom[3:10]),qnom[3:7],sum(qnom[8:10]))

  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T)
curve(dchisq(x,5),xlim=c(0,20),add=T)

出力乱数から、平均を mean(a)/10 と推定しているので、自由度は１ 減っているはずである。
自由度５のカイ二乗分布の理論曲線はヒストグラムによく当てはまっている。

3 正規分布.

3.1 乱数によるデータ.

a <- rnorm(100,50,10)#正規乱数
h <- hist(a,breaks=seq(0, 100, by = 10),xaxp=c(0, 100,10))

str(h)

## List of 6
##  $ breaks  : num [1:11] 0 10 20 30 40 50 60 70 80 90 ...
##  $ counts  : int [1:10] 0 1 3 9 31 36 16 4 0 0
##  $ density : num [1:10] 0 0.001 0.003 0.009 0.031 0.036 0.016 0.004 0 0
##  $ mids    : num [1:10] 5 15 25 35 45 55 65 75 85 95
##  $ xname   : chr "a"
##  $ equidist: logi TRUE
##  - attr(*, "class")= chr "histogram"

3.2 理論確率の分布.

#理論確率
pnom <- pnorm(h$breaks, 50, 10)#累積確率
qnom <-diff(pnom, lag = 1)#差分
y <- 100*c(1-sum(qnom[4:10]),qnom[4:7],sum(qnom[8:10]))
str(y)

##  num [1:6] 2.28 13.59 34.13 34.13 13.59 ...

sum(y<5)/6

## [1] 0.3333333

度数５未満の区分の割合は、コクランの規則を少し超えている。

3.3 シミュレーション.

1000回繰り返して、自由度５のカイ二乗分布に近づくかをみる。

K <- c()
for(i in 1:1000){
  a <- rnorm(100,50,10)#正規乱数
  h <- hist(a,breaks=seq(0, 100, by = 10), plot = FALSE)
  x <- c(sum(h$ counts[1:3]),h$ counts[4:7],sum(h$ counts[8:10]))
  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,xlim=c(0,20),prob=T,breaks="scott")
curve(dchisq(x,5),xlim=c(0,20),add=T)

mean(K)

## [1] 5.048954

自由度５のカイ二乗分布は、かなりよく当てはまっている。

3.3.1 推定母数があるとき.

• 理論累積確率　pnorm(h$breaks, 50, 10)　のところを.
• 推定値を用いて　pnorm(h$breaks, mean(a), sd(a)).

としてシミュレーションすると、自由度は２減るはずである.

K <- c()
for(i in 1:1000){
  a <- rnorm(100,50,10)
  h <- hist(a,breaks=seq(0, 100, by = 10), plot = FALSE)
  x <- c(sum(h$ counts[1:3]),h$ counts[4:7],sum(h$ counts[8:10]))
  
  pnom <- pnorm(h$breaks, mean(a), sd(a))#平均と標準偏差に推定値を用いる
  qnom <-diff(pnom, lag = 1)#差分
  y <- 100*c(1-sum(qnom[4:10]),qnom[4:7],sum(qnom[8:10]))

  k <- sum((x-y)^2/y)#カイ二乗値
  K <- c(K,k)
}
hist(K,prob=T)
curve(dchisq(x,3),add=T)

mean(K)

## [1] 3.219646

期待度数のところで、観測値から平均と標準偏差を推定して用いている。
従って、自由度は２減り 6-1-2=3 である。
カイ二乗値の平均も３に近くなっている。

4 分割表.

分割表の適合度について、当サイト 適合度検定 から引用する。

4.1 カイ二乗値の計算例

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

d <- rmultinom(n=1, size = 50, prob = c(1,5,3,1))  # 乱数を発生
t(d)#転置

##      [,1] [,2] [,3] [,4]
## [1,]    7   20   19    4

m <-50*c(1,5,3,1)/10;m #期待度数

## [1]  5 25 15  5

#(d-m)
#(d-m)^2
#(d-m)^2/m
sum((d-m)^2/m) #乖離度

## [1] 3.066667

4.2 カイ二乗値の分布.

ピアソンの定理によれば、
４水準の分割表から計算されるカイ二乗値は、自由度３のカイ二乗分布に従う。
size=50の実験を、3000回出力し、理論分布と比較する。

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)

size=20 としてみる

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)

やはりヒストグラムに大きな変化は見られない。
分割表の理論分布からの乖離度であるカイ二乗値は、

• 多項分布の確率分布に左右されない.
• 水準数にのみ依存する. 自由度 （水準数−１）のカイ二乗分布に従う。

水準の数はカテゴリーのまとめ方によって変わる。
サイコロの例で、カテゴリーのまとめ方を変えてシミュレーションする。

4.3 サイコロ60回

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)

驚くべきことに、カテゴリーのまとめ方に左右されない。水準数のみに依存する。

まとめ.
分割表のカイ二乗値の分布は、

• 確率分布に左右されず.
• 水準の数のみに左右される.

4.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.8413333

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

5 四分表.

四分表のカイ二乗値についてシミュレーションする。
詳しくは、当サイト 番外編 分割表の解析 を参照のこと。.

	肺がんあり	肺がんなし	計
A:喫煙あり	a	b	a+b
B:喫煙なし	c	d	c+d
合計	a+c	b+d	a+b+c+d

四分表の検定は、
　　　\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\)

が自由度１のカイ二乗分布 \(\chi^2(1)\)に従うことに基づいて行われる。.

この検定量は、観測値から周辺分布の確率を、
　　　\(\pi_1=\dfrac{a+b}{N}\) 　　　\(\pi_2=\dfrac{a+c}{N}\).　　によって推定し、

独立という仮説のもとで 理論度数と実現値との乖離度を計算したものである。
　

＜独立という仮説のもとでの確率分布＞.

確率分布	肺がんあり	肺がんなし	計
A：喫煙あり	\(\pi_1\pi_2\)	\(\pi_1(1-\pi_2)\)	\(\pi_1\)
B：喫煙なし	\((1-\pi_1)\pi_2\)	\((1-\pi_1)(1-\pi_2)\)	\(1-\pi_1\)
合計	\(\pi_2\)	\(1-\pi_2\)	1

この統計量 T は、「四分表の要因が独立である」 という状態からの乖離度を表している。 このことについて、乱数シミュレーションで実験してみよう。

5.1 例５−４.

市原清志先生の「バイオサイエンスの統計学」p125の例題５−４を参照.

白血病の患者140名をランダムに２群に分け、A,B２種の方法で治療するとき、治療後の結果を比較する。 これまでの説明に合わせ、次のような表にする。

	生存	死亡	計
治療 A	48	32	80
治療 B	24	36	60
合計	72	68	140

目的は「治療法A,Bに差があるか」を検定することである。

検定は、

• 帰無仮説「差がない」のもとで、調査を無限回行い統計量の分布を描いたとき、
  
• 上記の表から得られる統計量が、統計的に有意な位置にあるかどうかをみる。

無限回の乱数をどのように出力すれば良いだろうか？

\(N=140\) は固定しておこう。

5.1.1 二項乱数を２種類使用（１）.

周辺分布（合計の分布）から、生存確率を \(p_1=\dfrac{72}{140}\)　とし、.
治療A 80人、治療B 60人　を固定しておくことにすると、

a1 <- rbinom(n=1, size = 80, prob = c(72/140))
a2 <- rbinom(n=1, size = 60, prob = c(72/140))
x <- c(a1,80-a1,a2,60-a2)
x

## [1] 36 44 33 27

期待度数と乖離度.

#期待度数
y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
round(y,0)

## [1] 41 39 31 29

abs(x-y)^2/y

## [1] 0.6428571 0.6806723 0.1488095 0.1575630

#乖離度
sum(abs(x-y)^2/y)

## [1] 1.629902

これらは手計算で確かめることができる。

問題に与えられた調査データ (48,32,24,68) の、期待度数からの乖離度は.

x <- c(48,32,24,36)
s1 <- sum(abs(x-y)^2/y)
s1

## [1] 5.490196

シミュレーション.

1000回出力してみよう。

X = c()
for (i in 1:1000) {
  a1 <- rbinom(n=1, size = 80, prob = c(72/140))
  a2 <- rbinom(n=1, size = 60, prob = c(72/140))
  x <- c(a1,80-a1,a2,60-a2)
  #print(x)
  y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
}
mean(X)

## [1] 1.905077

hist(X,prob=T)
curve(dchisq(x,2),add=T)
abline(v = s1)

#1-pchisq(s1, 2)

自由度が２になっている理由は、

• ２種類の二項乱数からデータを出力している.
• 各二項乱数は、2つの区分を生成するので自由度は１.
• 検定量は各二項乱数の乖離度の和であるので、
  
• カイ二乗分布の再生性から、自由度は 1+1=2.

5.1.2 二項乱数を２種類使用（２）.

治療A 80人、治療B 60人　を固定しておくが、
生存確率を \(A:p_1=\dfrac{24}{60}~~~B:p_2=\dfrac{48}{80}\) のように異なる値に設定しても、得られる統計量の分布は \(\chi^2(2)\)である。.

a1 <- rbinom(n=1, size = 80, prob = c(24/60))
a2 <- rbinom(n=1, size = 60, prob = c(48/80))
x <- c(a1,80-a1,a2,60-a2)
x

## [1] 37 43 36 24

期待度数と乖離度.

#期待度数
y <- c(80*c(24/60,36/60),60*c(48/80,32/80))
round(y,0)

## [1] 32 48 36 24

abs(x-y)^2/y

## [1] 0.7812500 0.5208333 0.0000000 0.0000000

#乖離度
sum(abs(x-y)^2/y)

## [1] 1.302083

これらは手計算で確かめることができる。

調査データ (48,32,24,68) の乖離度は.

x <- c(48,32,24,36)
s1 <- sum(abs(x-y)^2/y)
s1

## [1] 23.33333

理論確率を変えたので、与えられたデータの乖離度は大きくなっている。

シミュレーション.

1000回出力してみよう。

X = c()
for (i in 1:1000) {
  a1 <- rbinom(n=1, size = 80, prob = c(24/60))
  a2 <- rbinom(n=1, size = 60, prob = c(48/80))
  x <- c(a1,80-a1,a2,60-a2)
  #print(x)
  y <- c(80*c(24/60,36/60),60*c(48/80,32/80))
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
}
mean(X)

## [1] 2.114809

hist(X,prob=T)
curve(dchisq(x,2),add=T)
abline(v = s1)

1-pchisq(s1, 2)

## [1] 8.574939e-06

２組の２分表のカイ二乗値の和なので、カイ二乗分布の再生性により自由度２のカイ二乗分布になる。

5.1.3 生存確率に推定値を用いる.

生存確率に推定値を用いれば自由度は１になるはずである。

x <- c(48,32,24,36)
a1=48;a2=24
y <- c(80*c((a1+a2)/140,1-(a1+a2)/140),60*c((a1+a2)/140,1-(a1+a2)/140))
s1 <- sum(abs(x-y)^2/y)
s1

## [1] 5.490196

X = c()
for (i in 1:1000) {
  a1 <- rbinom(n=1, size = 80, prob = c(72/140))
  a2 <- rbinom(n=1, size = 60, prob = c(72/140))
  x <- c(a1,80-a1,a2,60-a2)
  #print(x)
  y <- c(80*c((a1+a2)/140,1-(a1+a2)/140),60*c((a1+a2)/140,1-(a1+a2)/140))
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
}
mean(X)

## [1] 0.9358842

hist(X,prob=T)
curve(dchisq(x,1),add=T)
abline(v = s1)

#1-pchisq(s1, 2)

確かに、自由度は１になっている。

5.2 検定量 T を用いる.

四分表の検定量　\(T=\dfrac{(ad-bc)^2N}{(a+b)(c+d)(a+c)(b+d)}\) 　を用いると、自由度は１になる。.

注目すべきことは、

• 二項乱数を２種類用いて生成したデータでも.
• 多項乱数を用いて生成したデータでも.

自由度が１になっていることを確かめることができる。

x <- c(48,32,24,36)
a=x[1];b=x[2];c=x[3];d=x[4]
s1 <- (a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
s1

## [1] 5.490196

5.2.1 二項乱数２種類で生成.

X = c()
for (i in 1:1000) {
  a1 <- rbinom(n=1, size = 80, prob = c(72/140))
  a2 <- rbinom(n=1, size = 60, prob = c(72/140))
  x <- c(a1,80-a1,a2,60-a2)
  a=x[1];b=x[2];c=x[3];d=x[4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)

## [1] 1.062764

hist(X,prob=T)
curve(dchisq(x,1),add=T)
abline(v = s1)

1-pchisq(s1, 1)

## [1] 0.0191234

5.2.2 多項乱数で生成.

当サイト 計数値データの検定　参照.

p1=72/140;p2=80/140
#c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2))
d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2))) 
# 乱数を発生
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")
data[1,]

a	b	c	d
31	48	31	30

X = c()
for (i in 1:1000) {
  a=data[i,1]
  b=data[i,2]
  c=data[i,3]
  d=data[i,4]
  kai=(a*d-b*c)^2*140/((a+b)*(c+d)*(a+c)*(b+d))
  X <-c(X,kai)
}
mean(X)

## [1] 1.048524

同じデータでも、次のように母比率を用いると自由度は３になる。
これは、
検定量 T を用いた検定は、

• 母比率からの乖離度ではなく.
• 独立性からの乖離度を表している.

つまり、独立かどうかだけを検定していることになる。

p1=72/140;p2=80/140
d <- rmultinom(n=1000, size = 140, prob = c(p1*p2,p2*(1-p1),p1*(1-p2),(1-p1)*(1-p2))) 
data <- as.data.frame(t(d))
colnames(data) <- c("a", "b", "c", "d")

X = c()
for (i in 1:1000) {
  a=data[i,1]
  b=data[i,2]
  c=data[i,3]
  d=data[i,4]
  x <- c(a,b,c,d)
  #print(x)
  y <- c(80*c(72/140,68/140),60*c(72/140,68/140))
  s <- sum((abs(x-y))^2/y)
  X <-c(X,s)
}
mean(X)

## [1] 2.990447

hist(X,prob=T)
curve(dchisq(x,3),add=T)
abline(v = s1)

1-pchisq(s1, 3)

## [1] 0.1392262

自由度が３になっている.
しかも、乖離度は有意ではない。

6 感想.

• 推定した母数を用いて計算すると自由度が１減じるのは、不偏分散の計算で母平均を標本平均で推定すると自由度が１減じるのと似ている. このため不偏分散の計算では、nで割るところを n-1 で割るのであった。

• 帰無仮説のもとでは、二項乱数でも多項乱数でも簡易計算値は自由度１のカイ二乗分布に従っていることは驚きであり、びっくり ピアソン である。

参考サイト.

東北大学大学院 尾畑研究室.

高校情報Ⅱ 作成途中. http://hyoccori-toukei.la.coocan.jp/id-2/id-3/id/id1.html

看護応用 ichihara 四分表.

fitdist関数.

library(fitdistrplus)
a <- rnorm(100,50,10)#正規乱数
#h <- hist(a,breaks=seq(0, 100, by = 10),xaxp=c(0, 100,10))
#str(h)
fit <- fitdistr(a, densfun = "normal")
fit

##       mean          sd    
##   50.9384036   10.2296751 
##  ( 1.0229675) ( 0.7233473)