ブートストラップ信頼区間｜シミュレーション

🎓 レベル：発展　|　重要度：A（必須）

📎 前提：ブートストラップ法　|　統計：区間推定（母平均・母比率・母分散の信頼区間）

要点（BLUF）

ブートストラップ分布の分位点から信頼区間を作れます。最も簡単なのがパーセンタイル法（ブート統計量の2.5%・97.5%点をそのまま使う）。
ただし歪んだ分布ではパーセンタイル法は公称95%より過小被覆になります（実測93.4%）。
バイアスと歪みを補正する BCa 法（bias-corrected and accelerated）がより正確。過小被覆の実証で、なぜ補正が要るかを示します。

1. パーセンタイル法

ブートストラップで得た統計量の分布 $\{\hat{\theta}^*_1, \dots, \hat{\theta}^*_B\}$ から、信頼水準 $1-\alpha$ の区間を分位点で直接作ります。

\text{CI}_{1-\alpha} = \left[\ \hat{\theta}^*_{(\alpha/2)},\quad \hat{\theta}^*_{(1-\alpha/2)}\ \right]

95% なら、ブート統計量を並べた2.5パーセンタイルと97.5パーセンタイルを両端にするだけ。正規近似（ $\hat{\theta} \pm 1.96\,\text{SE}$ 、推定量の信頼区間）と違い、分布の非対称性をある程度反映でき、変換に対して不変という利点もあります。

2. 被覆率の検証：歪んだ分布での過小被覆

信頼区間が正しければ、95%区間は真値を約95%カバーするはず（被覆率）。歪んだ母集団（指数分布、平均2）でパーセンタイル法を検証します。

import numpy as np

# 乱数シードを固定
rng = np.random.default_rng(54)

true_mean = 2.0
n = 100; B = 2000; trials = 2000
covered = 0
for t in range(trials):
    data = rng.exponential(true_mean, n)
    idx = rng.integers(0, n, size=(B, n))
    boot_means = data[idx].mean(axis=1)
    lo, hi = np.percentile(boot_means, [2.5, 97.5])     # パーセンタイル区間
    if lo <= true_mean <= hi:
        covered += 1

print(f"パーセンタイル法 95%CI 被覆率 = {covered/trials:.4f}  (目標 0.95)")

出力：

パーセンタイル法 95%CI 被覆率 = 0.9340  (目標 0.95)

出力の意味：被覆率 0.934 で、目標の0.95を下回っています（過小被覆）。原因は母集団の歪み——指数分布は右に裾を引くので、標本平均の分布も非対称になり、パーセンタイル法の単純な分位点では真値を取りこぼしやすい。標本サイズ $n=100$ でも残るこのズレが、より洗練された区間（BCa）が必要な理由です。正規母集団なら、パーセンタイル法でも0.95付近になります。

3. BCa 法（バイアス・加速補正）

BCa（Bias-Corrected and accelerated）は、パーセンタイル法に2つの補正を加えます。

バイアス補正 $\hat{z}_0$ ：ブート分布が中央からずれている度合い（点推定より小さいブート値の割合から計算）。
加速 $\hat{a}$ ：統計量の分散が母数とともに変わる度合い（歪み）。ジャックナイフ（ジャックナイフ法）で推定。

これらを使って、使うパーセンタイルの位置を調整します：

\alpha_1 = \Phi\!\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{\alpha/2})}\right),\quad \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{\theta}^*_{(\alpha_1)}, \hat{\theta}^*_{(\alpha_2)}]$ を区間とします。 $\hat{z}_0 = 0$ （バイアスなし）かつ $\hat{a}=0$ （歪みなし）ならパーセンタイル法に一致します。BCa は変換不変性を保ちつつ被覆率を改善し、実用上の標準です（scipy.stats.bootstrap の既定）。理論的背景は統計へ。

4. 手法の使い分け

方法	計算	歪みへの強さ	備考
正規近似	SEだけ	弱い	対称前提（推定量の信頼区間）
パーセンタイル	分位点	中	簡単・変換不変
BCa	+ジャックナイフ	強い	実用標準・推奨
ブートストラップ-t	スチューデント化	強い	区間が外れることも

小標本・歪んだ統計量では BCa かブートストラップ-t を。区間の例として、上のコードで作ったパーセンタイル区間は標本平均2.097に対し $[1.708, 2.486]$ のような形になります。

数式の直観的意味

パーセンタイル法は「ブート分布の形をそのまま信頼区間にする」素朴さゆえに、分布が歪むと裾の片側を取りこぼします。BCa の補正は「ブート分布の歪みとバイアスを見て、左右のパーセンタイルの切り取り位置をずらす」操作。バイアス補正 $\hat{z}_0$ は「ブート分布が点推定に対してどちらに偏っているか」、加速 $\hat{a}$ は「標準誤差が母数の場所によってどう変わるか（裾の重さ）」を捉える。指数分布のように右に裾を引く統計量では、区間を右に少し広げ左を詰めることで被覆率が回復する。過小被覆 0.934 という数字は、「ブートストラップは万能ではなく、区間の作り方で精度が変わる」ことの教訓——リサンプリングしたら終わりではなく、分布の歪みまで補正して初めて正しい推測になります。

⚠️ よくある誤解・落とし穴

「パーセンタイル法はいつでも95%」ではない：歪んだ分布・小標本では過小被覆（本ノートで0.934）。BCa を検討。
「ブートストラップなら区間も自動的に正しい」ではない：区間の構成法で被覆率が変わります。BCa が実用標準。
「 $B$ を増やせば被覆率が改善する」ではない： $B$ はモンテカルロ誤差を減らすだけ。被覆率の偏りは構成法と $n$ の問題で、 $B$ では治りません。
「BCa は常にパーセンタイルより良い」ではない：多くの場合は改善しますが、極端な小標本では加速の推定が不安定になることも。
「区間が点推定に対して対称」ではない：ブート区間は非対称になりえます（むしろそれが利点）。対称を期待しないこと。

対応シミュレーション参照

本文のパーセンタイル法の被覆率検証（default_rng(54)、歪んだ分布で0.934の過小被覆）。