図1 10000個のブートストラップ標本平均値のヒストグラム
<< 目次へ
最終更新日： 2009.7.15
シミュレーション入門
■ モンテカルロシミュレーション
シミュレーションとは擬似実験のことである。シミュレーションの方法の1つにモンテカルロ法というものがあるが、これは乱数を用いて同じ実験（試行）を何度も繰り返し行うというものである。例えば、「コインを投げる」という試行について考えてみよう。
まずコイン投げをする簡単なプログラムを用意する。
coin <- function(){
x <- runif(1) # 0〜1までの乱数を発生させる
if(x <= 0.5) result <- 1 # 0.5以下なら表(1)とする
else result <- 0 # それ以外なら裏(0)とする
result # 結果を返す
}
このプログラムを実行することは、自分がコインをトスすることと同じだと考えればよい。
> coin() # 1回目の試行 [1] 1 > coin() # 2回目の試行 [1] 1 > coin() # 3回目の試行 [1] 1 > coin() # 4回目の試行 [1] 1 > coin() # 5回目の試行 [1] 0 > coin() # 6回目の試行 [1] 0
しかし、これではたくさんの試行を行うのには時間がかかりすぎる。そこでfor構文を使えば、難なく数百回、あるいは数千回という試行を瞬時に行うことができる。
my.repeat <- function(k){
a <- numeric(k) # k回の試行分のデータを格納する変数
for(i in 1:k) a[i] <- coin() # coin()をk回繰り返す
tbl <- table(a) # 表と裏の出た回数を集計する
tbl # 集計結果を返す
}
引数には繰り返しの数を指定する。例えば100回、1000回、10000回と行ってみる。
> my.repeat(100) # 100回の試行 a 0 1 47 53 > my.repeat(1000) # 1000回の試行 a 0 1 504 496 > my.repeat(10000) # 10000回の試行 a 0 1 4988 5012
このように同じ実験を繰り返し行うことによって、真の値（表が出る頻度と裏が出る頻度）へと近似させることができる。これがモンテカルロシミュレーションである。
■ ブートストラップ法によるシミュレーション
その平均値は正しいか？
N(mean=50, sd=10)に従う母集団から大きさ10（n=10）の標本値が得られているとしよう。
> n1 <- rnorm(10, mean=50, sd=10) > n1 [1] 35.51247 55.83346 62.82574 50.96489 53.52258 44.87819 52.66387 53.44522 [9] 61.14276 45.77619
ここで、標本から得られる平均値（標本平均）は真値（母平均）と等しいかどうかについて知りたいとしよう。古典的には「標本平均は母平均と同じである」という帰無仮説についてt検定が行われる。
# 「標本平均は母平均50と等しい」という帰無仮説について検定を行う
> t.test(n1, mu=50)
One Sample t-test
data: n1
t = 0.6517, df = 9, p-value = 0.5309
alternative hypothesis: true mean is not equal to 50
95 percent confidence interval:
45.90634 57.40673
sample estimates:
mean of x
51.65654
統計的仮説検定によれば、p = 0.5309なので「true mean is not equal to 50」という対立仮説は採択されることになる。そして95%の信頼区間は[45.91, 57.41]であることから、95%の確率で母平均がこの区間に含まれるということがわかる。
ブートストラップ標本
先ほどの例で得られた信頼区間を別のアプローチで求めることができる。それがブートストラップ法による信頼区間の推定なのであるが、それを説明する前にブートストラップについての基本を解説しておこう。
母集団から無作為に抽出された標本があるならば、そこから復元抽出によって大きさnの標本をk個抽出することができる（復元抽出なので理論的にnはいくらでも多く取れるし、何個の標本を抽出するかも自由である）。このような手順によって得られた2次的な標本をブートストラップ標本という。
先ほどのN(mean=50, sd=10)という母集団から抽出した大きさ10の標本からブートストラップ標本を作成し、その平均値を求めてみよう。なお、ブートストラップ標本の抽出はsample()という関数を利用すれば簡単に行える。
# 復元抽出を許すには、引数にreplace=TRUEと指定する > n2 <- sample(n1, 10, replace=TRUE) > n2 [1] 53.44522 45.77619 50.96489 53.52258 35.51247 53.44522 44.87819 52.66387 [9] 35.51247 52.66387 > mean(n2) [1] 47.8385
このような作業を、例えば10000回繰り返し行ったとしよう。すると10000個のブートストラップ標本が得られ、それぞれについて平均値を求めるのだから、結局は10000個のブートストラップ標本平均値が得られることになる（この10000回繰り返し行うというのは、すでに紹介したモンテカルロシミュレーションに他ならない）。そして10000個の平均値をヒストグラムで表してみよう（図1）。
boot.res <- numeric(10000)
for(i in 1:10000){
boot.n <- sample(n1, 10, replace=TRUE)
boot.res[i] <- mean(boot.n)
}
hist(boot.res)
図1 10000個のブートストラップ標本平均値のヒストグラム
さて図中に描かれている赤いラインは95%信頼区間の下側境界と上側境界である。これはquantile()という関数を用いることによって簡単に得ることができる。quantile()は四分位数を求めるための関数である。
> quantile(boot.res, c(0.025, 0.975))
2.5% 97.5%
46.77542 56.19036
回帰分析への応用
回帰分析に対してブートストラップ法を用いる方法はいくつかあるが、ここでは最も分かりやすい適用例をあげてみる。それは観測値（x, y）というデータを復元抽出したブートストラップ標本に対して回帰分析を行うというものである。
例としてRに組み込まれているデータセットcarsを使用してみる。これはdata()という関数を用いることによって呼び出すことができる。また、解析に用いる際にはデータセットを呼び出した後にはattach()を使ってパスを通しておく必要がある。
> data(cars) # carsを呼び出す > names(cars) # carsに含まれる変数を確認 [1] "speed" "dist" > attach(cars) # パスを通す
まずは通常通り回帰分析を行い、定数項と回帰係数の95%信頼区間を求めてみる。
> result <- lm(dist ~ speed)
> summary(result)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791 6.7584 -2.601 0.0123 *
speed 3.9324 0.4155 9.464 1.49e-12 ***
> -17.5791 - 1.96 * 6.7584 # 定数項の下側境界値
[1] -30.82556
> -17.5791 + 1.96 * 6.7584 # 上側境界値
[1] -4.332636
> 3.9324 - 1.96 * 0.4155 # 回帰係数の下側境界値
[1] 3.11802
> 3.9324 + 1.96 * 0.4155 # 上側境界値
[1] 4.74678
次にn=50のブートストラップ標本を10000個用意し、それぞれに対して回帰分析を行う。そうすると10000個の定数項と10000個の回帰係数が得られるので、それぞれヒストグラムを作成してみる。ただし、これをやるには少し時間がかかる（環境によるが2〜3分ほどかかる）。
intr <- numeric(10000)
coef <- numeric(10000)
for(i in 1:10000){
# speedとdistのペアデータを復元抽出する
a <- cars[sample(1:50, 50, replace=TRUE), ]
# 回帰分析の結果を代入する
res <- lm(a[,2] ~ a[,1])
intr[i] <- res$coefficients[1] # 定数項の取り出し
coef[i] <- res$coefficients[2] # 係数の取り出し
}
--- ---
# グラフウィンドウを1行2列に分割する
> par(mfrow=c(1, 2))
> hist(intr) # 定数項のヒストグラム
> quantile(intr, c(0.025, 0.975)) # 定数項の95%信頼区間
2.5% 97.5%
-29.371159 -6.883976
# 図中に95%の信頼区間の境界値を書き加える
> abline(v=c(-29.371159, -6.883976), col="red")
> hist(coef) # 係数のヒストグラム
> quantile(coef, c(0.025, 0.975)) # 係数の信頼区間
2.5% 97.5%
3.151969 4.760125
# 信頼区間の境界値を書く
> abline(v=c(3.151969, 4.760125), col="red")
図2 定数項と回帰係数のヒストグラム
| 定数項 | 回帰係数 | |
| 正規分布と仮定した信頼区間 | [-30.82556, -4.332636] | [3.11802, 4.74678] |
| ブートストラップによる信頼区間 | [-29.371159, -6.883976] | [3.151969, 4.760125] |
ブートストラップ法のコツ
ブートストラップ標本を得るとき、どのくらいの大きさの標本をいくつくらい用意すればよいだろうか、と考えるかもしれない。ブートストラップの特徴は、ある推定値（平均値や回帰係数など）をたくさん得ることによって中心極限定理を生じさせることにある。したがって、標本の大きさよりも、むしろ繰り返しの数を増やした方がより良いと考えられるだろう。つまりブートストラップ標本をたくさん用意するわけである。ブートストラップ標本の大きさについては扱う問題にもよるが、今回取り扱ったような例ではとりあえずn=20くらいにしておけば問題ないだろう。特に根拠はないが、元の標本と同じ大きさにしておけばよい様に思える。
ブートストラップ法の利点
