「統計学関連なんでもあり」の過去ログ--- 043

No.11047　誤差がポアソン分布に従うときの予測区間の算出法　　【Sai】　2009/10/09(Fri) 23:12

いつもお世話になっております。
過去ログ，ネット，書物を調べましたが，見当たらなかったので質問させてください。
今日は，予測区間(prediction confidence)の算出法についてご教授していただきたく思います。

例えばRで以下のような例を用意します。

set.seed(123)
x <- seq(1, 100)
y <- 2 * x + rpois(100, 20)
plot(x, y)

まず，誤差が正規分布に従うとした場合の予測区間は，

rr <- glm(y ~ x) # 回帰
v <- (x - mean(x)) ^ 2
ssx <- sum(v)
s <- sqrt((1 / (100 - 2)) * sum((rr$residuals - mean(rr$residuals)) ^ 2))
s * sqrt((1 / length(y)) + (v / ssx)) # 標準誤差
predict(rr, se.fit=TRUE)$se.fit # 確認

lines(x, predict(rr), col = 2) # 予測値
lines(x, predict(rr) + 1.96 * predict(rr, se.fit = TRUE)$se.fit, col=3) # 予測区間上限
lines(x, predict(rr) - 1.96 * predict(rr, se.fit = TRUE)$se.fit, col=3) # 予測区間下限

で表されると思います。

では，誤差がポアソン分布でリンク関数にlogを用いた場合，どのように計算されるのでしょうか。私としては，ポアソン分布は平均と分散が等しいため，以下のように考えました。

rr2 <- glm(y ~ x, family = poisson(link = "log")) # 回帰
x11()
plot(x, y)
lines(x, exp(predict(rr2)), col = 4) # 予測値
lines(x, exp(predict(rr2) + (predict(rr2) * predict(rr2, se.fit=TRUE)$se.fit)), col = 3) # 予測区間上限
lines(x, exp(predict(rr2) - (predict(rr2) * predict(rr2, se.fit=TRUE)$se.fit)), col = 3) # 予測区間下限

しかし，どうも合っているのか確証が持てません。もし，どなたか，ポアソン分布の場合の予測区間の算出法を知っていましたら，ご教授願えませんでしょうか。関連事項が記載された書物，Rパッケージの名前だけでも構いません。
どうぞよろしくお願い致します。

No.11048　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【後医は名医】　2009/10/10(Sat) 01:36

石村貞夫著「統計解析のはなし」（東京図書）p179にポアソン分布の母平均λの区間推定の計算式が記載されています。

No.11049　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【知ったかぶり】　2009/10/13(Tue) 10:12

予測区間と信頼区間は別物です．どちらを算出したいのでしょうか（信頼区間のようですが）．

質問内容とは直接関係ありませんが，例示されている乱数はちょっと不適当．「誤差がポアソン分布に従う」ということについて誤解があるのでは．
yがポアソン分布に従い，xとyの間にy=2xという関係があるのであれば，
y<-rpois(100,2*x)
とし，link関数にはidentityを指定するべきでしょう．link="log"とするのであれば，例えば，
y<-rpois(100,exp(0.05*x))

No.11062　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【Sai】　2009/10/14(Wed) 13:07

＞後医は名医さん
ありがとうございます。早速調べてみます。

＞知ったかぶりさん
やや，本当ですね。これは今まで大きな勘違いをしていたかもしれません。
これは，正規分布のときも一緒ですね。ただ，正規分布はどこでも分散が等しいので，上記のように
x <- seq(1, 100)
y <- 2*x+rnorm(100)
plot(x, y)
とやっても，問題なさそうですね。ポアソン分布だと場所によって分散が大きく違うようです。
y2 <- rpois(100, 2*x)
plot(x, y2)

そして確かに，これは予測区間ではなく信頼区間という名前で使われていますね。
ただ，自分の中では予測区間の定義が，データから推定した平均値が95%の信頼度で含まれる範囲と認識していましたので，定義が逆になっていただけのようです。

そして，改めていいますと，私が知りたいのは信頼区間のことです。
自分でも後医は名医さんの本を探していますが，近くにないようですので，取り寄せを考えています。もし，簡単な記述なのであれば，ご紹介していただけると幸いに思います。

No.11063　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【青木繁伸】　2009/10/14(Wed) 13:24

ポアソン分布のパラメータλの信頼区間
> poisson.conf <- function(x, alpha=0.05)
+ {
+ 	N <- length(x)
+ 	df <- 2*sum(x)+2
+ 	cat(sprintf("%f ≦ λ ≦ %f\n", qchisq(alpha/2, df)/2/N, qchisq(1-alpha/2, df)/2/N))
+ }
> poisson.conf(c(2,1,2,3,4,3))
1.524230 ≦ λ ≦ 4.123370
> set.seed(1)
> poisson.conf(rpois(50, lambda=5))
4.643444 ≦ λ ≦ 5.914431

No.11066　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【知ったかぶり】　2009/10/14(Wed) 17:58

線形予測子は最尤推定量なので，漸近正規であるとして，

rr2.pred<-predict(rr2, se.fit=TRUE)
lines(x, exp(rr2.pred$fit + 1.96 * rr2.pred$se.fit), col = 3) # 信頼区間上限
lines(x, exp(rr2.pred$fit - 1.96 * rr2.pred$se.fit), col = 3) # 信頼区間下限

で良いのではないかと思いますが.どなたかフォローを.

No.11094　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【Sai】　2009/10/18(Sun) 22:00

＞青木先生
むむ，カイ二乗分布ですか…。やはりこれは本を見ないとなぜこうなるのかわからないですね。しかし，これは知ったかぶりさんの出した信頼区間とは違う価になる…のでしょうか？？ちょっと自分で計算できませんでした。すいません，理解不足のようです。

＞知ったかぶりさん
なるほど，パラメータは最尤推定しており，それは漸近正規を仮定しているので，ポアソン分布ではなく，正規分布でよいということですか…。確かに以下のpptでも同じことを言っています。
http://bm.hus.osaka-u.ac.jp/~torii/poisson.ppt

なるほど…しかし，モデル当てはめの際には，データの誤差にポアソン分布を仮定しているのに，予測のときにはデータの誤差に正規分布を仮定しているみたいで若干気持ち悪いですね。このことを証明している本，URL（きちんとした著者名のあるもの）などあればご紹介していただけないでしょうか。

No.11095　Re: 誤差がポアソン分布に従うときの予測区間の算出法　　【青木繁伸】　2009/10/19(Mon) 08:35

> これは本を見ないとなぜこうなるのかわからないですね

その本には，式しか書かれていません。

また，「ポアソン分布のパラメータλの信頼区間」を求めるものですよ。λは，平均値です。
ようするに，あなたが必要としているものとは違います。