ブートストラップ法 は、標本サイズが大きい場合、 バイアス補正や信頼区間の近似精度が十分良いことが知られています。 しかし、標本サイズが小さい場合、 統計量分散の過小評価や、台の離散性及び外れ値の影響の拡大などの問題が起こります。
生物統計などの分野では十分な標本サイズが確保できないケースは多々あり、 統計量分散を過小評価してしまうと、 特に毒性試験では リスクの観点からとても危険です。 そこで今回は、 ブートストラップ法はどのくらい統計量分散を過小評価してしまうのか 数値実験してみます。
実験は以下の手順で行います。
この 3. により得られた 「ブートストラップ法による『統計量の分散』の平均」が、 理論値よりも下回っているかを確認したいと思います。 「『統計量の分散』の平均」…コンランしそうですね。
また、1. の各標本から統計量を計算し、その N 組の統計量の分散を「実験値」とします。
今回、標本サイズは \(n = \{10, 20, 50, 100, 200\}\), 標本数は \(N=1000\), ブートストラップの反復回数は \(B = \{50, 100, 1000, 2000\}\)として、実験します。 また、母集団分布は標準正規分布 \(N(0,1)\)です。
母集団が標準正規分布のとき、標本平均の分散の理論値はどうなるかというと、 \[ \begin{align} V{(\frac{1}{n} \sum_{i=1}^{n} x_i)} =& \frac{1}{n^2} V{(\sum_{i=1}^{n} x_i)}\\ =& \frac{1}{n} \end{align} \]
なので、この \(\frac{1}{n}\) よりも ブートストラップ標本平均の分散の 平均 が下回るかを確認します。
実験結果は次のようになりました。
| n | B = 50 | B = 100 | B = 1000 | B = 2000 | 実験値 | 理論値 |
|---|---|---|---|---|---|---|
| 10 | 0.09104 | 0.09143 | 0.09260 | 0.09255 | 0.10373 | 0.1 |
| 20 | 0.04673 | 0.04668 | 0.04737 | 0.04741 | 0.05001 | 0.05 |
| 50 | 0.01917 | 0.01937 | 0.01975 | 0.01970 | 0.02010 | 0.02 |
| 100 | 0.00984 | 0.00983 | 0.00993 | 0.00992 | 0.01062 | 0.01 |
| 200 | 0.00497 | 0.00499 | 0.00492 | 0.00483 | 0.00495 | 0.005 |
当然ですが、実験値は理論値になかなか近いですね。
しかし、ブートストラップ法による値は、おおむね理論値よりも小さな値をとっています。 特に標本サイズ n が小さいほど過小評価は顕著ですね。 ブートストラップ反復回数は多いほど過小評価が抑えられているようですが、 \(B=1000\) と \(B=2000\) にはあまり差はないようです。 反復回数を増やせば良いというものではなさそうです。
次は、分散の分散です。
母集団が標準正規分布のとき、不偏でない標本分散の分散の理論値は、 導出は省略しますが、
\[ \begin{align} V{(\frac{1}{n} \sum_{i=1}^{n} (x_i - \mu)^2)} =& \frac{2(n-1)}{n^2} \end{align} \]
となります。
この理論値よりもブートストラップ標本分散の分散の平均が下回るかを確認します。
実験結果は次のようになりました。
| n | B = 50 | B = 100 | B = 1000 | B = 2000 | 実験値 | 理論値 |
|---|---|---|---|---|---|---|
| 10 | 0.13081 | 0.13111 | 0.13220 | 0.13256 | 0.17812 | 0.18 |
| 20 | 0.07775 | 0.07883 | 0.07929 | 0.07932 | 0.09348 | 0.095 |
| 50 | 0.03583 | 0.03606 | 0.03661 | 0.03661 | 0.03804 | 0.0392 |
| 100 | 0.01894 | 0.01914 | 0.01931 | 0.01932 | 0.01961 | 0.0198 |
| 200 | 0.00953 | 0.00974 | 0.00978 | 0.00980 | 0.00925 | 0.00995 |
こちらも同様に、 ブートストラップ標本分散の分散は理論値よりも過小評価してしまっていることがわかります。 また、標本サイズ n が小さいほど過小評価が顕著なこと、 \(B=1000\) と \(B=2000\) にはあまり差はないことも同じですね。
ただ、理論値との差を見ると、 「平均の分散」よりも「分散の分散」の方が、 差が大きいように思えます。
例えば、\((n, B) = {10, 50}\)では、 「平均の分散」は理論値との差が約 \(0.01\) なのに対し、 「分散の分散」は理論値との差が約 \(0.05\) もあります。 分散は 2 次中心モーメントですから、 「高次モーメントの分散」になるほど 過小評価の度合いが大きくなるのでしょうか？
今回の結果は以下のコードから得ました。 無駄に R6 を使っているのは、使ってみたかったからです。 ちゃんと書けた気はしません。
library('R6')
BootstrapVarianceTest <- R6Class(
# ブートストラップ法による統計量の分散を調べるクラス
# 1,
# 母集団分布からサイズnの標本をN個発生させて、
# 「統計量の分散」を求める。
# 2,
# N個の標本それぞれからブートストラップ標本を
# b個ずつ作り「統計量の分散」を求め、
# N個の「統計量の分散」の平均を調べる。
# Args:
# dist: 母集団分布
# n: 標本サイズ
# N: 標本の数
# stat: 統計量
# b: 反復回数
'BootstrapVarianceTest',
private = list(
get_sample.var = function(x){
n <- length(x)
var(x) * (n - 1) / n
},
get_bootstrap.stat = function(x, b, stat){
# xからb個のブートストラップ標本を作り統計量を求める
replicate(b, {
x.boot <- sample(x, length(x), replace=TRUE)
stat(x.boot)
})
},
get_bootstrap.stats = function(xn, b, stat){
# xnのそれぞれの標本から得たブートストラップ標本の統計量を求める
apply(xn, 2, function(x){private$get_bootstrap.stat(x, b, stat)})
},
get_bootstrap.stat.var = function(xn, b, stat){
# ブートストラップ標本の統計量の分散を求める
apply(private$get_bootstrap.stats(xn, b, stat), 2, private$get_sample.var)
}
),
public = list(
xn = NA,
xn.stat.var = NA,
boot.stat.var.mean = NA,
initialize = function(dist, n, N){
self$set_sample(dist, n, N)
},
set_sample = function(dist, n, N){
# サイズnの標本をN個発生させる
self$xn <- replicate(N, dist(n))
},
set_sample.stat.var = function(stat){
# 標本の統計量の分散を求める
self$xn.stat.var <- private$get_sample.var(apply(self$xn,2,stat))
},
get_sample.stat.var = function(){
self$xn.stat.var
},
set_bootstrap.stat.var.mean = function(b, stat){
# ブートストラップ標本の統計量の分散の期待値を求める
self$boot.stat.var.mean <- mean(private$get_bootstrap.stat.var(self$xn, b, stat))
},
get_bootstrap.stat.var.mean = function(){
self$boot.stat.var.mean
}
)
)
###############################################
# Example
###############################################
N <- 1000
n10 <- BootstrapVarianceTest$new(rnorm, n=10, N=N)
n10$set_sample.stat.var(mean)
round(n10$get_sample.stat.var(),digits = 5)
li <- c(50,100,1000,2000)
lapply(li,function(x){n10$set_bootstrap.stat.var.mean(x, mean)})