ブートストラップ法 (Rでの書き方メモ)

ブートストラップ法に関して
今日学んだことについて、まとめておく。
初学者ですので間違えていたらすみません。

○ bootstrap
ブートストラップは、ある統計量が、どういう分布をとるか調べるのに使います。
なお、その統計量自体は、標本集団から不偏統計量を計算すればよいです。
分布がわかれば、標準誤差や、信頼区間がわかり、検定にも使えます。

方法は、標本集団を、1.標本サイズ同じ、2.要素の重複可 (Replace=True )でリサンプリング(＝ブートストラップして)
それぞれにおいて統計量を求め、それをその統計量の取りうる分布とします。

※ ブートストラップの注意点

ブートストラップでうまく分布が求められない統計量がある。
たとえば、最大値を求める max() 関数
理由1：現在の標本集団の最大値より、大きくならない。
理由2：標本集団に存在する値のところに、スパイクができてしまう。沢山ブートストラップしても消えないスパイク。

library(boot)

# 以下のように統計量を返す関数を用意します。
# data と indicies を受け取れるようにします。
# 第3引数以降も設定できます。 その場合は boot関数に、パラメーター引数として渡します。
# indiciesには、sample( n, replace=T ) のような感じのものが毎回渡されます。
statistic_func = function( data, indicies ){   
	quantile( data[indicies, ]$cost , probs = c(0.75) )
}


data(nuclear) # nuclearというdataを読み込みます


# 1000回 ブートストラップして、統計量の分布を生成します。
result = boot( nuclear, statistic_func , R=1000 )


# resultの中身は、boot クラスというもので、
# 今回、行った1000回のブートストラップ結果は、result$t に入っています。

# 分布のプロットは
plot(density(result$t))

# 分布のS.D.は
sd(result$t[,1])

# 分布の信頼区間は
boot.ci(result)