ブートストラップ法で信頼区間を求める方法はいくつかあります。
今回はこれらのうち、ブートストラップ-t 法を除いた4つの手法で、分散の推定値を求め、その性質を調べてみます。
まず、適当なサンプルを用意します。標準正規分布のサンプルにします。
x <- rnorm(10)
経験分布による分散の推定値を求める関数を定義します。
calcVar <- function(x) {
n <- length(x)
var(x) * (n - 1)/n
}
これを用いると、95%信頼区間は
var.boot <- replicate(10000, {
x.boot <- sample(x, length(x), replace = TRUE)
calcVar(x.boot)
})
# 正規分布近似法
norm.ci <- 2 * calcVar(x) - mean(var.boot) + qnorm(c(0.025, 0.975), sd = sd(var.boot))
# ベーシック法
basic.ci <- 2 * calcVar(x) - quantile(var.boot, probs = c(0.975, 0.025), names = FALSE)
# パーセンタイル法
perc.ci <- quantile(var.boot, probs = c(0.025, 0.975), names = FALSE)
ci <- matrix(c(norm.ci, basic.ci, perc.ci), ncol = 2, byrow = TRUE)
colnames(ci) <- c("Low", "High")
rownames(ci) <- c("Normal", "Basic", "Percentile")
ci
## Low High
## Normal 0.2662 0.9948
## Basic 0.2324 0.9538
## Percentile 0.1918 0.9133
のようにして求めることができます。(BCa 法はややこしいので割愛します)
ところで、R にはブートストラップを行うパッケージがいくつかあります。
simpleboot パッケージを使えば上と同じ計算が簡単にできます。
library(simpleboot)
boot.out <- one.boot(x, calcVar, 10000)
boot.ci(boot.out, conf = 0.95, type = c("norm", "basic", "perc", "bca"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot.out, conf = 0.95, type = c("norm", "basic",
## "perc", "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.2663, 0.9899 ) ( 0.2364, 0.9477 )
##
## Level Percentile BCa
## 95% ( 0.1979, 0.9092 ) ( 0.3021, 1.0908 )
## Calculations and Intervals on Original Scale
また、分散の信頼区間の理論値は下記で求められます。（参照）
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.9)
## [1] 0.3386 1.7228
理論値と比べるとなんだか結構違いますね。
ブートストラップ法により求められる信頼区間の性質を調べてみます。
まず、信頼区間は分布の両端を求めなければならず、端に行くほど精度が落ちるはずです。
求める信頼区間を 10%, 20%, …, 90% と変更してみて、信頼区間にどれだけ真値 (=1) が入るか確認してみます。
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")
}
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 は理論値です。
思ったとおり、端に行くほどブートストラップ信頼区間の精度は落ちるようです。
次に、サンプル数 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 が大きいほど精度は良くなるようです。
これは、サンプル数が多いほど、経験分布は確率分布に近づくためだと考えられます。
次に、ブートストラップ反復回数 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 で試してみます。
B = 2 ～ 40 では、B によってブートストラップ信頼区間の精度が向上しているのがわかります。
(BCa 法は B が小さいと求まらないので除外しています)
ところで、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)
経験分布と確率分布の違いが大きい場合はブートストラップでどんなに頑張っても信頼区間を求めるのは難しそうです。
ブートラップ信頼区間の性質を調べてみました。
今回わかったのは、ブートストラップ法で信頼区間を求めるには、前提として、ある程度のサンプル数が必要であり、ブートストラップ反復数 B をどんなに大きくしても、経験分布と確率分布がずれていれば意味がない、ということですね。
今回の結果から、B = 100 程度でもサンプル数が多ければ信頼区間の精度が良いことがわかりました。
しかし、教科書にはブートストラップ信頼区間を求めるには、B = 1000 ～ 2000 は必要と書かれています。
これはどういうことでしょうか？
下の図を見てください。
これは、同じサンプルに対してブートストラップ信頼区間を何度も求めたときの上側値(High)と下側値(Low)のばらつきです。
今回の結果では、B = 100 ～ 1000 では、真値が信頼区間に入る確率はほぼ同じでしたが、これを見ると、B = 200 では、同じサンプルに対して求まる信頼区間のばらつきが大きいことが見て取れます。
信頼区間のばらつきが安定するのは B = 1000 ～ 2000 の間のようです。
つまり、ブートストラップ信頼区間を安定して求めるためには、B = 1000 ～ 2000 が必要ということになります。