「適合度の検定について」に対し,中澤さんからコメントをいただいた。
> vcd パッケージの goodfit が便利です
はじめて使ってみました。
> goodfit(as.table(tab1))
Observed and fitted values for poisson distribution
with parameters estimated by `ML'
count observed fitted
0 0 2.150092e-03
1 0 2.805870e-02
2 0 1.830830e-01
3 1 7.964111e-01
4 1 2.598291e+00
5 3 6.781540e+00
6 12 1.474985e+01
7 27 2.749793e+01
8 42 4.485600e+01
9 63 6.504121e+01
10 80 8.487877e+01
11 108 1.006971e+02
12 129 1.095081e+02
13 108 1.099293e+02
14 106 1.024698e+02
15 87 8.914871e+01
16 79 7.271192e+01
17 45 5.581709e+01
18 42 4.046739e+01
19 31 2.779471e+01
20 18 1.813605e+01
21 8 1.127026e+01
22 5 6.685312e+00
23 1 3.793188e+00
24 1 2.062546e+00
25 2 1.076649e+00
26 1 5.403950e-01
> summary(goodfit(as.table(tab1)))
Goodness-of-fit test for poisson distribution
X^2 df P(> X^2)
Likelihood Ratio 19.73707 22 0.5994826
fitted の合計が 999.5218 となるのがちょっといやな感じです(x >= 27 に対する確率が考慮されていない)。
goodfit がやっていることは,以下の通り(冗長な部分を削除)。
x = 3:26
tab1 = c(1, 1, 3, 12, 27, 42, 63, 80, 108, 129, 108, 106, 87, 79, 45, 42, 31, 18, 8, 5, 1, 1, 2, 1)
names(tab1) = x
lambda = sum(x * tab1) / sum(tab1)
p = dpois(x, lambda)
cbind(x, tab1, p, p*sum(tab1))
ln = function(o, e) sum(ifelse(o == 0, 0, o*log(o/e)))
( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
( df = length(x) - 1 -1 )
( p.value = pchisq(G2, df, lower.tail=FALSE) )
結果
> ( G2 = sum(ln(tab1, p*sum(tab1)))*2 )
[1] 19.73707
> ( df = length(x) - 1 -1 )
[1] 22
> ( p.value = pchisq(G2, df, lower.tail=FALSE) )
[1] 0.5994826
G2(G スクエアー)を使う検定(独立性の検定でも,G2 を使う検定が定義されている)。
この検定では,観察値が 0 のセルは検定統計量に反映されないので,観察値が 0 のセルがあっても良いし,「期待値が小さいセルについての考慮」も不要ということにはなっている。
※コメント投稿者のブログIDはブログ作成者のみに通知されます