手元に 40 人分の年収データがあり、その 中央値 を報告することになったとします。「中央値は 420 万円でした」だけでは足りません。査読者や上司は必ずこう聞きます — 「その 420 万円は、どのくらいブレる数字なの？」。
標本平均であれば答えは簡単です。標準誤差は $s/\sqrt{n}$ という公式があり、電卓で 1 行です。ところが 中央値の標準誤差には、教科書に載っているような閉じた公式がありません。分位点も、変動係数も、四分位範囲も、2 群の相関係数の差も同じです。統計学の教科書に出てくる「きれいな公式」は、ごく一部の統計量にしか用意されていないのです。
この壁を、たった 1 つのアイデアで壊したのが ブートストラップ法（bootstrap method） です。1979 年に Bradley Efron が提案したこの手法は、統計量が何であろうと、分布が何であろうと、まったく同じ手続きで標準誤差と信頼区間を出します。手元のデータからサイコロを振り直すように再標本を作り、そのばらつきを見るだけです。
図の左側が「本当はやりたいけれど、できないこと」です。母集団から標本を何度も取り直せれば、統計量がどれだけばらつくかは実際に見れば分かります。しかし現実に手に入るのは 1 組の標本だけで、この経路は最初から閉ざされています。右側がブートストラップの発想で、母集団の代わりに手元の標本そのものを「母集団の代役」として使い、そこから復元抽出を繰り返します。左の破線の矢印が右の実線の矢印に置き換わっている、この 1 点がブートストラップのすべてです。
ブートストラップが効いてくる場面は統計学の内側だけではありません。機械学習の バギングやランダムフォレストの OOB 推定はブートストラップそのものですし、A/B テストで「コンバージョン率の比」や「95 パーセンタイル応答時間」の信頼区間を出すとき、実務で使われているのはたいていブートストラップです。センサ計測の分野でも、非線形な校正式を通した後の測定量の不確かさ評価に使われます。
本記事では次の順で見ていきます。
- 復元抽出が実際に何をしているか（1 回の再標本に元データの 63.2% しか入らない理由）
- プラグイン原理 — ブートストラップが「何を何で置き換えているのか」
- 公式が存在しない統計量への適用（本手法の最大の価値）
- 信頼区間の作り方 4 種と、その被覆確率の実測比較
- 反復回数 $B$ の決め方、バイアス補正、ジャックナイフとの違い
- ブートストラップが失敗する場面 — 極値統計・極小標本・時系列
- 機械学習のバギング／OOB とのつながり
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ブートストラップ法とは — 標本を母集団の代わりに使う
統計的推論が本当に欲しがっているもの
推定値そのもの（$\bar{x} = 3.29$、中央値 $= 2.77$）は、データを見ればすぐ計算できます。統計的推論が本当に欲しいのはその先、「同じ母集団からもう一度データを取り直したら、この推定値はどれくらい違う値になるか」 という情報です。これを表すのが 統計量の標本分布 で、その標準偏差が標準誤差、そのパーセンタイルが信頼区間の材料になります。
標本平均については、この標本分布が理論的に分かっています。中心極限定理により $\bar{X}$ は $n$ が大きければ近似的に正規分布に従い、その標準偏差は $\sigma/\sqrt{n}$ です。$\sigma$ を標本標準偏差 $s$ で置き換えれば標準誤差 $s/\sqrt{n}$ が得られます。
ところがここから先が問題です。中央値の標本分布は？ 90 パーセンタイルは？ 「上位 10% の平均と下位 10% の平均の比」は？ これらを解析的に導こうとすると、母集団の密度関数 $f$ の中央値における値 $f(m)$ のようなそもそも観測できない量が式に現れてしまい、実務では手が止まります。
「もう一度取り直せたら」を強引に実現する
ブートストラップの核心は、次の一文に尽きます。
母集団 $F$ が分からないなら、$F$ の最良の推定値である経験分布 $\hat{F}_n$ を母集団の代役にすればいい。
経験分布 $\hat{F}_n$ とは、観測された $n$ 個の値それぞれに $1/n$ の確率を置いた離散分布のことです。そして「$\hat{F}_n$ からサイズ $n$ の標本を取る」という操作は、実際にやってみると 元のデータから $n$ 個を復元抽出する ことにほかなりません。母集団からのサンプリングという「できないこと」が、手元のデータからのリサンプリングという「いくらでもできること」にすり替わったわけです。
アルゴリズム
元の標本を $\bm{x} = (x_1, x_2, \dots, x_n)$、関心のある統計量を $\hat{\theta} = s(\bm{x})$ とします。
- $\bm{x}$ から大きさ $n$ の 復元抽出（with replacement） を行い、ブートストラップ標本 $\bm{x}^{*b} = (x_1^{*b}, x_2^{*b}, \dots, x_n^{*b})$ を得る
- ブートストラップ標本に対して統計量を計算: $\hat{\theta}^{*b} = s(\bm{x}^{*b})$
- ステップ 1〜2 を $B$ 回繰り返す（$b = 1, 2, \dots, B$）
- 得られた $\hat{\theta}^{*1}, \hat{\theta}^{*2}, \dots, \hat{\theta}^{*B}$ の分布が ブートストラップ分布
注目すべきは、この手続きのどこにも「統計量が何であるか」が入っていないことです。ステップ 2 の $s(\cdot)$ を平均から中央値に差し替えても、相関係数に差し替えても、機械学習モデルの交差検証スコアに差し替えても、手続きはまったく変わりません。これがブートストラップの汎用性の正体です。
ただし、ここまでは「復元抽出をする」と言葉で言っているだけです。次節では、その復元抽出が実際にどんなデータの集まりを作っているのかを、具体的に覗いてみましょう。
復元抽出は何をしているのか — 63.2% という数字
1 回の再標本の中身を見る
「$n$ 個から $n$ 個を復元抽出する」と聞くと、結局ほとんど同じデータが返ってくるように思えるかもしれません。実際には、ある点は 3 回選ばれ、別の点は 1 回も選ばれない、というかなり偏った標本になります。
上段のヒートマップは、$n = 20$ のデータから 4 回の復元抽出を行い、各データ点が何回選ばれたかを数えたものです。0 のマスが「その回では一度も使われなかったデータ点」で、かなりの数あることが分かります。左下の棒グラフのとおり、この 4 回の再標本に含まれた元データの割合は 0.60 / 0.70 / 0.60 / 0.70 で、赤い破線 $1 – 1/e = 0.6321$ の周りに散らばっています。右下は $n$ を変えて 2000 回ずつ平均を取ったもので、$n = 20$ で理論値 0.6415・実測 0.6393、$n = 100$ で理論 0.6340・実測 0.6345、$n = 1000$ で理論 0.6323・実測 0.6322 と、理論曲線と実測がぴたりと重なっています。
63.2% の導出
なぜ 63.2% なのかは、簡単な計算で出ます。まず、ある特定のデータ点 $x_i$ が「1 回の抽選で選ばれない」確率は、$n$ 個から等確率に 1 個引くので
$$ 1 – \frac{1}{n} $$
です。復元抽出では各抽選が独立なので、$n$ 回の抽選すべてで $x_i$ が選ばれない確率は、これを $n$ 乗して
$$ P(x_i \text{ が再標本に入らない}) = \left(1 – \frac{1}{n}\right)^{n} $$
となります。ここで、指数関数の定義 $\lim_{n\to\infty}(1 + a/n)^n = e^{a}$ に $a = -1$ を当てはめると、右辺は $n \to \infty$ で $e^{-1}$ に収束します。したがって $x_i$ が 入る確率は
$$ 1 – \left(1 – \frac{1}{n}\right)^{n} \;\xrightarrow[n \to \infty]{}\; 1 – \frac{1}{e} = 0.6321\dots $$
各データ点について「入る／入らない」の期待値を足し合わせれば、再標本に含まれる元データの割合の期待値がそのまま $1 – 1/e$ になります。この 63.2%（と、裏返しの 36.8%）は本記事の後半で OOB（out-of-bag）推定として再登場するので、覚えておいてください。
つまりブートストラップ標本は、元データの約 3 分の 2 を使い、残り 3 分の 1 を捨て、使う分は重複させた「少しずつ違う標本」の集まりです。では、この人工的に作った標本たちのばらつきは、本物の標本分布をどれだけ再現しているのでしょうか。
プラグイン原理 — 何を何で置き換えているのか
経験分布関数
$n$ 個のデータ $x_1, \dots, x_n$ の 経験分布関数（empirical distribution function） は、各 $x_i$ に質量 $1/n$ を置いた離散分布です。
$$ \hat{F}_n(t) = \frac{1}{n}\sum_{i=1}^{n}\mathbf{1}(x_i \leq t) $$
これが母集団の代役として使える根拠が、グリベンコ・カンテリの定理です。この定理は $\hat{F}_n$ が真の分布 $F$ に一様収束することを保証します。
$$ \sup_t |\hat{F}_n(t) – F(t)| \xrightarrow{a.s.} 0 \quad (n \to \infty) $$
「一様」というのがポイントで、どこか一点で近いのではなく、分布関数の全域で同時に近づきます。つまり $n$ が十分大きければ、$\hat{F}_n$ は $F$ の代わりとして相当に良い、ということです。
プラグイン原理
母集団の特性が分布関数 $F$ の汎関数 $\theta = T(F)$ で表されるとき、その推定量として経験分布 $\hat{F}_n$ を「プラグイン（差し込む）」した
$$ \hat{\theta} = T(\hat{F}_n) $$
を用いるのが プラグイン原理 です。実は私たちが普段使っている推定量の多くは、これを無意識にやっています。
- 母平均 $\mu = T(F) = \int x \, dF(x)$ → $\hat{\mu} = T(\hat{F}_n) = \frac{1}{n}\sum x_i = \bar{x}$
- 母分散 $\sigma^2 = T(F) = \int (x – \mu)^2 \, dF(x)$ → $\hat{\sigma}^2 = T(\hat{F}_n) = \frac{1}{n}\sum(x_i – \bar{x})^2$
- 母中央値 $m = F^{-1}(1/2)$ → $\hat{m} = \hat{F}_n^{-1}(1/2) = $ 標本中央値
ブートストラップは「分布そのもの」にプラグインする
ブートストラップの独創性は、推定量ではなく標本分布にプラグイン原理を適用した点にあります。統計量 $\hat{\theta} = s(\bm{X})$ の標本分布は $F$ に依存するので、これを $G_F$ と書きます。
$$ G_F(\cdot) = P_F(s(\bm{X}) \leq \cdot) $$
私たちが本当に知りたいのはこの $G_F$ です。ブートストラップは $F$ を $\hat{F}_n$ に置き換えて、
$$ \hat{G}_n(\cdot) = P_{\hat{F}_n}(s(\bm{X}^*) \leq \cdot) $$
を代わりに使います。$\bm{X}^*$ は $\hat{F}_n$ からの復元抽出、すなわち元のデータからのリサンプリングです。
誤差は 2 種類ある — この区別が実務で効く
ここで重要な整理をしておきます。ブートストラップの結果には、性質のまったく違う 2 つの誤差が混ざっています。
- 統計的近似の誤差: $G_F$ を $\hat{G}_n$ で置き換えたことによる誤差。これは $\hat{F}_n \neq F$ であることに由来し、データを増やす以外に減らす方法がありません。
- モンテカルロ誤差: $\hat{G}_n$ を厳密に計算するのは組合せ爆発（$\binom{2n-1}{n}$ 通り）で不可能なので、$B$ 回のリサンプリングで近似します。この誤差は $B$ を増やせばいくらでも小さくできます。
$B = 10^6$ にすれば信頼区間が正確になる、というのは誤解です。$B$ を増やして消えるのは 2 番目だけで、1 番目は $n$ が決めます。この点は後で図でも確認します。
では実際、$\hat{G}_n$ は $G_F$ をどれだけ再現できているのでしょうか。答え合わせができるように、母集団を人工的に決めて実験してみます。
母集団を指数分布（scale = 3）とし、赤い塗りつぶしが「母集団から $n$ 個を 2 万回サンプリングして得た真の標本分布」、色つきの折れ線が「たまたま手に入った 1 組の標本からのブートストラップ分布」です。手元の標本を 3 通り変えて重ねてあります。読み取れるのは 2 点で、位置（中心）は標本ごとに大きくずれるが、幅（ばらつき）はよく再現されているということです。実測では、標本平均の真の SD が $n=20$ で 0.6718 に対しブートストラップ SD の平均が 0.6186（真の 92%）、$n=100$ では 0.3008 に対し 0.2963（真の 98%）でした。中央値でも $n=20$ で 103%、$n=100$ で 101% です。
中心がずれるのは当然で、手元の標本の中心はそもそも真の値からずれているからです。ブートストラップが推定しているのは中心の位置ではなく 中心のまわりのばらつき であって、それは狙いどおり再現されています。なお下段の中央値のグラフがギザギザなのは、$n$ 個のデータから作る中央値が取りうる値が有限個しかないためで、リサンプリング特有の離散性です。
モンテカルロ近似としての実装
$\hat{G}_n$ を解析的に求めるのは一般に困難ですが、$B$ 回のリサンプリングによるモンテカルロシミュレーションで近似できます。前述のとおり、この近似誤差は $B$ を増やせば消えます。
ここまでで「ブートストラップが何をしているか」が固まりました。次はいよいよ、この手続きから標準誤差を取り出します。
ブートストラップ標準誤差
定義
統計量 $\hat{\theta}$ の標準誤差は、その標本分布の標準偏差です。
$$ \text{SE}_F(\hat{\theta}) = \sqrt{\text{Var}_F[s(\bm{X})]} $$
ブートストラップ標準誤差は、この式の $F$ を $\hat{F}_n$ で置き換えたものにほかなりません。
$$ \text{SE}_{\text{boot}}(\hat{\theta}) = \sqrt{\text{Var}_{\hat{F}_n}[s(\bm{X}^*)]} $$
実際には $B$ 個の値の標本標準偏差として計算します。
$$ \hat{\text{SE}}_{\text{boot}} = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\left(\hat{\theta}^{*b} – \bar{\theta}^*\right)^2} $$
ここで $\bar{\theta}^* = \frac{1}{B}\sum_{b=1}^{B}\hat{\theta}^{*b}$ はブートストラップ推定値の平均です。
標本平均で答え合わせをする
新しい手法を信用してよいかは、「答えが分かっている問題」に当てて確かめるのが定石です。$\hat{\theta} = \bar{X}$ の場合、真の標準誤差は $\text{SE} = \sigma/\sqrt{n}$ です。
ブートストラップ標本 $\bm{X}^*$ は $\hat{F}_n$ からのサイズ $n$ の復元抽出です。$\hat{F}_n$ は各 $x_i$ に確率 $1/n$ を置いた分布なので、その分散は定義どおり
$$ \text{Var}(\hat{F}_n) = \frac{1}{n}\sum_{i=1}^{n}(x_i – \bar{x})^2 = \hat{\sigma}^2 $$
です。ここで $\bm{X}^*$ の各成分は $\hat{F}_n$ からの独立同分布な抽出なので、標本平均の分散公式 $\text{Var}(\bar{X}) = \text{Var}(X)/n$ がそのまま使えて、
$$ \text{Var}_{\hat{F}_n}(\bar{X}^*) = \frac{\hat{\sigma}^2}{n} \quad\Longrightarrow\quad \text{SE}_{\text{boot}} = \frac{\hat{\sigma}}{\sqrt{n}} $$
が得られます。つまり真の $\sigma/\sqrt{n}$ を $\hat{\sigma}/\sqrt{n}$ で推定していることになり、教科書どおりの標準誤差の公式が自動的に再現されました。ブートストラップは既存の公式を否定する手法ではなく、公式があるところでは公式に一致する、その意味で安全な一般化なのです。
本領発揮 — 公式が存在しない統計量
答え合わせが済んだので、公式が存在しない領域に踏み込みます。ここがブートストラップの最大の価値です。
対数正規分布から取った $n = 40$ のデータ 1 組に対して、6 つの統計量へまったく同じコードを適用した結果です。左上の標本平均だけは検算ができて、ブートストラップ SE が 0.3266、教科書公式 $s/\sqrt{n}$ が 0.3281 と小数第 2 位まで一致しました。残りの 5 つ — 中央値（SE = 0.4190）、90% 分位点（SE = 0.6813）、変動係数（SE = 0.1075）、四分位範囲（SE = 0.4353）、相関係数（SE = 0.0821）— には、そもそも比較すべき初等的な公式がありません。
読み取ってほしいのは 3 点です。第一に、中央値の SE 0.4190 は平均の SE 0.3266 より 28% 大きいこと。正規分布なら中央値の SE は平均の約 1.25 倍というよく知られた関係がありますが、ブートストラップはその種の知識を一切使わずに同じ結論に到達します。第二に、中央値と 90% 分位点のヒストグラムが櫛状に離散化していること — 分位点は元データの値そのものしか取れないためで、$n$ が小さいときの本質的な制約です。第三に、変動係数の分布が二山になっていること。これは 1 個の大きな外れ値が「再標本に入るか入らないか」で結果が二分されるためで、「この推定値は特定の 1 点に強く依存している」という診断情報がタダで手に入っています。公式では絶対に見えない情報です。
右下の相関係数も示唆的です。正規性を仮定した近似式 $(1-r^2)/\sqrt{n-3}$ は 0.0920 を返しますが、ブートストラップは 0.0821 と、やや小さい値を出しました。データが正規分布から外れている以上、正規性を前提にした公式より、データ自身に語らせたほうが素直です。
標準誤差が出れば、次は信頼区間です。ところがここには落とし穴があり、「ブートストラップ信頼区間」と一口に言っても作り方が何通りもあって、性能がはっきり違います。
ブートストラップ信頼区間 — 4 つの流儀
パーセンタイル法
最も単純な方法です。ブートストラップ分布の $\alpha/2$ 分位点と $1 – \alpha/2$ 分位点をそのまま区間の端点にします。
$$ CI_{\text{percentile}} = \left[\hat{\theta}^*_{(\alpha/2)}, \hat{\theta}^*_{(1-\alpha/2)}\right] $$
95% 信頼区間なら、$B$ 個のブートストラップ推定値を小さい順に並べ、2.5 パーセンタイルと 97.5 パーセンタイルを取るだけです。
この方法には見逃せない長所があります。それは 変換不変性 です。$\theta$ の区間が $[L, U]$ なら、単調変換 $g$ を通した $g(\theta)$ の区間は自動的に $[g(L), g(U)]$ になります。相関係数のように $[-1, 1]$ に収まるべき量や、分散のように非負であるべき量を扱うとき、この性質は「あり得ない区間」を出さないという安心につながります。
基本（pivotal）法
パーセンタイル法とは逆向きの発想です。$\hat{\theta} – \theta$ というずれの分布を、$\hat{\theta}^* – \hat{\theta}$ の分布で近似します。$\hat{\theta} – \theta$ が区間 $[a, b]$ に入る確率が 95% なら、$\theta$ は $[\hat{\theta} – b, \hat{\theta} – a]$ に入る、という理屈で、結果として
$$ CI_{\text{basic}} = \left[2\hat{\theta} – \hat{\theta}^*_{(1-\alpha/2)},\; 2\hat{\theta} – \hat{\theta}^*_{(\alpha/2)}\right] $$
となります。パーセンタイル法の区間を $\hat{\theta}$ を中心に左右反転した形です。ブートストラップ分布が右に歪んでいれば、基本法の区間は左に長くなります。
正規近似法
ブートストラップは標準誤差の計算だけに使い、区間は正規分布で作る折衷案です。
$$ CI_{\text{normal}} = \left[\hat{\theta} – z_{1-\alpha/2}\,\hat{\text{SE}}_{\text{boot}},\; \hat{\theta} + z_{1-\alpha/2}\,\hat{\text{SE}}_{\text{boot}}\right] $$
必ず左右対称になるため、歪んだ分布には向きません。
BCa 法
BCa（Bias-Corrected and accelerated）法は、パーセンタイル法の切り取り位置を 2 つの補正でずらします。
$$ CI_{\text{BCa}} = \left[\hat{\theta}^*_{(\alpha_1)}, \hat{\theta}^*_{(\alpha_2)}\right] $$
$$ \alpha_1 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{\alpha/2}}{1 – \hat{a}(\hat{z}_0 + z_{\alpha/2})}\right), \qquad \alpha_2 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{1-\alpha/2}}{1 – \hat{a}(\hat{z}_0 + z_{1-\alpha/2})}\right) $$
一見複雑ですが、$\hat{z}_0 = 0, \hat{a} = 0$ を代入すると分子・分母が簡単になり $\alpha_1 = \Phi(z_{\alpha/2}) = \alpha/2$、$\alpha_2 = 1 – \alpha/2$ となって、パーセンタイル法にぴたりと戻ります。つまり BCa は「補正が必要ないときは何もしない」拡張です。
バイアス補正係数 $\hat{z}_0$ は、ブートストラップ分布の中心が $\hat{\theta}$ からどれだけずれているかを測ります。
$$ \hat{z}_0 = \Phi^{-1}\left(\frac{\#\{\hat{\theta}^{*b} < \hat{\theta}\}}{B}\right) $$
$\hat{\theta}$ にバイアスがなければ、ちょうど半分の $\hat{\theta}^{*b}$ が $\hat{\theta}$ より小さいので $\Phi^{-1}(0.5) = 0$、すなわち $\hat{z}_0 = 0$ です。
加速度定数 $\hat{a}$ は、標準誤差が $\theta$ の値によって変わる度合い（分布の歪み）を、ジャックナイフ推定値の 3 次モーメントで測ります。
$$ \hat{a} = \frac{\sum_{i=1}^{n}(\bar{\theta}_{(\cdot)} – \hat{\theta}_{(-i)})^3}{6\left[\sum_{i=1}^{n}(\bar{\theta}_{(\cdot)} – \hat{\theta}_{(-i)})^2\right]^{3/2}} $$
ここで $\hat{\theta}_{(-i)}$ は $i$ 番目のデータを除いて計算した ジャックナイフ推定値、$\bar{\theta}_{(\cdot)} = \frac{1}{n}\sum_i \hat{\theta}_{(-i)}$ はその平均です。分子が 3 乗和なので歪度に相当する量になり、対称なら 0 になります。
4 つを同じデータで比べる
言葉で並べても差は分かりません。相関係数を例に、同じブートストラップ分布から 4 つの区間を作ってみます。相関係数を選んだのは、$[-1, 1]$ という明確な範囲があるので「区間が正しいか」が一目で判定できるからです。
真の相関 $\rho = 0.7$ の二変量正規分布から $n = 20$ を取り、標本相関は $\hat{r} = 0.750$ でした。ブートストラップ分布（上段）は上限 1 に押しつぶされて左に長い裾を引いています。下段の結果は決定的で、正規近似法は上限 1.023、基本法は上限 1.154 と、相関係数としてあり得ない値を出しました。一方パーセンタイル法 [0.345, 0.895] と BCa 法 [0.427, 0.912] は範囲内に収まっています。前述の変換不変性が効いた形です。
さらに BCa の中身を見ると、$\hat{z}_0 = 0.007$、$\hat{a} = 0.076$ で、切り取り位置が 2.5%/97.5% から 4.52%/98.98% へずれています。その結果、区間幅は 0.485 とパーセンタイル法の 0.550 より 12% 狭くなりました。歪みを補正した分だけ無駄が削れています。基本法がパーセンタイル法の鏡像（下側 0.146 / 上側 0.404 と、パーセンタイル法の下側 0.404 / 上側 0.146 がちょうど入れ替わっている）になっている点も、定義どおりです。
決着 — 被覆確率を実測する
とはいえ「1 例で範囲を守った」だけでは推奨の根拠になりません。信頼区間の性能を測る正しい物差しは 被覆確率（coverage probability）、すなわち手続きを何度も繰り返したとき真値を含んだ割合が名目の 95% に一致するか、です。
右に強く歪んだ対数正規分布の母平均について、各 $n$ で 1500 回ずつ実験した結果です。まず全体として、すべての手法が名目 95% を下回っています。$n = 10$ ではパーセンタイル法 0.835、基本法 0.819、BCa 法 0.849 と、95% どころか 85% 前後です。$n = 25$ で 0.865〜0.899、$n = 60$ で 0.909〜0.929、$n = 150$ でようやく 0.932〜0.944 まで回復します。歪んだ母集団では、ブートストラップ信頼区間は小標本で確実に楽観的すぎるということです。
手法間の順位ははっきりしていて、BCa 法がすべての $n$ でブートストラップ勢の最良（0.849 / 0.899 / 0.929 / 0.944）、基本法が一貫して最悪（0.819 / 0.865 / 0.909 / 0.932）でした。右のパネルは $n = 25$ での外し方の内訳で、基本法は 12.7% が「真値が区間より上」側の失敗、つまり片側に偏って外していることが分かります。BCa はこの偏りを 7.1% まで減らしています。歪んだ分布では左右対称に外れないので、対称な補正では足りないのです。
なお $n = 10$ では正規性を仮定した t 区間が 0.889 と最も高い被覆でしたが、これは平均幅 4.430 と、BCa の 3.929 より 13% も広いことの結果です。被覆確率は幅とセットで見る必要があります。
実務上の指針は明快です。迷ったら BCa。計算資源が厳しければパーセンタイル法。基本法は積極的に選ぶ理由が乏しい。そして小標本＋強い歪みなら、どの手法も名目水準を信用しすぎない。
ここまで手法の選択を見てきましたが、もう 1 つ決めるべきパラメータが残っています。反復回数 $B$ です。
反復回数 $B$ はいくつにすべきか
$B$ を大きくすれば結果が良くなる、と考えたくなりますが、先ほど整理したとおり $B$ が減らすのはモンテカルロ誤差だけです。
同一のデータ（$n = 40$）に対して乱数シードだけを変え、中央値の 95% 信頼区間を 60 回ずつ作り直した結果です。左のパネルで、$B = 50$ では下限が 2.05〜2.40 の範囲でばらついているのに対し、$B = 10000$ では点がほぼ 1 本の線に収束しています。中央のパネルは端点の標準偏差の対数プロットで、$B = 50$ での 0.1043 から $B = 10000$ での 0.0059 まで、赤い点線で示した $1/\sqrt{B}$ の傾きに沿って減少しています。
決定的なのは右のパネルです。ブートストラップ SE の平均値は $B = 50$ で 0.4261、$B = 10000$ で 0.4204 とほとんど動かず、動いているのはエラーバー（ばらつき）だけです。$B$ を 200 倍にしても推定値そのものは 1% ちょっとしか変わりません。$B$ は「乱数の粗さ」を決めるパラメータであって、データが持つ情報量には触れないのです。
実務的な目安は、標準誤差だけなら $B = 200$ 程度、信頼区間なら分布の裾を推定する必要があるため $B = 1000$〜$10000$、論文に載せる数値なら $B = 10000$ 以上、といったところです。上の実験でも $B = 1000$ で端点のばらつきは 0.0196 と実用上十分小さくなっています。
$B$ の話が片付いたので、ブートストラップのもう 1 つの使い道 — ばらつきではなく偏りの推定 — に進みます。
バイアスの推定と補正
ブートストラップは標準誤差だけでなく、推定量のバイアス（偏り）も推定できます。真のバイアスは $\text{Bias} = E_F[\hat{\theta}] – \theta$ ですが、この式の $F$ を $\hat{F}_n$ に、$\theta$ を $\hat{\theta}$ に置き換えると
$$ \widehat{\text{Bias}}_{\text{boot}} = \bar{\theta}^* – \hat{\theta} $$
という、拍子抜けするほど簡単な推定量が得られます。ブートストラップ推定値の平均が元の推定値からずれていれば、その分だけ偏っているだろう、というわけです。補正した推定量は
$$ \hat{\theta}_{\text{corrected}} = \hat{\theta} – \widehat{\text{Bias}}_{\text{boot}} = 2\hat{\theta} – \bar{\theta}^* $$
となります。答え合わせに使えるのは、$1/n$ で割る版の標本分散 $\hat{\sigma}^2 = \frac{1}{n}\sum(x_i – \bar{x})^2$ です。これは真のバイアスが $-\sigma^2/n$ であることが理論的に分かっています。
$\sigma^2 = 4$ の正規母集団から $n = 25$ を 3000 回取った結果です。左のパネルで、補正前の推定値の平均は 3.826（真値 4.0 より 0.174 小さい）ですが、補正後は 3.979 とほぼ真値に乗りました。中央のパネルは各標本でのバイアス推定値の分布で、平均 $-0.1528$ が理論値 $-\sigma^2/n = -0.1600$ とよく一致しています。標本 1 組だけから「自分がどれだけ偏っているか」を当てられている、という点が面白いところです。
ただし右のパネルが重要な警告になっています。バイアスは 0.174 から 0.021 へほぼ消えましたが、RMSE は 1.1066 から 1.1367 へむしろ悪化しました。バイアス補正は $\bar{\theta}^*$ という推定量を引き算するので、その分だけ分散が増えるからです。参考として不偏分散（$1/(n-1)$ 版）の RMSE は 1.1385 で、これも補正後とほぼ同じです。「偏りをなくす」ことと「誤差を小さくする」ことは別物であり、バイアス補正は前者だけを保証します。バイアスが標準誤差に比べて小さいとき（目安として $|\widehat{\text{Bias}}| < 0.25\,\hat{\text{SE}}$）は補正しないほうが賢明です。
ところで、BCa の加速度定数のところで「ジャックナイフ」という別のリサンプリング手法が登場しました。ブートストラップとどう違うのか、ここで整理しておきます。
ジャックナイフとの関係
ジャックナイフ（jackknife） はブートストラップより古く、1949 年に Quenouille が提案した手法です。復元抽出の代わりに、$i$ 番目のデータを 1 つだけ抜いた標本を $n$ 通り作ります。
$$ \hat{\theta}_{(-i)} = s(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n), \qquad i = 1, \dots, n $$
標準誤差の推定量は
$$ \hat{\text{SE}}_{\text{jack}} = \sqrt{\frac{n-1}{n}\sum_{i=1}^{n}\left(\hat{\theta}_{(-i)} – \bar{\theta}_{(\cdot)}\right)^2} $$
です。乱数を使わず計算量も $n$ 回で済むという利点があり、実際 BCa の $\hat{a}$ はこれを使っています。しかし決定的な弱点があります。
標準正規母集団から各 $n$ につき 600 標本を取り、「推定 SE ÷ 真の SE」の分布を描いたものです。点が平均、ひげが 5〜95% の範囲を表します。左の平均では両者の挙動がほぼ完全に一致し、変動係数（CV）はどちらも $n=15$ の 0.20 から $n=125$ の 0.07 へ、$n$ とともにきれいに縮んでいきます。
決定的な差は右の中央値に現れます。平均値だけ見るとジャックナイフも 0.95〜1.00 で悪くないのですが、ひげの長さがまったく縮みません。ジャックナイフの CV は $n=15$ で 0.65、$n=125$ でも 0.67 と、標本サイズを 8 倍以上にしてもほぼ不変です。5〜95% 幅も 1.94 → 2.08 と、むしろ広がってさえいます。対するブートストラップは CV 0.32 → 0.21、幅 1.11 → 0.71 と着実に縮んでいます。
これは偶然ではなく、ジャックナイフは中央値の標準誤差に対して一致性を持たないという理論的事実の現れです。理由は直感的で、$n$ 個から 1 個抜いても中央値は隣の順序統計量へ「1 段」ずれるだけで、その変化量が $n$ を増やしても比例して小さくならないためです。一方ブートストラップは中央値がもっと大きく動く再標本も作れるので、分布全体を捉えられます。
滑らかな統計量（平均、回帰係数など）にはジャックナイフでも十分だが、非滑らかな統計量（中央値、分位点、最大値など）にはブートストラップを使う — これが実用的な使い分けです。
さて、ここまでは推定の話でした。ブートストラップは検定にも使えます。
ブートストラップ検定
基本的な枠組み
帰無仮説 $H_0: \theta = \theta_0$ を検定するとき、観測された検定統計量を $t_{\text{obs}}$ とすると、ブートストラップ p 値は次のように計算します。
- 帰無仮説のもとでのリサンプリング: データを $H_0$ が成り立つように変換してからリサンプリングする（この変換が肝心です）
- 各ブートストラップ標本から検定統計量 $t^{*b}$ を計算
- $p = \dfrac{\#\{|t^{*b}| \geq |t_{\text{obs}}|\} + 1}{B + 1}$
3 番目の式で分子・分母に 1 を足しているのは、$B$ が有限のときに p 値が 0 になるのを防ぐためで、実務ではこの形が推奨されます。
2 標本問題での実装
2 群の平均が等しいかを検定する場合、帰無仮説のもとでは 2 群のデータは同じ分布から来ているはずです。そこで、
- 2 群のデータをプールする
- プールされたデータから復元抽出で元と同じサイズの 2 群に分け直す
- 2 群の平均差を計算する
これを $B$ 回繰り返して得られる分布が、$H_0$ が正しいときに平均差が取りうる値の分布です。観測された差がその裾にあれば $H_0$ を棄却します。
信頼区間から検定を行うこともできますが、検定では必ず「$H_0$ が成り立つ世界」を作ってからリサンプリングするという点を外さないでください。単に元データをリサンプリングして得た分布は $H_0$ の分布ではありません。
Permutation test との違い
よく似た手法に並べ替え検定（permutation test）があります。
- ブートストラップ検定: 復元抽出（with replacement）。同じデータ点が複数回選ばれうる
- Permutation test: 非復元（without replacement）でラベルを並べ替える。各データ点はちょうど 1 回ずつ使われる
Permutation test は「2 群が交換可能」という仮定のもとで 有限標本でも厳密に正確な p 値を与えます。一方ブートストラップ検定は大標本で漸近的に正確です。したがって、2 群の分布が完全に同一かどうかを検定したいなら permutation test が優れています。しかし「平均だけが等しく分散は違うかもしれない」といった緩い帰無仮説や、区間推定と検定を同じ枠組みで扱いたい場合はブートストラップに分があります。
ここまでブートストラップの万能ぶりを見てきましたが、当然ながら限界があります。むしろここからが実務で最も大事な部分です。
ブートストラップが失敗するとき
一致性の条件
ブートストラップの一致性、すなわち
$$ \sup_t \left|P^*\left(\frac{\sqrt{n}(\bar{X}^* – \bar{X}_n)}{\hat{\sigma}} \leq t\right) – P\left(\frac{\sqrt{n}(\bar{X}_n – \mu)}{\sigma} \leq t\right)\right| \xrightarrow{p} 0 $$
が成り立つ（$P^*$ はブートストラップ確率）ためには、おおむね次の条件が必要です。
- 統計量が経験分布の摂動に対して滑らかであること（Hadamard 微分可能性）
- 母集団に有限な分散が存在すること
- 統計量の極限分布が正規分布であること
裏を返せば、これらが破れる場面ではブートストラップは動作するのに答えが間違う、という最も厄介な失敗をします。エラーが出ないので気づきにくいのです。
失敗例 1: 極値統計
代表的な破綻例が最大値 $X_{(n)} = \max(X_1, \dots, X_n)$ です。
一様分布 $U(0, \theta)$（$\theta = 10$）から $n = 50$ を取り、最大値をブートストラップした結果です。左のパネルで、赤い塗りつぶしの「真の標本分布」は $\theta = 10$ の直下に連続的に広がっているのに対し、オレンジのブートストラップ分布は元データの値の上にしか立っていません。しかも 64.1%（理論値 63.6%）が標本最大値 9.756 ちょうどに集中しています。これは先ほどの 63.2% の裏返しで、「標本最大値が再標本に少なくとも 1 回は入る確率」がそのまま $1-(1-1/n)^n$ になるためです。
結果は中央のパネルのとおり壊滅的で、パーセンタイル区間の被覆確率は $n = 20, 50, 100, 200$ のすべてで 0.000 でした。$n$ を 10 倍にしても改善しません。理由は明白で、ブートストラップ分布のどの値も標本最大値以下であり、標本最大値は必ず真の $\theta$ より小さいからです。区間の上限は構造的に $\theta$ に届きません（実測でも $n=200$ の平均上限は 9.949 で、$\theta = 10$ に 0.051 届いていない）。
最大値の極限分布は正規分布ではなく極値分布（一様分布の場合は指数型）なので、上の一致性条件に真正面から違反しています。極値統計を扱うなら $m$ out of $n$ ブートストラップ（$m \ll n$ の部分標本を使う）やサブサンプリングが必要です。
失敗例 2: 極端に小さい標本
$\hat{F}_n$ が $F$ の代役として機能するのは $n$ がそれなりに大きいときだけです。右のパネルは一様分布の母平均に対するパーセンタイル区間の被覆確率で、$n = 3$ で 0.739、$n = 5$ で 0.846、$n = 8$ で 0.895 と、名目 95% を大きく下回ります。$n = 15$ で 0.929、$n = 30$ で 0.937 とようやく近づきますが、それでも 95% には届いていません。
極端な話 $n = 2$ なら、ブートストラップ標本は $\{x_1,x_1\}, \{x_1,x_2\}, \{x_2,x_1\}, \{x_2,x_2\}$ の 4 通りしかなく、そこから 95% 区間を切り出すことに意味はありません。$n$ が 1 桁のときにブートストラップを使うのは避けるべきです。
そのほかの注意すべき場面
- 重い裾の分布: 母分散が存在しない（コーシー分布など）場合、一致性が破れます
- パラメータが定義域の境界にある場合: 「分散成分がゼロか」の検定など
- 独立でないデータ: これは頻度が高いので、次節で独立に扱います
そして、実務で最も遭遇しやすい失敗が「独立でないデータ」です。
時系列データ — ブロックブートストラップ
ブートストラップは「データが独立同分布である」ことを暗黙に仮定しています。復元抽出は観測の順序を完全にシャッフルするので、時間的な依存関係は破壊されます。
$\phi = 0.8$ の AR(1) 過程（$n = 200$）で標本平均の標準誤差を推定した結果です。左と中央のパネルのとおり、この系列は隣接時点が強く相関しており、自己相関は理論値 $\phi^k$ に沿ってゆっくり減衰します。ところが通常の（iid）ブートストラップが返す標準誤差は 0.1157 で、真の 0.3481 のわずか 33.2% でした。3 倍近く過小評価しています。
この 1/3 という数字は偶然ではありません。iid ブートストラップは周辺分散 $\sigma^2/(1-\phi^2)$ を $n$ で割った値を返すのに対し、真の分散は長期分散 $\sigma^2/(1-\phi)^2$ を $n$ で割った値です。比を取ると
$$ \frac{\text{SE}_{\text{iid}}}{\text{SE}_{\text{true}}} = \sqrt{\frac{1/(1-\phi^2)}{1/(1-\phi)^2}} = \sqrt{\frac{(1-\phi)^2}{(1-\phi)(1+\phi)}} = \sqrt{\frac{1-\phi}{1+\phi}} $$
となり、$\phi = 0.8$ を代入すると $\sqrt{0.2/1.8} = 1/3$。実測の 33.2% とぴたり一致します。「有効サンプルサイズが実際には 200 ではなく約 22 しかない」ということを、iid ブートストラップは見落としているのです。
対処法が ブロックブートストラップです。1 点ずつではなく、長さ $L$ の連続した区間（ブロック）を単位としてリサンプリングし、ブロック内の依存構造を保存します。右のパネルで、ブロック長を $L=1$（通常のブートストラップ）から伸ばしていくと、SE は 0.1157 → 0.1554（$L=2$）→ 0.2188（$L=5$）→ 0.2647（$L=10$）→ 0.2892（$L=20$、真の 83.1%） と改善します。ただし $L=40$ で 0.2799、$L=80$ で 0.2248 と、長くしすぎると再び悪化します。ブロックが長いと独立に取れるブロックの本数が減り、リサンプリングの多様性が失われるためです。
つまりブロック長には最適値があり、理論的には $L \propto n^{1/3}$ 程度が推奨されます。今回の $n = 200$ なら $L \approx 6$〜20 が目安で、実測のピーク位置と整合しています。なお $L = 20$ でも真の 83% にとどまり 100% には届きませんが、これはブロック境界で依存が切れることによるブロックブートストラップ固有の下向きバイアスです。「時系列にそのままブートストラップをかけると標準誤差が数分の 1 になる」 という事実だけは、絶対に覚えておいてください。
同じ議論は空間データやクラスタ構造を持つデータ（同じ学校の生徒、同じ患者からの複数測定）にも当てはまり、そこではクラスタブートストラップ（個体ではなくクラスタを丸ごとリサンプリング）が使われます。
最後に、統計学の外側で最も広く使われているブートストラップの応用を見ておきます。
機械学習との接点 — バギングと OOB 推定
バギング（bagging = bootstrap aggregating） は、ブートストラップ標本ごとにモデルを学習し、その予測を平均（分類なら多数決）する手法です。名前のとおりブートストラップそのもので、ランダムフォレストの基礎になっています。
面白いのは、ここで 63.2% / 36.8% の分割が思わぬ副産物を生むことです。各ブートストラップ標本には元データの約 63.2% しか含まれないので、残りの約 36.8% はその木にとって「一度も見ていないデータ」になります。これを OOB（out-of-bag） と呼び、各データ点をそれが OOB になっている木だけで予測すれば、検証データを切り出さずに汎化性能を推定できます。
右のパネルが実測で、訓練データ $n = 720$ に対して 1 本の木が使うデータの割合は 0.6318、OOB は 0.3682 と、理論値 $1 – 1/e = 0.6321$ と $1/e = 0.3679$ に小数第 3 位まで一致しました。左のパネルは OOB 誤り率とテスト誤り率の推移で、木が 5 本のときは OOB 0.2736 に対しテスト 0.1750 と大きく食い違いますが（各データ点を評価できる木が少なすぎるため）、木を増やすと 600 本で OOB 0.1306 / テスト 0.1208 まで接近します。OOB がわずかに悲観的なのは、各点を全体の約 36.8% の木だけで予測しているためで、これは既知の性質です。
灰色の点線が示すとおり、単一の決定木の平均誤り率は 0.2291 で、バギングによって 0.121 まで半減しています。同じデータを少しずつ違う形で使い回すだけで汎化性能が上がるのは、木の予測誤差が互いに独立に近く、平均によって分散が打ち消されるためです。
実装の骨格
以上をまとめると、実装は驚くほど短くなります。統計量を関数として受け取る形にしておけば、あとは何にでも使えます。
import numpy as np
from scipy import stats
def bootstrap(data, statistic, B=10000, alpha=0.05, seed=42):
"""統計量 statistic の SE・バイアス・各種95%信頼区間を返す。"""
rng = np.random.default_rng(seed)
data = np.asarray(data)
n = len(data)
theta_hat = statistic(data)
# --- ブートストラップ分布 ---
boot = np.array([statistic(data[rng.integers(0, n, n)]) for _ in range(B)])
se = boot.std(ddof=1)
bias = boot.mean() - theta_hat
# --- 区間3種 ---
ci_pct = np.percentile(boot, [100 * alpha / 2, 100 * (1 - alpha / 2)])
ci_basic = np.array([2 * theta_hat - ci_pct[1], 2 * theta_hat - ci_pct[0]])
# BCa: バイアス補正 z0 と加速度 a
p = np.clip(np.mean(boot < theta_hat), 1 / (2 * B), 1 - 1 / (2 * B))
z0 = stats.norm.ppf(p)
jack = np.array([statistic(np.delete(data, i)) for i in range(n)])
d = jack.mean() - jack
a = np.sum(d ** 3) / (6 * np.sum(d ** 2) ** 1.5)
zl, zu = stats.norm.ppf(alpha / 2), stats.norm.ppf(1 - alpha / 2)
a1 = stats.norm.cdf(z0 + (z0 + zl) / (1 - a * (z0 + zl)))
a2 = stats.norm.cdf(z0 + (z0 + zu) / (1 - a * (z0 + zu)))
ci_bca = np.percentile(boot, [100 * a1, 100 * a2])
return dict(estimate=theta_hat, se=se, bias=bias, z0=z0, a=a,
ci_percentile=ci_pct, ci_basic=ci_basic, ci_bca=ci_bca)
statistic に np.median を渡せば中央値、lambda x: np.percentile(x, 90) を渡せば 90 パーセンタイル、lambda x: x.std(ddof=1) / x.mean() を渡せば変動係数の推論ができます。関数を差し替えるだけで対象が変わるというのが、ブートストラップの実装上の最大の美点です。
なお SciPy には scipy.stats.bootstrap が用意されており、実務ではこちらを使えば十分です（method='BCa' が指定できます）。上のコードは中身を理解するためのものと考えてください。
まとめ
本記事では、ブートストラップ法の考え方から限界までを実測つきで見てきました。
- 基本原理: 未知の母集団 $F$ を経験分布 $\hat{F}_n$ で置き換え、復元抽出で統計量の標本分布を近似する（プラグイン原理）。手続きに統計量の種類が入らないので、何にでも同じコードが使える
- 63.2%: 1 回の再標本に入る元データの割合は $1-(1-1/n)^n \to 1-1/e = 0.6321$。実測でも $n=1000$ で 0.6322 と一致。この数字は OOB 推定として機械学習で再登場する
- 最大の価値: 中央値（SE 0.4190）、90% 分位点（0.6813）、変動係数（0.1075）、四分位範囲（0.4353）など、閉じた公式が存在しない統計量に標準誤差を与えられること。標本平均では教科書公式 0.3281 に対し 0.3266 と一致し、既存の理論と矛盾しない
- 信頼区間: 歪んだ母集団での被覆確率の実測では、$n = 25$ で BCa 0.899 > パーセンタイル 0.890 > 基本法 0.865。BCa がすべての $n$ で最良、基本法が最悪。ただし小標本では全手法が名目 95% を下回る（$n=10$ で 0.82〜0.85）
- 反復回数 $B$: 増やして減るのはモンテカルロ誤差だけ（$1/\sqrt{B}$）。SE の平均値は $B=50$ の 0.4261 から $B=10000$ の 0.4204 とほとんど動かない。SE なら $B=200$、区間なら $B=1000$〜$10000$
- バイアス補正: バイアス推定 $-0.1528$ は理論値 $-0.1600$ と一致し、補正で平均が 3.826 → 3.979 と真値 4.0 に寄る。ただし RMSE は 1.107 → 1.137 と悪化するので、補正は常に得とは限らない
- ジャックナイフ: 平均のような滑らかな統計量では同等だが、中央値では CV が $n=15$ の 0.65 から $n=125$ の 0.67 と縮まず、一致性を持たない
- 失敗する場面: 最大値では被覆確率が $n$ によらず 0.000。$n=3$ の平均でも 0.739。時系列では標準誤差が真の 33.2% と 3 倍の過小評価になり、ブロックブートストラップ（$L=20$ で 83.1% まで回復）が必要
ブートストラップは「分布の仮定から自由になる」魔法ではなく、「$\hat{F}_n \approx F$ という 1 つの仮定に賭ける」手法です。その仮定が成り立つ限り驚くほど広く使えますが、極値・極小標本・依存データという 3 つの場面では静かに間違った答えを返します。この線引きを持っておくことが、ブートストラップを使いこなす条件です。
次のステップとして、以下の記事も参考にしてください。
