ブートストラップ法
諸般の事情でブートストラップ法を利用する可能性が高いので復習をかねて書きます。
こちらがすごくまとまっていたので、参考にしました。
Web上であまり情報が見つからなかったのは探し方が悪かったのかな??
とりあえず
- パラメトリック・ブートストラップ法と、ノンパラメトリック・ブートストラップ法がある。
- 有名、というかよく使われるのはノンパラメトリック・ブートストラップ法の方で、今回書いているのもたぶんノンパラメトリック・ブートストラップの方
ブートストラップ法とは
標本集団からリサンプリングを繰り返し(重複を許す)、得られた新たな標本集団(ブートストラップ標本)の統計量の分布が、母集団の分布に近いものになる、という性質を利用して、母集団に対する事前知識なし(確率密度関数を使わず)に、母集団の統計量を推定する手法です。
確認
Rのリハビリをかねつつ、確認してみます。
適当に与えたデータ(1〜10)までを標本とし、それを重複を許しながらリサンプリングして、結果をヒスト
グラムで表示するプログラムです。
# 標本データ x <- c(1,2,3,4,5,6,7,8,9,10) # とりあえず統計量を求めてみる x.stat <- c(mean(x), var(x)) # boot.num回のリサンプリングによる推定 boot.num <- 1000 x.boot <- numeric(boot.num) # 10個の標本データから重複を許して10個をサンプリング # 平均を求める for (i in 1:boot.num) { x.boot[i] <- mean(sample(x, 10, replace=T)) } # ヒストグラムによるグラフの描画 hist(x.boot) # 標本平均と真の平均値の誤差 # 実はこんな計算をしなくても、平均は一定なので # 見かけの分布はx.bootのヒストグラムでも確認できる x.clt <- x.boot - mean(x) # ヒストグラムによるグラフの描画 hist(x.clt) # パーセンタイル点の計算 # 95%信頼区間を計算 print(quantile(x.boot, c(0.025, 0.975))) # 通常のhist関数では平滑化し過ぎていてあまりよろしくない # らしいので、オプションをつけて使う。 hist(x.boot, breaks="Scott", xlim=c(2.5,8.5), ylim=c(0,0.5), prob=T, ann=F) par(new=T) plot(density(x.boot), xlim=c(2.5,8.5), ylim=c(0,0.5), xlab="", ylab="", main="", col="red")
当たり前ですが、いい感じでグラフが表示されます。
ついでにRの勉強もかねて、hist関数をちょちょいと変形させたり、確率密度の推定曲線を重ね合わせて書いています。
コードの後半部分です。
参考書に書いてあることが喜楽に試せるのがRのいいところですよね。
理解が深まります。
参考
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/49.html
中心極限定理 - Wikipedia
信頼区間 - Wikipedia
このPDFが一番しっかりしていると思うけど、まだ理解しきれていない
Rによるやさしい統計学 | |
山田 剛史 杉澤 武俊 村井 潤一郎 オーム社 2008-01-25 売り上げランキング : 3880 Amazonで詳しく見る by G-Tools |