読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ポアソン分布のカイ二乗適合度検定. R による練習.

カイ二乗適合度検定の枠組み

カイ二乗検定は分割表に対して行うもの, というイメージがある. たぶんそのイメージはまちがっていない.

タイプ 観測値の個数 生起確率
A_1 X_1 p_1
A_2 X_2 p_2
: : :
A_k X_k p_k
  • ある実験を行った結果が  A_1, \ldots , A_kk 通りに分類されるとする.
  •  P(A_i) = p_i,  i=1, \ldots, k かつ p_i + \cdots + p_k = 1 とする.
  •  A_i に分類される観測値の個数を  X_i とし, X_1 + X_2 + \cdots + X_k =n と表す.

これは  (X_1, X_2, X_k) が, パラメータ  n, p=(p_1,\ldots,p_k) の多項分布(多項分布 - Wikipedia)に従っているということだ.

そしてカイ二乗適合度検定は,

  • 帰無仮説: (p_{o1}, p_{o2}, \ldots, p_{ok})= (p_{1}, p_{2}, \ldots, p_{k})
  • 対立仮説: (p_{o1}, p_{o2}, \ldots, p_{ok)} \neq (p_{1}, p_{2}, \ldots, p_{k})

について次の検定統計量,

 \displaystyle Q(X) = \sum^{k}_{i=1}\frac{(X_i- np_{i})^2}{np_{i}}

(観測度数引く理論度数の二乗割る理論度数の総和)

を評価する.

有意確率( p 値)カイ二乗分布の上側確率で求められる.

ケース・スタディ

ラザフォードとガイガーのポロニウムデータ

Rutherford and Geiger (1910) による原子核崩壊で表れたポロニウムポロニウム - Wikipedia)の数を観測したデータは, ポアソン分布によく適合することが知られている.

統計ソフト R では VGAMパッケージにこのデータが入っている.

R: Rutherford-Geiger Polonium Data

data(ruge,package="VGAM")
> ruge
   counts number
1      57      0
2     203      1
3     383      2
4     525      3
5     532      4
6     408      5
7     273      6
8     139      7
9      45      8
10     27      9
11     10     10
12      4     11
13      0     12
14      1     13
15      1     14

f:id:abrahamcow:20150114160729p:plain

(N = with(ruge, sum(counts)))
lambdahat = with(ruge, weighted.mean(number, w=counts))
fit1 <-  with(ruge, cbind(number, counts,
                          N * dpois(number, lambda=lambdahat)))
plot(fit1[,1:2],type="h",lwd=2)
points(fit1[,c(1,3)], pch=4, col="red3",lwd=2,cex=1.5)

#nihongo()
legend("topright", c("観測度数", "理論度数"), col = c("black","red3"),
         lty = c(1, -1), pch = c(-1, 4),lwd=2, merge = TRUE)

ポアソン・モデルの適合度検定

ruga データを下表のようにまとめ, カイ二乗適合度検定を行ってみる.

number k counts x_k
0 57
1 203
2 383
3 525
4 532
5 408
6 273
7 139
8 45
9 27
10+ 16

 k \ge 10 をまとめて,  n = 2608 の11項分布だと考える.
カイ二乗分布でよく近似されるためには, 各セルの期待度数が 3 以上であることが望ましいらしい. その理由は知らない. )

tmp <- with(ruge,counts[number>=10])
cnt <- c(ruge$counts[ruge$number<10],sum(tmp))
#length(cnt)
ruge2 <- data.frame(number=c(0:9,"10+"),counts =cnt)

各項の生起確率がポアソン分布に従うという帰無仮説は,

\displaystyle p_k = \frac{e^{-\lambda} \lambda^k}{k !},

 k=0,1, \ldots, k;  p_{10} =1-p_0-p_1- \cdots, -p_0

という仮定である.

ポアソン分布のパラメータ λ の有効推定量(推定量と分散について 2 (クラメル・ラオの不等式) | singular point)は, 重み付き平均,

\displaystyle \hat{\lambda} = \sum^{\infty}_{k=0} \frac{k x_k}{n}

である. 今回 k \ge 11x_k はぜんぶ 0.

( lambdahat <- sum(0:10*cnt)/ N) ) #10086 / 2608
#[1] 3.871549
#lambdahat <- with(ruge2, weighted.mean(0:10, w=cnt))

カイ二乗検定の検定統計量を求めると,

k <- 0:9
p_k <- ( exp(-lambdahat)*lambdahat^k )/ factorial(k)

p_k[11] <- 1- sum(p_k)

(chi2 <- sum(((cnt - fitted)^2)/fitted))
#[1] 12.91345

自由度は,

(項の数-1)- 当てはめたモデルのパラメータの数
=11 - 1 - 1
= 9

ポアソン分布のパラメータは λ 一個)

だから,

pchisq(chi2,9,lower.tail=FALSE)
#[1] 0.1665617

p 値は約0.17.

ポアソンモデルの適合は, 通常の有意水準では棄却されない.

ポアソン分布を当てはめることとこの結果は矛盾しない.

別解

ポアソン分布の場合, 平均と分散が等しい(ともに λ )ため, 標本(不偏)分散  s^2 と標本平均  \bar{x} を用いて,

 \displaystyle \frac{(n-1) s^2}{\bar{x}}

を検定統計量とすることができる.

7.2節 適合度検定を参照.

X <- rep(with(ruge, rep(number, counts)))
chi2 <- (var(X)*(N-1))/mean(X)
pchisq(chi2,df=N-1,lower.tail=FALSE)
#[1] 0.9506711

結び

交代再生過程のある時刻における状態がポアソン分布に従うか否かシミュレーションしてみた. - 廿TT で, 適合度の検定まちがえまくったので, 改めて勉強中です.

他人には「A/Bテストで仮説検定をブラックボックス的に使うな!」とえらそうに再三主張しておきながら, 自分は R をわけもわからず使ってるんだから世話はない.

お恥ずかしい限りです.

参考文献

入門・演習 数理統計

入門・演習 数理統計

(主にpp.283-284 を参照した.)

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

(主にpp.149-151 を参照した.)