先日の投稿で紹介した、ポアソン過程の観測値がxの時の期待値λの分布が形状母数x+1, 尺度母数1のガンマ分布に従う性質の応用編。
計数データの誤差を推定するため、データがポアソン過程で得られていると仮定して、観測値の信頼区間を求めることにした。
Rの場合、簡易的にはpoisson.test()で求めることができるらしい(奥村先生のページ)。
たとえば観測値が1の時の期待値の95%信頼区間は以下のようにして0.025 - 5.572とされる。
> poisson.test(x = 1, conf.level = 0.95)
Exact Poisson test
data: 1 time base: 1
number of events = 1, time base = 1, p-value = 1
alternative hypothesis: true event rate is not equal to 1
95 percent confidence interval:
0.02531781 5.57164339
sample estimates:
event rate
1
しかし、poisson.test()は複数の観測値に対して信頼区間を計算するようにできていない。
また、信頼区間以外にもいろいろ計算していて信頼区間だけを取り出す手間が必要になる。
そこで、ポアソン過程の観測値がxの時の期待値λの分布が形状母数x+1, 尺度母数1のガンマ分布に従うことを利用し、qgamma()を用いる。
qgamma()はガンマ分布の左側p%のquantileを求めてくれる関数。
poisson.test()内部でも用いられている。
x <- 0:10 #観測値
x.L <- qgamma(0.025, x + 1) #95%信頼区間の下側
x.H <- qgamma(0.975, x + 1) #95%信頼区間の上側
knitr::kable(data.frame(x, x.L, x.H)) #表に出力
| x | x.L | x.H |
|---|---|---|
| 0 | 0.0253178 | 3.688880 |
| 1 | 0.2422093 | 5.571643 |
| 2 | 0.6186721 | 7.224688 |
| 3 | 1.0898654 | 8.767273 |
| 4 | 1.6234864 | 10.241589 |
| 5 | 2.2018943 | 11.668332 |
| 6 | 2.8143631 | 13.059474 |
| 7 | 3.4538322 | 14.422675 |
| 8 | 4.1153731 | 15.763189 |
| 9 | 4.7953887 | 17.084803 |
| 10 | 5.4911604 | 18.390356 |
観測値が0 - 10の場合の95%信頼区間を推定できた!
と、思いきや、先のpoisson.test(1)の結果
| x | x.L | x.H |
|---|---|---|
| 1 | 0.02531781 | 5.57164339 |
と比較すると何かおかしい。
poisson.test()で推定されたx = 1の時のx.Lは、qgamma()で推定されたx = 0の時のx.Lと一致している。
x.Hは問題ない。
どなたか理由ご存知の方いらっしゃれば教えて下さい。