統計学で信頼区間を構成したり検定を行ったりする際、多くの場合は統計量の分布に関する理論的な結果（正規近似など）を利用します。しかし、分布の仮定が満たされない場合や、複雑な統計量の分布が理論的に求められない場合にはどうすればよいでしょうか。
この問いに対する強力な回答が ブートストラップ法（bootstrap method） です。1979年にBradley Efronが提案したこの手法は、手元のデータから復元抽出（リサンプリング）を繰り返すことで、統計量の分布を近似的に推定します。分布の仮定を必要とせず、ほぼすべての統計量に適用できる汎用性の高い方法です。
本記事の内容
- ブートストラップの基本原理（プラグイン原理）
- ブートストラップ標準誤差の理論
- パーセンタイル法・BCa法による信頼区間
- ブートストラップ検定（permutation test との違い）
- 大標本での一致性
- Pythonでのスクラッチ実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ブートストラップ法とは
統計的推論の本質
統計的推論の核心は、統計量の標本分布 を知ることにあります。例えば、標本平均 $\bar{X}$ の標準誤差を求めたいとき、もし母集団の分布がわかっていれば、$\bar{X}$ の分布を理論的に導出できます。
しかし、複雑な統計量（中央値、四分位範囲、相関係数の差など）の標本分布を理論的に導出するのは困難なことが多いです。
ブートストラップの基本アイデア
ブートストラップの核心的なアイデアは次のとおりです。
「もし母集団から何度も標本を取り直せたら、統計量の分布がわかるのに」
実際には母集団から繰り返し標本を取ることはできません。しかし、手元のデータ自体を母集団の代理として扱い、そこからの復元抽出で擬似的な再標本を生成する ことで、統計量の分布を近似できます。
アルゴリズム
元の標本を $\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}$ の分布がブートストラップ分布
プラグイン原理
経験分布関数
$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) $$
プラグイン原理
母集団の特性が分布関数 $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$
ブートストラップとプラグイン原理
統計量 $\hat{\theta} = s(\bm{X})$ の標本分布は $F$ に依存します。真の標本分布を $G_F$ と書くと、
$$ G_F(\cdot) = P_F(s(\bm{X}) \leq \cdot) $$
ブートストラップは $F$ を $\hat{F}_n$ に置き換えて、
$$ \hat{G}_n(\cdot) = P_{\hat{F}_n}(s(\bm{X}^*) \leq \cdot) $$
を計算します。$\bm{X}^*$ は $\hat{F}_n$ からの復元抽出（= 元のデータからの復元抽出）です。
$\hat{G}_n$ を解析的に求めるのは一般に困難ですが、モンテカルロシミュレーション（$B$ 回のリサンプリング）で近似できます。
ブートストラップ標準誤差
理論
統計量 $\hat{\theta}$ の標準誤差は
$$ \text{SE}_F(\hat{\theta}) = \sqrt{\text{Var}_F[s(\bm{X})]} $$
ブートストラップ標準誤差は、これを $\hat{F}_n$ で置き換えたものです。
$$ \text{SE}_{\text{boot}}(\hat{\theta}) = \sqrt{\text{Var}_{\hat{F}_n}[s(\bm{X}^*)]} $$
モンテカルロ近似として、
$$ \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$ の分散は
$$ \text{Var}(\hat{F}_n) = \frac{1}{n}\sum_{i=1}^{n}(x_i – \bar{x})^2 = \hat{\sigma}^2 $$
ブートストラップ標準誤差は
$$ \text{SE}_{\text{boot}} = \hat{\sigma}/\sqrt{n} $$
真の $\sigma/\sqrt{n}$ を $\hat{\sigma}/\sqrt{n}$ で推定しており、これは標準的なプラグイン推定に一致します。
ブートストラップ信頼区間
パーセンタイル法（Percentile Method）
最も単純なブートストラップ信頼区間です。ブートストラップ分布の $\alpha/2$ パーセンタイルと $1 – \alpha/2$ パーセンタイルを端点とします。
$$ CI_{\text{percentile}} = \left[\hat{\theta}^*_{(\alpha/2)}, \hat{\theta}^*_{(1-\alpha/2)}\right] $$
ここで $\hat{\theta}^*_{(q)}$ はブートストラップ推定値の $q$ 分位点です。
95%信頼区間であれば、$B$ 個のブートストラップ推定値を小さい順に並べ、2.5パーセンタイルと97.5パーセンタイルを取ります。
パーセンタイル法の問題点
パーセンタイル法は、$\hat{\theta}$ がバイアスを持つ場合や、$\hat{\theta}$ の分散が $\hat{\theta}$ の値に依存する場合に、被覆確率（coverage probability）が名目水準から乖離します。
BCa法（Bias-Corrected and Accelerated Method）
BCa法は、パーセンタイル法のバイアスと歪みを補正した方法です。
$$ 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) $$
$$ \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$:
$$ \hat{z}_0 = \Phi^{-1}\left(\frac{\#\{\hat{\theta}^{*b} < \hat{\theta}\}}{B}\right) $$
これは、ブートストラップ推定値のうち元の推定値 $\hat{\theta}$ より小さいものの割合を、正規分布の分位点に変換したものです。$\hat{\theta}$ にバイアスがなければ、ちょうど半分の $\hat{\theta}^{*b}$ が $\hat{\theta}$ より小さいので $\hat{z}_0 = 0$ になります。
加速度定数 $\hat{a}$（acceleration constant）:
$$ \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$ 番目のデータを除いた ジャックナイフ推定値 で、
$$ \hat{\theta}_{(-i)} = s(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n) $$
$\bar{\theta}_{(\cdot)} = \frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{(-i)}$ はジャックナイフ推定値の平均です。
加速度定数 $\hat{a}$ は、統計量の分布の 歪度（skewness） を捕捉します。歪みがなければ $\hat{a} = 0$ です。
$\hat{z}_0 = 0$ かつ $\hat{a} = 0$ の場合、BCa法はパーセンタイル法に一致します。
ブートストラップ検定
基本的な枠組み
帰無仮説 $H_0: \theta = \theta_0$ を検定するとします。観測された検定統計量を $t_{\text{obs}}$ とすると、ブートストラップ p 値は次のように計算します。
- 帰無仮説のもとでのブートストラップ: データを帰無仮説が成り立つように変換してからリサンプリング
- 各ブートストラップ標本から検定統計量 $t^{*b}$ を計算
- p値 = $\frac{\#\{t^{*b} \geq t_{\text{obs}}\}}{B}$
2標本問題での実装
2群の平均が等しいかを検定する場合、帰無仮説のもとでは2群のデータは同じ分布から来ています。
- 2群のデータをプールする
- プールされたデータから復元抽出で2群に分け直す
- 2群の平均差を計算
これを $B$ 回繰り返し、観測された平均差以上に大きい差が何回得られたかの割合がp値です。
Permutation test との違い
ブートストラップ検定: 復元抽出（with replacement）を行う。同じデータ点が複数回選ばれうる
Permutation test: 非復元（without replacement）でデータのラベルを並べ替える。各データ点はちょうど1回ずつ使われる
Permutation test は帰無仮説のもとでの正確な分布を生成する（有限標本で正確）のに対し、ブートストラップ検定は大標本で漸近的に正確です。
大標本での一致性
ブートストラップの一致性
正則条件のもとで、ブートストラップ分布は真の標本分布に漸近的に一致します。
具体的には、標本平均 $\hat{\theta}_n = \bar{X}_n$ に対して、
$$ \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^*$ はブートストラップ確率（$\hat{F}_n$ のもとでの確率）です。
この結果は、ブートストラップ分布が真の標本分布を一様に近似することを示しています。
一致性が成り立つ条件
ブートストラップの一致性が成り立つには、以下が必要です。
- 統計量が十分に「滑らか」であること（経験分布の小さな変化に対して連続的に変化する）
- 母集団に有限な分散が存在すること
一致性が破れる例
ブートストラップが一致しない例も存在します。例えば、最大値 $X_{(n)} = \max(X_1, \dots, X_n)$ の分布のブートストラップは、一般に一致しません。$X_{(n)}$ は経験分布の端に強く依存するため、プラグイン原理が適切に機能しないためです。
具体例
数値例: 中央値の信頼区間
20個のデータ $\{3.2, 5.1, 2.8, 7.3, 4.5, 6.1, 3.9, 8.7, 5.6, 4.2, 9.1, 2.5, 6.8, 3.7, 7.9, 5.3, 4.8, 6.5, 3.1, 8.2\}$ に対する中央値の95%信頼区間を、ブートストラップ法で求めます。
標本中央値: $\hat{\theta} = 5.2$
$B = 10000$ 回のリサンプリングを行い、各ブートストラップ標本の中央値を計算します。ブートストラップ中央値の2.5パーセンタイルと97.5パーセンタイルが信頼区間の端点です。
中央値のような統計量には理論的な標準誤差の公式がないため、ブートストラップは特に有用です。
Pythonでの実装
ブートストラップの基本実装
import numpy as np
from scipy import stats
def bootstrap(data, statistic, B=10000, confidence=0.95, random_state=42):
"""
ブートストラップ法の基本実装
Parameters
----------
data : array-like
元のデータ
statistic : callable
統計量を計算する関数（配列を受け取りスカラーを返す）
B : int
ブートストラップの反復回数
confidence : float
信頼水準
random_state : int
乱数のシード
Returns
-------
dict : ブートストラップの結果
"""
rng = np.random.RandomState(random_state)
data = np.array(data)
n = len(data)
# 元のデータでの統計量
theta_hat = statistic(data)
# ブートストラップ
boot_stats = np.zeros(B)
for b in range(B):
# 復元抽出
boot_sample = data[rng.randint(0, n, size=n)]
boot_stats[b] = statistic(boot_sample)
# 標準誤差
se_boot = np.std(boot_stats, ddof=1)
# パーセンタイル法による信頼区間
alpha = 1 - confidence
ci_percentile = np.percentile(boot_stats, [100*alpha/2, 100*(1-alpha/2)])
# 正規近似法による信頼区間
z = stats.norm.ppf(1 - alpha/2)
ci_normal = [theta_hat - z * se_boot, theta_hat + z * se_boot]
# バイアスの推定
bias = np.mean(boot_stats) - theta_hat
return {
'estimate': theta_hat,
'se': se_boot,
'bias': bias,
'ci_percentile': ci_percentile,
'ci_normal': ci_normal,
'boot_stats': boot_stats
}
# --- 具体例: 中央値の信頼区間 ---
data = np.array([3.2, 5.1, 2.8, 7.3, 4.5, 6.1, 3.9, 8.7, 5.6, 4.2,
9.1, 2.5, 6.8, 3.7, 7.9, 5.3, 4.8, 6.5, 3.1, 8.2])
result = bootstrap(data, np.median, B=10000)
print("=== ブートストラップ法: 中央値の推定 ===")
print(f"標本中央値: {result['estimate']:.4f}")
print(f"ブートストラップ標準誤差: {result['se']:.4f}")
print(f"バイアス推定値: {result['bias']:.4f}")
print(f"95%信頼区間（パーセンタイル法）: [{result['ci_percentile'][0]:.4f}, {result['ci_percentile'][1]:.4f}]")
print(f"95%信頼区間（正規近似法）: [{result['ci_normal'][0]:.4f}, {result['ci_normal'][1]:.4f}]")
# 標本平均との比較
result_mean = bootstrap(data, np.mean, B=10000)
theoretical_se = np.std(data, ddof=1) / np.sqrt(len(data))
print(f"\n=== 比較: 標本平均の推定 ===")
print(f"標本平均: {result_mean['estimate']:.4f}")
print(f"ブートストラップSE: {result_mean['se']:.4f}")
print(f"理論的SE (s/√n): {theoretical_se:.4f}")
BCa法の実装
import numpy as np
from scipy import stats
def bootstrap_bca(data, statistic, B=10000, confidence=0.95, random_state=42):
"""
BCa法（Bias-Corrected and Accelerated）ブートストラップ信頼区間
Parameters
----------
data : array-like
元のデータ
statistic : callable
統計量を計算する関数
B : int
ブートストラップの反復回数
confidence : float
信頼水準
random_state : int
乱数のシード
Returns
-------
ci_bca : tuple
BCa信頼区間の下限と上限
"""
rng = np.random.RandomState(random_state)
data = np.array(data)
n = len(data)
# 元のデータでの統計量
theta_hat = statistic(data)
# ブートストラップ
boot_stats = np.zeros(B)
for b in range(B):
boot_sample = data[rng.randint(0, n, size=n)]
boot_stats[b] = statistic(boot_sample)
# --- バイアス補正係数 z0 ---
# θ*b < θ_hat となる割合
prop_less = np.mean(boot_stats < theta_hat)
# 0や1を避ける
prop_less = np.clip(prop_less, 1/(2*B), 1 - 1/(2*B))
z0 = stats.norm.ppf(prop_less)
# --- 加速度定数 a （ジャックナイフ法）---
jackknife_stats = np.zeros(n)
for i in range(n):
# i番目を除いたデータ
jack_data = np.delete(data, i)
jackknife_stats[i] = statistic(jack_data)
jack_mean = np.mean(jackknife_stats)
diff = jack_mean - jackknife_stats
# 歪度に基づく加速度
a_hat = np.sum(diff**3) / (6 * (np.sum(diff**2))**1.5)
# --- BCa信頼区間の計算 ---
alpha = 1 - confidence
z_alpha_lower = stats.norm.ppf(alpha / 2)
z_alpha_upper = stats.norm.ppf(1 - alpha / 2)
# 補正されたパーセンタイル
alpha1 = stats.norm.cdf(z0 + (z0 + z_alpha_lower) / (1 - a_hat * (z0 + z_alpha_lower)))
alpha2 = stats.norm.cdf(z0 + (z0 + z_alpha_upper) / (1 - a_hat * (z0 + z_alpha_upper)))
ci_bca = np.percentile(boot_stats, [100*alpha1, 100*alpha2])
return ci_bca, z0, a_hat, boot_stats
# --- 具体例: 歪んだ分布でのBCa法の効果 ---
np.random.seed(42)
# 対数正規分布からのデータ（右に歪んだ分布）
data_skewed = np.random.lognormal(mean=1.0, sigma=0.8, size=30)
print(f"データの歪度: {stats.skew(data_skewed):.4f}")
# 平均のブートストラップ信頼区間
theta_hat = np.mean(data_skewed)
ci_bca, z0, a_hat, boot_stats = bootstrap_bca(data_skewed, np.mean, B=10000)
ci_pct = np.percentile(boot_stats, [2.5, 97.5])
se_boot = np.std(boot_stats, ddof=1)
z_crit = stats.norm.ppf(0.975)
ci_normal = [theta_hat - z_crit * se_boot, theta_hat + z_crit * se_boot]
print(f"\n=== 歪んだ分布での信頼区間比較 ===")
print(f"標本平均: {theta_hat:.4f}")
print(f"バイアス補正係数 z₀ = {z0:.4f}")
print(f"加速度定数 â = {a_hat:.4f}")
print(f"95%信頼区間:")
print(f" 正規近似法: [{ci_normal[0]:.4f}, {ci_normal[1]:.4f}]")
print(f" パーセンタイル: [{ci_pct[0]:.4f}, {ci_pct[1]:.4f}]")
print(f" BCa法: [{ci_bca[0]:.4f}, {ci_bca[1]:.4f}]")
ブートストラップ分布の可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# データ生成
data = np.random.exponential(scale=3.0, size=25)
B = 10000
# 異なる統計量のブートストラップ
statistics = {
'Mean': np.mean,
'Median': np.median,
'Std Dev': lambda x: np.std(x, ddof=1),
'IQR': lambda x: np.percentile(x, 75) - np.percentile(x, 25)
}
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()
for idx, (name, stat_func) in enumerate(statistics.items()):
# ブートストラップ
boot_vals = np.zeros(B)
for b in range(B):
boot_sample = data[np.random.randint(0, len(data), size=len(data))]
boot_vals[b] = stat_func(boot_sample)
theta_hat = stat_func(data)
se_boot = np.std(boot_vals, ddof=1)
ci = np.percentile(boot_vals, [2.5, 97.5])
# ヒストグラム
axes[idx].hist(boot_vals, bins=50, density=True, alpha=0.7,
color='steelblue', edgecolor='black')
axes[idx].axvline(theta_hat, color='red', linewidth=2.5,
label=f'Original: {theta_hat:.3f}')
axes[idx].axvline(ci[0], color='orange', linestyle='--', linewidth=2,
label=f'95% CI: [{ci[0]:.3f}, {ci[1]:.3f}]')
axes[idx].axvline(ci[1], color='orange', linestyle='--', linewidth=2)
axes[idx].set_title(f'{name} (SE = {se_boot:.3f})', fontsize=13)
axes[idx].set_xlabel(f'Bootstrap {name}', fontsize=11)
axes[idx].set_ylabel('Density', fontsize=11)
axes[idx].legend(fontsize=9)
axes[idx].grid(True, alpha=0.3)
plt.suptitle(f'Bootstrap Distributions (n = {len(data)}, B = {B})',
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
ブートストラップ検定の実装
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
def bootstrap_two_sample_test(x, y, B=10000, random_state=42):
"""
2標本ブートストラップ検定
H0: μ_x = μ_y
Parameters
----------
x, y : array-like
2群のデータ
B : int
ブートストラップの反復回数
Returns
-------
t_obs : float
観測された検定統計量（平均差）
p_value : float
ブートストラップ p値
boot_diffs : ndarray
ブートストラップ平均差の配列
"""
rng = np.random.RandomState(random_state)
x, y = np.array(x), np.array(y)
n1, n2 = len(x), len(y)
# 観測された平均差
t_obs = np.mean(x) - np.mean(y)
# H0のもとでリサンプリング: プールしたデータから復元抽出
pooled = np.concatenate([x, y])
boot_diffs = np.zeros(B)
for b in range(B):
# プールデータからの復元抽出
boot_pool = pooled[rng.randint(0, len(pooled), size=len(pooled))]
boot_x = boot_pool[:n1]
boot_y = boot_pool[n1:]
boot_diffs[b] = np.mean(boot_x) - np.mean(boot_y)
# 両側p値
p_value = np.mean(np.abs(boot_diffs) >= np.abs(t_obs))
return t_obs, p_value, boot_diffs
# --- 具体例 ---
group_a = np.random.normal(50, 10, 25)
group_b = np.random.normal(55, 10, 30)
t_obs, p_boot, boot_diffs = bootstrap_two_sample_test(group_a, group_b)
# t検定との比較
t_stat, p_t = stats.ttest_ind(group_a, group_b)
print("=== 2標本ブートストラップ検定 ===")
print(f"群A: n = {len(group_a)}, x̄ = {np.mean(group_a):.4f}")
print(f"群B: n = {len(group_b)}, x̄ = {np.mean(group_b):.4f}")
print(f"観測された平均差: {t_obs:.4f}")
print(f"ブートストラップ p値: {p_boot:.4f}")
print(f"t検定 p値: {p_t:.4f}")
# 可視化
fig, ax = plt.subplots(figsize=(10, 6))
ax.hist(boot_diffs, bins=60, density=True, alpha=0.7, color='steelblue',
edgecolor='black', label='Bootstrap null distribution')
ax.axvline(t_obs, color='red', linewidth=2.5,
label=f'Observed diff = {t_obs:.3f}')
ax.axvline(-t_obs, color='red', linewidth=2.5, linestyle='--')
# 棄却域を塗りつぶし
boot_sorted = np.sort(boot_diffs)
mask_right = boot_sorted >= np.abs(t_obs)
mask_left = boot_sorted <= -np.abs(t_obs)
ax.fill_between(boot_sorted[mask_right],
0, stats.gaussian_kde(boot_diffs)(boot_sorted[mask_right]),
color='red', alpha=0.3, label=f'p-value region (p = {p_boot:.4f})')
ax.fill_between(boot_sorted[mask_left],
0, stats.gaussian_kde(boot_diffs)(boot_sorted[mask_left]),
color='red', alpha=0.3)
ax.set_xlabel('Mean Difference', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Bootstrap Two-Sample Test (Null Distribution)', fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
被覆確率のシミュレーション
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n_simulations = 1000
B = 2000
sample_sizes = [10, 20, 50, 100]
confidence = 0.95
# 対数正規分布（歪んだ分布）のパラメータ
mu_ln = 1.0
sigma_ln = 0.8
true_mean = np.exp(mu_ln + sigma_ln**2 / 2) # 対数正規分布の平均
methods = ['Normal', 'Percentile']
colors = {'Normal': 'blue', 'Percentile': 'red'}
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()
for idx, n in enumerate(sample_sizes):
coverage = {m: 0 for m in methods}
for sim in range(n_simulations):
data = np.random.lognormal(mu_ln, sigma_ln, n)
theta_hat = np.mean(data)
# ブートストラップ
boot_means = np.zeros(B)
for b in range(B):
boot_sample = data[np.random.randint(0, n, size=n)]
boot_means[b] = np.mean(boot_sample)
se = np.std(boot_means, ddof=1)
# 正規近似法
z = stats.norm.ppf(0.975)
ci_normal = [theta_hat - z * se, theta_hat + z * se]
if ci_normal[0] <= true_mean <= ci_normal[1]:
coverage['Normal'] += 1
# パーセンタイル法
ci_pct = np.percentile(boot_means, [2.5, 97.5])
if ci_pct[0] <= true_mean <= ci_pct[1]:
coverage['Percentile'] += 1
# 被覆確率
for m in methods:
coverage[m] /= n_simulations
# 棒グラフ
x_pos = np.arange(len(methods))
bars = axes[idx].bar(x_pos, [coverage[m] for m in methods],
color=[colors[m] for m in methods], alpha=0.7,
edgecolor='black')
axes[idx].axhline(y=confidence, color='gray', linestyle='--', linewidth=1.5,
label=f'Nominal = {confidence}')
for bar, m in zip(bars, methods):
axes[idx].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005,
f'{coverage[m]:.3f}', ha='center', fontsize=11)
axes[idx].set_xticks(x_pos)
axes[idx].set_xticklabels(methods, fontsize=11)
axes[idx].set_ylabel('Coverage Probability', fontsize=11)
axes[idx].set_title(f'n = {n}', fontsize=13)
axes[idx].set_ylim(0.85, 1.0)
axes[idx].legend(fontsize=10)
axes[idx].grid(True, alpha=0.3, axis='y')
plt.suptitle('Bootstrap CI Coverage for Lognormal Mean',
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
ブートストラップ反復回数 $B$ の影響
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# データ
data = np.random.exponential(scale=5.0, size=30)
theta_hat = np.median(data)
B_values = [50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000]
n_trials = 50 # 各Bでの試行回数
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
se_means = []
se_stds = []
ci_width_means = []
for B in B_values:
ses = []
widths = []
for trial in range(n_trials):
boot_vals = np.zeros(B)
for b in range(B):
boot_sample = data[np.random.randint(0, len(data), size=len(data))]
boot_vals[b] = np.median(boot_sample)
ses.append(np.std(boot_vals, ddof=1))
ci = np.percentile(boot_vals, [2.5, 97.5])
widths.append(ci[1] - ci[0])
se_means.append(np.mean(ses))
se_stds.append(np.std(ses))
ci_width_means.append(np.mean(widths))
# SEのばらつき
ax1.errorbar(B_values, se_means, yerr=se_stds, fmt='o-', capsize=5,
linewidth=2, color='steelblue', markersize=6)
ax1.set_xlabel('Number of Bootstrap Replicates ($B$)', fontsize=12)
ax1.set_ylabel('Bootstrap SE (mean $\\pm$ std)', fontsize=12)
ax1.set_title('Stability of Bootstrap SE', fontsize=13)
ax1.set_xscale('log')
ax1.grid(True, alpha=0.3)
# 信頼区間の幅
ax2.plot(B_values, ci_width_means, 'ro-', linewidth=2, markersize=6)
ax2.set_xlabel('Number of Bootstrap Replicates ($B$)', fontsize=12)
ax2.set_ylabel('Mean CI Width', fontsize=12)
ax2.set_title('Stability of CI Width', fontsize=13)
ax2.set_xscale('log')
ax2.grid(True, alpha=0.3)
plt.suptitle('Effect of Bootstrap Replicates $B$ on Inference',
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
print("Bが大きいほどブートストラップの結果は安定します。")
print("一般に、SEの推定にはB=200程度、信頼区間にはB=1000〜10000程度が推奨されます。")
まとめ
本記事では、ブートストラップ法の理論と実装について解説しました。
- 基本原理: 手元のデータから復元抽出を繰り返し、統計量の標本分布を近似する（プラグイン原理）
- 標準誤差: ブートストラップ推定値のばらつきから標準誤差を推定。標本平均の場合は $\hat{\sigma}/\sqrt{n}$ と一致
- 信頼区間: パーセンタイル法（単純）、BCa法（バイアスと歪みを補正、推奨）
- ブートストラップ検定: 帰無仮説のもとでリサンプリングし、p値を計算。Permutation test とは復元/非復元の違い
- 一致性: 正則条件のもとで、ブートストラップ分布は真の標本分布に漸近的に一致
- 実用: $B = 1000$〜$10000$ 程度の反復が一般的。分布の仮定を置けない場合に特に有用
次のステップとして、以下の記事も参考にしてください。