Rでブートストラップ法

「ブートストラップ」という言葉は聞いたことはあるけど、具体的にどんなことなのかはよく知らない、というのが私のレベルです。

が、「統計学入門」（東京大学教養学部統計学教室 編）という教科書↓

統計学入門 (基礎統計学)


の中の演習問題に出てきた（でも、教科書の本編には解説はない）ので、Rでちょっとやってみました。

第3章 練習問題
3.4　＜ブートストラップ＞　以下の手法をパソコンで行え。

i)　　[1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

ii) 　i)を11回繰返して実行し、11個のランダムな番号を得たのち、表3.1、図3.1のデータからその番号のデータを取り出し（ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える）、相関係数rを計算せよ。

iii)　ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップという。

「表3.1、図3.1」というのは、教科書の本編のところに出てきた表と図で、相関関係が認められる例として紹介されていたものです。

統計学入門 (基礎統計学) より

統計学入門 (基礎統計学) より

このデータを使う必要があるので、まずは入力しましょう。

> bro <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66)
> sis <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62)

書籍に載っていた図と同じような感じでプロットしてみました↓

> plot(bro, sis, xlim=c(68,70), ylim=c(59,70), asp=1, pch=16)

書籍に載っていたデータをRでプロットしたもの

書籍に出てくる図とほぼ同じものができました。ただ、点が10個しかなく1個少ないです。これは(70, 65)という全く同じデータが2つあって同じ座標にプロットされているためです。書籍の方では少しずらして、分かるようになっていました。

あと、書籍だと(64,62)あたりに点がありますが（一番左側のもの）、これは間違いでしょうね。点を数えてみると12個あるし・・・

使うデータは入力したので、演習問題をやっていきましょう。

i)　　[1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

乱数の従う分布について特に言及していないので、一様分布でしょうね。

まず思いつくのは、

> ceiling(runif(1, 0, 11))   # 1～11の整数の乱数

runif関数が返す値をceiling関数で整数化してやる感じですね。

でも、Rで整数の乱数を発生させるんだったら、sample関数を使った方がシンプルな気がします↓

> sample(1:11, 1, replace=T) # 1～11の整数の乱数

1, 2, ..., 11 の中から1個選ぶ、復元抽出（replace=T）で、という感じですね。

ii) 　i)を11回繰返して実行し、11個のランダムな番号を得たのち、

これはsample関数の引数を変えるだけでOKですね。

> rnd <- sample(1:11, 11, replace=T)
> rnd
 [1]  8  1  8  2  4 11  5  3  9  8 10

表3.1、図3.1のデータからその番号のデータを取り出し（ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える）、

番号に応じたデータを取り出すには、先ほど発生させたrndをインデックスに指定ればいいですね。

> bro[rnd]
 [1] 73 71 73 68 67 66 70 66 72 73 65

> sis[rnd]
 [1] 64 69 64 64 63 62 65 65 66 64 59

相関係数rを計算せよ。

> cor(bro[rnd], sis[rnd])
[1] 0.5357521

iii)　ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップ（太字）という。

> r <- numeric(200) # 200回分の格納用
>
> for ( i in 1:200 ) {
+   rnd <- sample(1:11, 11, replace=T)
+   r[i] <- cor(bro[rnd], sis[rnd])
+ }
>
> hist(r)

ブートストラップ法で作成したヒストグラム

これで宿題は終わりました。単位は取れそうです。

元のデータの相関係数をそのまま出すと、

> cor(bro, sis)
[1] 0.5580547

になるのですが、ヒストグラムのモードは違う場所（0.65あたり）に出てますね。分布も左に歪んでいるし。

以上の統計手法をブートストラップという。

とのことですが、ヒストグラムの意味（モードのずれや歪み）や使い方については、もうちょっと勉強しないといけなさそうです。