ブートストラップ(bootstrap)は統計学者Efronが(考案した|流行らせた)計算方法です。
コンピュータによるシミュレーションの一種ですが,データを乱数で生成するのではなく,実際のデータに基づいた推論ができるのが,ブートストラップの特徴です。
例えば平均値と中央値(メディアン)のどちらが安定か(ぶれが少ないか)を調べるために,通常のシミュレーションでは「乱数でデータを生成して平均値と中央値を計算する」という作業を何度も繰り返し,その分布を比較します。これに対して,ブートストラップでは,「実際のデータ(n 個)からランダムに n 個を復元抽出(重複を許した抽出)し,その平均値と中央値を計算する」という作業を何度も繰り返します。例えば実際のデータが (2, 3, 5, 7) であれば,それから復元抽出した (2, 3, 3, 7), (3, 5, 7, 7), (2, 5, 5, 5) などについて統計量(平均値,中央値など)を計算し,その分布を調べます。
この驚くほど単純な方法で,従来は正規分布など数学的に理想化された分布についてしか求めることができなかった統計量の分布(例えば中央値の分布)を,現実のデータとほぼ同じ分布について求めることができるのです。
ブートストラップについてはたくさんの本が出ていますが,考案者Efronと弟子のTibshiraniが書いた入門書 An Introduction to the Bootstrap (1993) がわかりやすいと思います。三重大学の図書館にもあります。
ブートストラップ用のライブラリを使わない簡単な例です。
データとしては,あまり意味がありませんが,昔の鳩山内閣の資産データが私のブログ 閣僚資産:平均かメジアンか にありますので,これを使って,閣僚本人の資産の中央値の標準偏差を求めてみましょう。
X = read.csv("shisan.csv")
N = 1000
m = numeric(N)
for (i in 1:N) {
k = sample(18, replace=TRUE)
m[i] = median(X$本人[k])
}
これで長さ1000の配列 m にブートストラップされた中央値が入ります。この標準偏差を sd(m) で求めてみれば,だいたい1000程度の値が得られます。
上のプログラムは,わかりやすいようにループで書きましたが,
m = replicate(1000, { k = sample(18, replace=TRUE); median(X$本人[k]) } )
または
m = replicate(1000, median(X$本人[sample(18,replace=TRUE)]))
のような書き方もできます。
Rのブートストラップ用ライブラリとしては,上で紹介した入門書
An Introduction to the Bootstrap
に基づいた bootstrap
や,DavisonとHinkleyの Bootstrap Methods and Their Applications (1997)
に基づいた boot があります。後者がより一般的です。
> library(boot)
> b = boot(X, function(d,k) median(d$本人[k]), R=100000)
> b
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = X, statistic = function(d, k) median(d$本人[k]),
R = 1e+05)
Bootstrap Statistics :
original bias std. error
t1* 2681 -32.10432 1005.341
この2行目 b = ... は,次のように分けて書くほうがわかりやすいかもしれません:
foo = function(d, k) {
median(d$本人[k])
}
b = boot(X, foo, R=100000)
boot() の第1引数は全データを含むベクトルまたは行列またはデータフレームで,第2引数で与える関数の第1引数に引き継がれます。第2引数で与える関数は統計量を求めるためのもので,その第2引数はデータの行インデックスです(ライブラリを使わない例で sample(18,replace=TRUE)
を入れたものです)。引数 R=... がブートストラップ・サンプルの数です。boot() の第2引数で与えた関数はR回呼び出されます。Rは通常1000回あれば十分です。
最初の original = 2681 が元のサンプルについての統計量(ここでは中央値)で,これはブートストラップとは関係せず単に median(X$本人) を求めたものです。
次の bias = -32.10432 は,中央値の期待値が中央値の真の値からどれだけ外れているかをブートストラップで推測したものです。こんなに有効桁数があるわけではありませんが,中央値の期待値が真の中央値より数十程度小さい可能性があることを示しています。具体的には,ブートストラップによる統計量の平均値から,元のサンプルの統計量を引いたものです。ここでは,このバイアスの値は,次に述べる標準誤差より十分小さいので,無視することができますが,そうでない場合は統計量に問題があるのかもしれません。
最後の std. error が中央値の標準誤差(standard error)をブートストラップで推定したもので,中央値はサンプルごとにこれくらいの標準偏差のばらつきがあることを示しています。具体的には,ブートストラップによる統計量の標準偏差(n - 1 で割る方式)を求めているだけです。
結果を変数(ここでは b)に入れておくと,後でより詳しい結果を出力することができます。例えば plot(b) とすればブートストラップされた統計量の分布がプロットされます。
par(mgp=c(2,0.8,0)) plot(b)
Last modified: 2011-11-12 22:00:53