ブートストラップについての調べもの

Tags:
Categories:

Rでブートストラップについて調べたときの記録。

ブートストラップとは

大雑把に言えば, 統計量とかの分布を調べたりするために1つの標本から新しい標本（ブートストラップ標本）を何個も作ることのようだ。パラメトリックとノンパラメトリックがある。

今回扱うのはノンパラメトリックの方。

自作関数によるブートストラップ

ブートストラップ標本を作るときには, もとの標本から同じサイズで復元抽出を行う。復元抽出とは要するに, くじ引きでくじを引くたびにもとの箱にくじを戻すことである。たとえば, 【太郎さん、花子さん、和夫さん】の3人からなるもとの標本があった場合に, あるブートストラップ標本では【太郎さん, 太郎さん, 花子さん】になるかもしれないし, 別のブートストラップ標本では【和夫さん, 和夫さん, 和夫さん】となるかもしれない。Rで復元抽出を行う際にはsamaple()関数の引数でreplace=TRUEを指定すると良い。試しにやってみる。Rの組み込みのデータセットであるstateからアメリカ50州の収入のデータを使って, 平均値と中央値がどのように分布するかを見てみる。

# もとの標本
original <- state.x77[,"Income"]

# 元の標本から復元抽出する回数
num.bs.sample <- 1000

# 結果を保存するようベクトル
median.vector <- numeric(num.bs.sample)
mean.vector <- numeric(num.bs.sample)

for (i in 1:num.bs.sample){
  median.vector[i] <- median(sample(original, length(original),
                                    replace=TRUE))
  mean.vector[i] <- mean(sample(original, length(original),
                                replace=T))
}

# 結果の図示
boxplot(median.vector, mean.vector, names=c("Median","Mean"),
        horizontal = T, las=T, xlab="Income")

結果は次のとおりである。これを使うと中央値の経験的な分布なども求められる。

boot1.png

bootパッケージによるブートストラップ

毎回, 自分で繰り返しのループを書くのも面倒であろうということでbootという便利なパッケージがある。boot(もとの標本, 統計量を求める関数, ブートストラップ標本数)のように使う使うらしい。同じデータでやってみよう。

library(boot)

foo <- function(d,k){
  median.Income <- median(d[k])
  mean.Income <- mean(d[k])
  return(c(median.Income,mean.Income))
}

res.boot <- boot(original, foo, R=1000)

res.bootはリスト形式で見てみるとこんな結果が出てくる。

> res.boot

ORDINARY NONPARAMETRIC BOOTSTRAP

Call:
boot(data = original, statistic = foo, R = 1000)

Bootstrap Statistics :
    original   bias    std. error
t1*   4519.0 -20.9665   105.96233
t2*   4435.8  -1.6806    85.47084

オリジナルというのがもとの標本で求めた統計量（ここでは中央値と平均値）である。biasというのは, ブートストラップ標本の統計量の期待値からもとの標本の統計量を引いたもので, 真の値から外れている程度を表すらしい。次にあるブートストラップにより推定された標準誤差と大きさを比べて判断していくようである。

この結果を箱ひけにすると以下の通り。

boot2.png

【参考にしたサイト】

ブートストラップ