ポアソン分布に従うか？
100日間、1日の事故の回数を計測した観測結果はポアソン分布に従うと考えられるか？
# グラフ
# 事故回数100日分
y <- c(rep(0,43),rep(1,31),rep(2,14),
rep(3,8),rep(4,3),rep(5,1))
y
(n <- length(y))
# histgramの横軸
(x <- 0:max(y))
# ポアソン分布のパラメータ（標本平均）
(theta <- sum(y)/n)
hist(y, xlim=range(x),ylim=c(0,1),probability=T,nclass=5,col="yellow")
# ポアソン分布の確率密度
lines(x, dpois(x,theta),col="red", lwd=3)
ポアソン分布
確率変数のとりうる値が離散かつ無限個
単位時間あたり、平均$\lambda$回発生する事象が
単位時間に$x$回（$0$〜)発生する確率は
$$P(x)=\frac{e^{-\lambda}\lambda^x}{x!}$$
$p$が微小で$n$がきわめて大きいとき、$\lambda=np$
ポアソン分布との適合度
$H_0:$ポアソン分布に従う
$H_1:$ポアソン分布に従うとはいえない
$H_0$を仮定した場合、理論的には データ数×ポアソン分布の確率密度 になる。
ただ、実際には偶然を考えて完全なポアソン分布にはならない。
では、偶然性を考慮した上で、ポアソン分布に従っているか？
観測したデータを観測度数のテーブル
$H_0$に基づいた期待度数のテーブル
この２つを作る。
期待度数の最小値 $<$ 5 の場合、近似がうまくいかないので注意
今回は期待度数の最小値が5より小さくなるセルは集約する
観測度数と期待度数のズレの和を考える
このズレの和はカイ二乗値と等しくなる。
適合度検定
$H_0$が正しければ、nが大きい場合、
$$\chi^2=\sum\frac{(O_i-E_i)^2}{E_i}$$
- パラメータλを与える（推定しない）場合
自由度：$k-1$ のカイ二乗分布に従う。 - パラメータλを推定るす場合
自由度：$k-1-1$ のカイ二乗分布に従う。
$O_i$:観測度数テーブル$O$のi番目
$E_i$:期待度数テーブル$E$のi番目
$k$：属性数
$H_0$のもとでこの値は小さいはず。
小さいならば$H_0$を支持
大きいなら$H_0$を棄却し、$H_1$を支持。
R
# ポアソン分布に従うか
(event_count = c(0,1,2,3,4,5))
(observe_count = c(43,31,14,8,3,1))
d <-data.frame(
observe = c(rep(0,43),rep(1,31),rep(2,14),rep(3,8),rep(4,3),rep(5,1))
)
# 格納確認
dim(d)
head(d)
names(d)
summary(d)
# table関数からクロス集計表を作る
(obs_tabs <- table(d))
# xtabsを使ってもtableを作れる
# (obs_tabs <- xtabs(~observe,d))
# 合計付きの集計表へ
f_tabs <- addmargins(obs_tabs)
f_tabs
# 期待度数のtableを作る
(Expected <- obs_tabs) # 観測度数の表とdimensionが同じになるので、いったん格納して
(Expected <- obs_tabs*0) # 期待度数を入れるために初期化
# ポアソン分布に従うとして期待度数のtableを作る
# sumを拾う
(n <- f_tabs["Sum"])
# ポアソン分布のパラメータλを推定する場合(標本平均)
(lambda_hat <- sum(event_count*observe_count)/sum(observe_count))
# n*ポアソン分布の確率密度でExpectedに格納
for (i in 1:dim(obs_tabs))
Expected[i] <- n*dpois(i-1, lambda = lambda_hat)
# ポアソン分布のパラメータを推定しない場合
#lambda_given <- 0.8
# n*ポアソン分布の確率密度でExpectedに格納
# for (i in 1:dim(obs_tabs))
# Expected[i] <- n*dpois(i-1, lambda = lambda_given)
# 以降もlambda_hatをlambda_givenに変更
# 確率密度で計算しているので、合計値がこの段階では異なる
sum(obs_tabs)
sum(Expected)
# 対策1
# これを解消するために、一番端は一個手前までの累積をで引く。
# 一個前までの累積
sum(Expected[1:dim(Expected)-1])
# 累積を引く
Expected[dim(Expected)] <- sum(obs_tabs) - sum(Expected[1:dim(obs_tabs)-1])
Expected
# ここがTRUEで返れば正しい期待度数が作れている
sum(Expected) == sum(obs_tabs)
# 対策2
# 一番端は、端から2番目までの累積確率を引いた確率に訂正
Expected[dim(Expected)] <- n*(1-ppois(dim(Expected)-2, lambda= lambda_hat))
Expected
# ここがTRUEで返れば正しい期待度数が作れている
sum(Expected) == sum(obs_tabs)
# (期待度数の最小値) ≧ 5 を満たしてないので3以上で集約する
# データフレームで4以上を3に書き換える
# 期待度数の最小値を満たしていればこのまま
# chisq = sum((obs_tabs-Expected)^2/Expected)
# chisq
# これで終了
d_fix <- d
# 4以上を3に書き換えることで3以上を一番端にする
d_fix$observe[d_fix$observe >= 4] <- 3
d_fix
########## ここから下はデータフレームを書き換えて
########## 同じことをしているだけです。
# d_fixにして再度期待度数のテーブルを作る
# d_fixをもう一度dに返して再利用
d <- d_fix
(obs_tabs <- table(d))
# xtabsを使ってもtableを作れる
# (obs_tabs <- xtabs(~observe,d))
# 合計付きの集計表へ
f_tabs <- addmargins(obs_tabs)
f_tabs
# 期待度数のtableを作る
(Expected <- obs_tabs) # 観測度数の表とdiimensionが同じになるので、いったん格納して
(Expected <- obs_tabs*0) # 期待度数を入れるために初期化
# n*ポアソン分布の確率密度でExpectedに格納
for (i in 1:dim(obs_tabs)) Expected[i] <- n*dpois(i-1, lambda = lambda_hat)
Expected
# 確率密度で計算しているので、合計値がこの段階では異なる
sum(obs_tabs)
sum(Expected)
# 対策1
# これを解消するために、一番端は一個手前までの累積を使う。
# 一個前までの累積
sum(Expected[1:dim(Expected)-1])
Expected[dim(Expected)] <- sum(obs_tabs) - sum(Expected[1:dim(obs_tabs)-1])
Expected
# ここがTRUEで返れば正しい期待度数が作れている
sum(Expected) == sum(obs_tabs)
# 対策2
# 一番端は、端から2番目までの累積確率を引いた確率に訂正
Expected[dim(Expected)] <- n*(1-ppois(dim(Expected)-2, lambda= lambda_hat))
Expected
# ここがTRUEで返れば正しい期待度数が作れている
sum(Expected) == sum(obs_tabs)
# 集約によって期待度数の最小値が5以上になりクリア
# 観察度数と期待度数のズレの和
chisq = sum((obs_tabs-Expected)^2/Expected)
chisq
$\chi^2$ = 4.971824
$H_0$：ポアソン分布に従っている、について
両側$5$%水準で検討する。
自由度：$4−1-1=2$のカイ二乗分布の上側$5$%点が$5.99$
棄却限界値$5.99$よりもカイ二乗値が小さいので
$H_0$は棄却されず、$H_0$が支持される。
パラメータの推定あり、なし
$\lambda$推定なしの場合、適合度検定の結果は与えられたlambdaによるポアソン分布に従うかどうかを検定することになる。
$\lambda$推定ありの場合は、（一応）ポアソン分布に最も当てはまるだろうと思われる状況でポアソン分布に従うかどうかを検定することになる。
カテゴリーまとめの基準について
期待度数が5以上になるように設定しましたが、その基準がそもそも怪しい可能性があり、勉強中です。