ブートストラップの適切なサンプル数 -「データ解析のための統計モデリング入門」第5章 #みどりぼん

勉強会 みどりぼん

7/8 に開催された「データ解析のための統計モデリング入門」、通称「みどりぼん」の第 5 回読書会にのこのこ参加。主宰のやまかつさん、発表者＆参加者の皆さん、会場を提供してくださったドワンゴさん、ありがとうございました。


今回は第５章、尤度比検定。
みどりぼんは、この前の回の AIC も含めたモデル選択について比較的ゆるふわな説明で、これはこれでありなんだけど、しっかりしゃっきりしたものも読んでみたいなあ（でも数式少なめでよろしく！）という人は、伊庭さんの「モデル選択とその周辺」がおすすめ度高し。


検定は（有用だろうけど）あんまり好きになれなくて、気分的な盛り上がりは今ひとつなのだが、知らなかったら好きも嫌いもないよね！　ということでまじまじ手を動かす。
勉強会の質疑応答も、前回参加した第３回よりおとなしめだったかな。モデル選択って自明な例見せられても実感わかないよねー。


そんな中、パラメトリック・ブートストラップについて、p104 の脚注には「実はサンプルの個数は 10^3 ぐらいでは十分なサイズではありません」とあるが、適切なサンプルサイズはいくつ？　そりゃあ増やすほどいいんだろうけどどこまで増やす？　とかいう質問が上がっていた。
サンプルを増やすほど推定量の精度が上がる（＝分散が下がる）ことは、直感的にわかるだろうし、事実その通りなのだが、残念なことに分散の下がり方は遅い。
実際、一番シンプルな問題として「サンプルの平均」を平均の推定量とする問題を考えるとき、簡単な計算ですぐわかるのだが、サンプル数を N 倍するとその推定量の分散は 1/N に、つまり標準偏差は 1/√N となる（分布の形が漸近的に正規分布に一致することを言うには中心極限定理が必要だけど、ここではそこまでの話は求めていない）。
つまり頑張ってサンプルを 100倍にしても精度は 10倍にしか上がらず、どこかでコストが見合わなくなるわけで、そこが「どこまで増やすか」の答えとなる。

と、いうのが机上のお話なわけだが、「1000 じゃあ足りないよね〜」と訳知り顔でうなずいている人たちは常にそういう事を意識してるかというと必ずしもそういうわけではなく、もっと単純に 1000 では足りない事を経験として知っているのだ。


みどりぼんの 5.4.1 では、パラメトリック・ブートストラップで 1000 個の逸脱度の差を得て、その分布の 95% 点を 3.953957 と推定するお話が出てくる。
ランダムサンプリングしている以上、実行するたびに得られる値は異なる（つまりこの推定値は確率変数）わけだが、さてはて、実際どれくらい結果に幅があるのだろうか。
コードが全て掲載されているのだから、実行して確かめてみればいい。

d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv") 
get.dd <- function(d) {
  n.sample <- nrow(d);
  y.mean <- mean(d$y);
  d$y.rnd <- rpois(n.sample, lambda=y.mean)
  fit1 <- glm(y.rnd~1, data=d, family=poisson);
  fit2 <- glm(y.rnd~x, data=d, family=poisson);
  fit1$deviance - fit2$deviance;
}
pb <- function(d, n.bootstrap) {
  replicate(n.bootstrap, get.dd(d));
}


ここまで準備をしておけば、あとは quantile(pb(d, 1000), 0.95) と叩くだけで毎回異なる推定値が得られる。
ぜひ、答えを見る前に、みどりぼんの中の値 3.953957 と比べて、どれくらい推定値に幅が出るか予想してみよう。
テキトーでいいよ〜。



OK？　では実際に5回実行して得られた5つの推定値がこちら。

> quantile(pb(d, 1000), 0.95)
     95% 
4.165724 

> quantile(pb(d, 1000), 0.95)
     95% 
4.312189 

> quantile(pb(d, 1000), 0.95)
     95% 
3.513216 

> quantile(pb(d, 1000), 0.95)
     95% 
4.064914 

> quantile(pb(d, 1000), 0.95)
     95% 
3.751253


おそらく多くの人が「そんなに変わるんかー……」と目を丸くしてくれたのではと期待しているのだが、どうだろう。
これでもう次からは「1000 じゃあ足りないよね〜」といっしょにうなずく仲間入り。


ただし、みどりぼん 5.4.1 の目的は「逸脱度の差 4.5 の p 値が 0.05 より小さいか」を確認することだったのを思い出すと、上の推定値は最大でも 4.3 ちょっとと、p<0.05 という結論には影響ない。
つまり、今回の問題にとって必要な精度は 1000 サンプルでも確保できてそう（まあ、もう少しくらい精度高いほうが安心できるだろうけど）。
しかし逆に、もうちょっときわどい値だったら、例えば「逸脱度の差 4.0 の p 値が 0.05 より小さいか」という問題だったときに、1000サンプルで推定した値がたまたま 3.953957 で、「あーよかった、ギリギリ超えたから棄却できた！」なぁんて言っちゃうたら……。
そういったことを適切に判断することが一番大事で大変だったりする。