1. はじめに
ブートストラップ信頼区間について調べていたんですが、理論的な求め方は教科書などに載っているのですが、実践的な情報が少ないように思います。
今回、少し調査してみて、実際に適用する際に注意が必要だなと感じたことについて書いておきます。
2. ブートストラップ信頼区間
ブートストラップ法は、理論的に求めるのが難しい統計量を、経験分布からのシンプルなリサンプリングによって推定できるという手法です。
ブートストラップ法では、推定された統計量に対して、その信頼区間を求めることもできます。
このような信頼区間をブートストラップ信頼区間といいます。
ブートストラップ信頼区間を求める方法については色々議論があるようですが、主な手法は次の5つです。
- 正規分布近似法
- ベーシック法
- パーセンタイル法
- BCa法(bias-corrected and accelerated percentile method)
- ブートストラップ-t 法
ここでは、ブートストラップ-t 法を除く4つの手法について、分散の推定値に対する信頼区間を求め、R によるシミュレーションによりその性質を調べてみることにします。
3. R でブートストラップ信頼区間を求める
まず、適当なサンプルを用意します。標準正規分布のサンプルにします。
x <- rnorm(100)
経験分布による分散の推定値を求める関数を定義します。
calcVar <- function(x) { n <- length(x) var(x) * (n - 1)/n }
R にはブートストラップを行うパッケージがいくつかあります。ここでは、simpleboot パッケージを使って信頼区間を求めてみます。
library(simpleboot) boot.out <- one.boot(x, calcVar, 10000) boot.ci(boot.out, conf = 0.95, type = c("norm", "basic", "perc", "bca"))
Intervals : Level Normal Basic 95% ( 0.6309, 1.0187 ) ( 0.6273, 1.0104 ) Level Percentile BCa 95% ( 0.6188, 1.0018 ) ( 0.6508, 1.0582 )
このように、正規分布近似法(Normal)、ベーシック法(Basic)、パーセンタイル法(Percentile)、BCa 法について、分散のブートストラップ信頼区間を求めることができました。
ところで分散の信頼区間の理論値は下記で求められるため*1、ブートストラップ法で求める必然性はありません。
calcTheoreticalCI <- function(x, conf) { n <- length(x) alpha <- (1 - conf)/2 var(x) * (n - 1)/qchisq(c(1 - alpha, alpha), df = n - 1) } calcTheoreticalCI(x, 0.95)
[1] 0.6610565 1.0572708
ここでは、理論値と比較するために、あえて分散の信頼区間をブートストラップ法で計算しています。
4. ブートストラップ信頼区間の性質
ブートストラップ法により求められる信頼区間の性質を調べてみます。
まず、信頼区間は分布の両端を求めなければならず、端に行くほど精度が落ちるはずです。求める信頼区間を 10%, 20%, …, 90% と変更してみて、信頼区間にどれだけ真値 (=1) が入るか確認してみます。*2
conf <- seq(0.1, 0.9, by = 0.1) params <- data.frame(conf = conf, N = 10, B = 1000) result <- data.frame(conf, calcResult(params)) library(ggplot2) ggplot(result, aes(x = conf, y = prob, col = group)) + geom_line(size = 1) + guides(col = guide_legend(title = NULL)) + ylim(0, 1)
expect は設定した値、theoretical は理論値です。
計算結果が expect に近いほど、推定精度が良いということになります。
理論値はさすがに推定精度が良く、他は同じぐらいかなと思います。
思ったとおり、端に行くほどブートストラップ信頼区間の精度は落ちるようです。*3
「性質1:求める信頼区間の区間幅が大きいほど推定精度が悪くなる」
次に、サンプルサイズ N を変化させてみます。N = 10, 20, …, 100 まで変えてみます。
N <- seq(10, 100, by = 10) params <- data.frame(conf = 0.95, N = N, B = 1000) result <- data.frame(N, calcResult(params)) library(ggplot2) ggplot(result, aes(x = N, y = prob, col = group)) + geom_line(size = 1) + guides(col = guide_legend(title = NULL)) + ylim(0, 1)
N が大きいほど精度は良くなるようです。
これは、サンプルサイズが多いほど、経験分布は確率分布に近づくためだと考えられます。
「性質2:サンプルサイズが多いほど推定精度は良くなる」
次に、ブートストラップ反復回数 B を変化させてみます。B = 100, 200, …, 1000 としてみます。
B <- seq(100, 1000, by = 100) params <- data.frame(conf = 0.95, N = 100, B = B) result <- data.frame(B, calcResult(params)) library(ggplot2) ggplot(result, aes(x = B, y = prob, col = group)) + geom_line(size = 1) + guides(col = guide_legend(title = NULL)) + ylim(0, 1)
どうやら精度は B によらないようです。教科書には、ブートストラップ信頼区間を求めるには、反復回数 B = 1000 〜 2000 は必要と書かれていますが、おかしいですね。*4
もう少し小さい B で試してみます。
B がごく小さい場合では、B によってブートストラップ信頼区間の精度が向上しているのがわかります。*5
「性質3:ブートストラップ反復回数はごく小さい場合を除き、推定精度には関係しない」
ところで、N が小さいときに、B を増やせば精度が上がる、ということは無いようです。
B <- seq(200, 2000, by = 200) params <- data.frame(conf = 0.95, N = 10, B = B) result <- data.frame(B, calcResult(params)) library(ggplot2) ggplot(result, aes(x = B, y = prob, col = group)) + geom_line(size = 1) + guides(col = guide_legend(title = NULL)) + ylim(0, 1)
経験分布と確率分布の違いが大きい場合はブートストラップでどんなに頑張っても信頼区間を求めるのは難しそうです。
「性質4:サンプルサイズが小さいときに、ブートストラップ反復数を増やせば推定精度が上がる、ということは無い」
5. まとめ
今回わかったブートストラップ信頼区間の性質をまとめると、次の通りです。
- 求める信頼区間の区間幅が大きいほど推定精度が悪くなる
- サンプルサイズが多いほど推定精度は良くなる
- ブートストラップ反復回数はごく小さい場合を除き、推定精度には関係しない
- サンプルサイズが小さいときに、ブートストラップ反復数を増やせば推定精度が上がる、ということは無い
この結果から、ブートストラップ信頼区間を求める際にまず注意すべきは、サンプルサイズが十分あるか、という点だと思います。
今回、分散の信頼区間を求めるにあたって、サンプルサイズは 100 程度は必要だということがわかりました。
私は以前 N = 10 に対してブートストラップ信頼区間を計算したことがあり、これは間違いだったと反省せざるをえません。
さらに、サンプルサイズが少ないときに、ブートストラップ反復数を増やしても意味がないこともわかりました。
サンプルサイズが十分でない場合に、どうしても信頼区間を求めたい場合は、区間幅を小さくしてみるといいかもしれません。
以上です。補足も読んでくださいね。
6. 補足
今回の結果から、B = 100 程度でもサンプルサイズが多ければ信頼区間の精度が良いことがわかりました。
しかし、教科書にはブートストラップ信頼区間を求めるには、B = 1000 〜 2000 は必要と書かれています。
これはどういうことでしょうか? 下の図を見てください。
これは、同じサンプルに対してブートストラップ信頼区間を何度も求めたときの上側値(High)と下側値(Low)のばらつきです。
今回の結果では、B = 100 〜 1000 では、真値が信頼区間に入る確率はほぼ同じでしたが、これを見ると、B = 200 では、同じサンプルに対して求まる信頼区間のばらつきが大きいことが見て取れます。
信頼区間のばらつきが安定するのは B = 1000 〜 2000 の間のようです。
つまり、ブートストラップ信頼区間を安定して求めるためには、B = 1000 〜 2000 が必要ということになります。
補足2: calcResult() のコード
library(parallel) library(reshape2) contains <- function(ci, real = 1) { ci[1] < real & real < ci[2] } calcResult <- function(params, repNum = 1000) { cl <- makePSOCKcluster(rep("localhost", detectCores())) clusterSetRNGStream(cl, 314) clusterExport(cl, varlist = c("params", "calcTheoreticalCI", "calcVar", "contains")) result <- sapply(seq_len(nrow(params)), function(i) { success <- parSapply(cl, rep(i, repNum), function(i) { conf <- params[i, "conf"] N <- params[i, "N"] B <- params[i, "B"] x <- rnorm(N) library(simpleboot) boot.out <- one.boot(x, calcVar, B) boot.ci <- boot.ci(boot.out, type = c("norm", "basic", "perc", "bca"), conf = conf) theo.ci <- calcTheoreticalCI(x, conf) c(contains(boot.ci$normal[2:3], 1), contains(boot.ci$basic[4:5], 1), contains(boot.ci$percent[4:5], 1), contains(boot.ci$bca[4:5], 1), contains(theo.ci, 1)) }) conf <- params[i, "conf"] c(conf, apply(success, 1, sum)/repNum) }) stopCluster(cl) result <- data.frame(t(result)) colnames(result) <- c("expect", "normal", "basic", "percentile", "BCa", "theoretical") melt(result, value.name = "prob", variable.name = "group") }