ブートストラップ

数理統計　第11章 §11.2　★ 母集団分布を経験分布で置き換え、推定量の標準誤差・偏り・信頼区間を見積もる。標本平均や中央値なら経験分布から厳密に計算できる

経験分布とプラグイン原理

$X_1,\dots,X_n$ が分布関数 $F$ から独立に得られたとする。実現値を $x_1,\dots,x_n$ とする。推定量 $\hat\theta=t(X_1,\dots,X_n)$ の標準誤差 $\sqrt{V_F[\hat\theta]}$ は母集団の分布に依存する。そのため $F$ が未知なら計算できない。

標本平均なら標準誤差は $\sigma/\sqrt n$ なので、$\sigma$ を推定すればよい。$\hat\theta$ が中央値や相関係数のような複雑な統計量だと、標準誤差の式自体が手に入らない。ブートストラップ法は、この $F$ を次の経験分布で置き換える。

各 $x$ で $n\hat F_n(x)\sim\mathrm{Bin}(n,F(x))$ である。よって $\hat F_n(x)$ は $F(x)$ の不偏推定量である。大数の法則により、これは $F(x)$ に確率収束する。実は一様収束 $\sup_x\lvert \hat F_n(x)-F(x)\rvert \to0$ も確率1で成り立つ（Glivenko–Cantelli の定理）。

母数が分布の汎関数 $\theta=T(F)$ で表されるとする。その分布を経験分布関数に置き換えたもの $T(\hat F_n)$ をプラグイン推定量と呼ぶ。例を2つ挙げる。

• 平均 $T(F)=\int x\,\mathrm{d} F(x)$ のプラグイン推定量は $\bar x$ である。
• 分散 $T(F)=\int(x-\mu)^2\,\mathrm{d} F(x)$ のプラグイン推定量は

$$\hat\sigma^2=\frac1n\sum_{i=1}^n(x_i-\bar x)^2\tag{11.7}$$

である。除数が $n$ なので、不偏分散 $S^2$ とは違う。両者には $\hat\sigma^2=\frac{n-1}nS^2$ の関係がある。

ブートストラップ推定

推定量 $\hat\theta$ の分散 $V_F[\hat\theta]$ も、分布の汎関数である。そのプラグイン推定量 $V_{\hat F_n}[\hat\theta]$ が作れる。これがブートストラップ推定である。

偏りの定義 $E_F[\hat\theta]-T(F)$ に現れる分布を、全部 $\hat F_n$ に置き換えたものが (11.8) の第2式である。$\hat\theta=T(\hat F_n)$ のとき、偏り修正推定量は $2\hat\theta-E_\ast[\hat\theta^\ast]$ になる。

ブートストラップ分布は $n^n$ 通りの（順序を区別した）復元抽出に等確率 $n^{-n}$ を置いた離散分布なので、原理的には厳密に計算できる。$n$ が大きいと数えきれないので、実際には次のモンテカルロ近似を使う。

1. 各 $b=1,\dots,B$ について、観測値から復元抽出でブートストラップ標本を作り、$\hat\theta^{\ast b}$ を計算する。
2. $\hat\theta^{\ast1},\dots,\hat\theta^{\ast B}$ の標本標準偏差で、$\widehat{\mathrm{se}}_{\mathrm B}$ を近似する。その標本平均と $\hat\theta$ の差で、偏りを近似する。

この手続きには誤差が2つ重なっている。

• 統計的な誤差：$F$ を $\hat F_n$ で置き換えたことによる。データを増やさない限り減らない。
• モンテカルロ誤差：$V_\ast$ を $B$ 個の平均で近似したことによる。$B$ を増やせば、§11.1 の議論どおり $1/\sqrt B$ で減る。

標準誤差の推定なら $B$ は数百、信頼区間の端点なら数千とすることが多い。

標本平均と標本分散での厳密計算

$\hat\theta$ が標本の積率で書けるとき、ブートストラップ分布の平均と分散は経験分布の積率から手で計算できる。1級で問われるとすればこの形である。ブートストラップ標本の各 $X_i^\ast$ は独立に経験分布に従うので

$$E_\ast[X_i^\ast]=\frac1n\sum_{j=1}^nx_j=\bar x,\qquad V_\ast[X_i^\ast]=\frac1n\sum_{j=1}^n(x_j-\bar x)^2=\hat\sigma^2\tag{11.9}$$

である。

証明

前の2つは (11.9) と独立性から、$E_\ast[\bar X^\ast]=\bar x$、$V_\ast[\bar X^\ast]=\frac1{n^2}\cdot n\hat\sigma^2=\hat\sigma^2/n$ である。第3式は $\hat\sigma^{\ast2}=\frac1n\sum X_i^{\ast2}-\bar X^{\ast2}$ と書き、各項の期待値をとる。$E_\ast[X_i^{\ast2}]=V_\ast[X_i^\ast]+\bar x^2=\hat\sigma^2+\bar x^2$、$E_\ast[\bar X^{\ast2}]=V_\ast[\bar X^\ast]+\bar x^2=\hat\sigma^2/n+\bar x^2$ なので

$$E_\ast[\hat\sigma^{\ast2}]=(\hat\sigma^2+\bar x^2)-\Bigl(\frac{\hat\sigma^2}n+\bar x^2\Bigr)=\frac{n-1}n\hat\sigma^2 .$$

偏りの推定は $E_\ast[\hat\sigma^{\ast2}]-\hat\sigma^2=-\hat\sigma^2/n$ で、偏り修正推定量は $\hat\sigma^2-(-\hat\sigma^2/n)=\hat\sigma^2(1+1/n)$ である。

真の偏りは $E_F[\hat\sigma^2]-\sigma^2=-\sigma^2/n$ である。ブートストラップ推定は、その $\sigma^2$ を $\hat\sigma^2$ に置き換えた形になっている。偏り修正推定量 $\hat\sigma^2(n+1)/n$ は、不偏分散 $S^2=\hat\sigma^2\cdot n/(n-1)$ と一致はしない。しかし $n/(n-1)=1+1/n+1/n^2+\cdots$ なので、差は $\mathcal{O}(n^{-2})$ である。

一般の奇数 $n=2m-1$ でも同じ考え方が使える。ブートストラップ標本の中央値が $x_{(j)}$ 以下になる条件を考える。それは、$x_{(1)},\dots,x_{(j)}$ が選ばれた回数の合計が $m$ 以上になることである。この合計は $\mathrm{Bin}(n,j/n)$ に従う。したがって

$$P(\hat\theta^\ast\le x_{(j)})=P(\mathrm{Bin}(n,j/n)\ge m)$$

で、ブートストラップ分布は二項分布の和で厳密に書ける。

ジャックナイフ

ブートストラップより前から使われている偏りの推定法にジャックナイフがある。$i$ 番目の観測値を除いた残りから計算した推定量を $\hat\theta_{(i)}$ と書く。その平均を $\hat\theta_{(\cdot)}=\frac1n\sum_i\hat\theta_{(i)}$ と書き、

$$\widehat{\mathrm{bias}}_{\mathrm J}=(n-1)(\hat\theta_{(\cdot)}-\hat\theta),\qquad \hat\theta_{\mathrm J}=\hat\theta-\widehat{\mathrm{bias}}_{\mathrm J}=n\hat\theta-(n-1)\hat\theta_{(\cdot)}\tag{11.11}$$

を偏りの推定と偏り修正推定量とする。係数 $n-1$ の理由は次のとおりである。偏りが次の左の形なら、$n-1$ 個から作る推定量は右の形になる。

$$E[\hat\theta]-\theta=\frac an+\mathcal{O}(n^{-2}),\qquad E[\hat\theta_{(\cdot)}]-\theta=\frac a{n-1}+\mathcal{O}(n^{-2})$$

主要項の差は

$$\frac a{n-1}-\frac an=\frac a{n(n-1)}$$

である。これを $(n-1)$ 倍すると $a/n$ に戻る。したがって $\hat\theta_{\mathrm J}$ の偏りは $\mathcal{O}(n^{-2})$ になる。

証明

1つの観測値 $x_i$ を除いた残りについて、平均とプラグイン分散を次のように書く（除数は残りの個数 $n-1$）。

$$\bar x_{(i)}=\frac{n\bar x-x_i}{n-1},\qquad \hat\sigma^2_{(i)}=\frac1{n-1}\sum_{j\ne i}(x_j-\bar x_{(i)})^2$$

また $Q=\sum_j(x_j-\bar x)^2=n\hat\sigma^2$ とおく。$x_i$ を除いた偏差平方和について、次が成り立つ。

$$\sum_{j\ne i}(x_j-\bar x_{(i)})^2=Q-\frac{n}{n-1}(x_i-\bar x)^2$$

これは、平均を $\bar x$ から $\bar x_{(i)}$ に動かす平行移動の補正から分かる。実際、次の3つを使う。

$$\begin{aligned}\sum_{j\ne i}(x_j-\bar x)^2&=Q-(x_i-\bar x)^2,\\ \bar x_{(i)}-\bar x&=-\frac{x_i-\bar x}{n-1},\\ \sum_{j\ne i}(x_j-\bar x)&=-(x_i-\bar x).\end{aligned}$$

これらを代入すると

$$\sum_{j\ne i}(x_j-\bar x_{(i)})^2=Q-(x_i-\bar x)^2-2(\bar x_{(i)}-\bar x)\sum_{j\ne i}(x_j-\bar x)+(n-1)(\bar x_{(i)}-\bar x)^2 =Q-(x_i-\bar x)^2\Bigl(1+\frac2{n-1}-\frac1{n-1}\Bigr)$$

となり、括弧は $n/(n-1)$ である。$i$ について平均すると $\sum_i(x_i-\bar x)^2=Q$ より

$$\hat\theta_{(\cdot)}=\frac1{n-1}\Bigl(Q-\frac{Q}{n-1}\Bigr)=\frac{(n-2)Q}{(n-1)^2},\qquad \hat\theta_{\mathrm J}=n\cdot\frac Qn-(n-1)\frac{(n-2)Q}{(n-1)^2}=Q\Bigl(1-\frac{n-2}{n-1}\Bigr)=\frac{Q}{n-1}=S^2 .$$

偏りがちょうど $-\sigma^2/n$ で $\mathcal{O}(n^{-2})$ の項を持たないので、ジャックナイフは偏りを完全に除き、不偏分散に一致する。

ブートストラップ信頼区間

ブートストラップ分布は、信頼区間の構成にも使える。ブートストラップ分布の下側 $100\gamma\%$ 点を $\hat\theta^\ast_{(\gamma)}$ と書く。

• パーセンタイル区間：$[\hat\theta^\ast_{(\alpha/2)},\ \hat\theta^\ast_{(1-\alpha/2)}]$。この区間は単調変換に対して不変である。変換後の母数 $\phi=g(\theta)$ のパーセンタイル区間は、もとの区間の端点を $g$ で写したものと一致する。
• 基本区間：$\hat\theta-\theta$ の分布を $\hat\theta^\ast-\hat\theta$ の分布で近似する。区間は

$$[2\hat\theta-\hat\theta^\ast_{(1-\alpha/2)},\ 2\hat\theta-\hat\theta^\ast_{(\alpha/2)}]$$

である。第8章のピボットによる区間の作り方を、ピボットの分布をブートストラップで近似して行うものである。

• ブートストラップ $t$ 区間：次の2つの量について、前者の分布を後者の分布で近似する。

$$T=\frac{\hat\theta-\theta}{\widehat{\mathrm{se}}},\qquad T^\ast=\frac{\hat\theta^\ast-\hat\theta}{\widehat{\mathrm{se}}^\ast}$$

その分布の下側 $100\gamma\%$ 点を $t^\ast_{(\gamma)}$ とし、区間を

$$[\hat\theta-t^\ast_{(1-\alpha/2)}\widehat{\mathrm{se}},\ \hat\theta-t^\ast_{(\alpha/2)}\widehat{\mathrm{se}}]$$

とする。$T$ が近似的にピボット（分布が $\theta$ によらない）であれば、パーセンタイル区間や基本区間より被覆確率の誤差が小さくなる。

いずれも $\alpha/2$ 点の近似には $B$ を数千とる。

母集団分布が、母数 $\theta$ を持つ分布族 $F_\theta$ に属すると分かっていることもある。そのときは、経験分布の代わりに推定した$F_{\hat\theta}$ からブートストラップ標本を作ってもよい。これをパラメトリック・ブートストラップと呼び、§11.2 の最初に書いた方法（経験分布から作る方法）をノンパラメトリック・ブートストラップと呼んで区別する。

確認テスト

オリジナル

1. 11.2-1 相異なる $n$ 個の観測値から作るブートストラップ標本で、順序を区別しない（多重集合として）異なるものは何通りあるか。$n=5$ のときの数を求めよ。また、ブートストラップ標本の $n$ 個がすべて異なる観測値になる確率を $n=5$ について求めよ。
2. 11.2-2 特定の観測値 $x_1$ がブートストラップ標本に1回も現れない確率を求めよ。$n=10$ のときの値と $n\to\infty$ の極限も求めよ。また、ブートストラップ標本に現れる異なる観測値の個数の期待値を求めよ。
3. 11.2-3 $\mu^2$ の推定量 $\bar X^2$ を考える。偏りのブートストラップ推定と偏り修正推定量を、$\bar x$ と $\hat\sigma^2$ で表せ。それを $\mu^2$ の不偏推定量 $\bar X^2-S^2/n$ と比べよ。
4. 11.2-4 観測値が 0 か 1 で、1 の割合が $\hat p$ であるとする。$\hat p^\ast$ のブートストラップ分布を求めよ。また、標準誤差のブートストラップ推定が$\sqrt{\hat p(1-\hat p)/n}$ であることを示せ。$\mathrm{Bin}(1,\hat p)$ からのパラメトリック・ブートストラップと結果が一致する理由を述べよ。