ブートストラップ法により、パラメータの不偏推定量を求めることができます。
一般に、求めたいパラメータを \(\theta\)、経験分布により求められる推定値を \(\hat{\theta}\)、ブートストラップサンプルにより求められる推定値を \(\hat{\theta}^*\) とおくと、バイアスは \[
\rm{Bias} = E(\hat{\theta}^*) - \hat{\theta}
\] 不偏推定量 \(\hat{\theta_c}\) は \[
\hat{\theta_c} = \hat{\theta} - \rm{Bias}
\] で求めることができます。
不偏分散は良く知られているように、 \[ \sigma^2_c = \frac{1}{n - 1} \sum_{i=1}^n (x_i - \mu)^2 \] で求められるため、ブートストラップを使う必要はありませんが、今回はブートストラップ法で不偏推定値が求められることを確認するために、あえてやってみます。
まず、適当にデータを用意します。
標準正規分布からのサンプルにしましょう。
x <- rnorm(10)
n <- length(x)経験分布による分散の推定値は、標本分散と一致します。(参照: http://rpubs.com/hoxo_m/63120 ) \[
\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \mu)^2
\] R の var() は不偏分散を算出するため、まどろっこしいですが、次のように計算します。
s2.hat <- var(x) * (n - 1) / n
s2.hat## [1] 0.5728373ブートストラップにより求められる推定値は
s2.boot <- replicate(10000, {
x.boot <- sample(x, length(x), replace=TRUE)
var(x.boot) * (n - 1) / n
})
head(s2.boot)## [1] 0.1788716 0.6226761 0.3832046 0.2587838 0.4661823 0.4611383となり、バイアスは
bias <- mean(s2.boot) - s2.hat
bias## [1] -0.0576353と計算されます。
したがって、ブートストラップによる不偏分散は
s2.unbias <- s2.hat - bias
s2.unbias## [1] 0.6304726と求まりました。
これは、不偏分散の理論値に近くなります。
var(x)## [1] 0.6364859上記で求められた推定値が不偏推定値になっていることを確認します。すなわち
\[
E(\hat{\sigma}^2_c) = \sigma^2
\] となることを確認します。
s2.unbias <- replicate(500, {
x <- rnorm(10)
n <- length(x)
s2.hat <- var(x) * (n - 1) / n
s2.boot <- replicate(500, {
x.boot <- sample(x, length(x), replace=TRUE)
var(x.boot) * (n - 1) / n
})
bias <- mean(s2.boot) - s2.hat
s2.unbias <- s2.hat - bias
s2.unbias
})
mean(s2.unbias)## [1] 0.9974269だいたい合っているんじゃないでしょうか。
以上