ブートストラップ：コンピュータ時代の統計革命
ブートストラップリサンプリングは、観測データから繰り返しリサンプリングすることで、任意の統計量の標本分布を推定する強力な統計手法です。1979年にブラッドレイ・エフロンが提唱し、数学的な公式や分布の仮定に頼らずに複雑な統計量の分析を可能にすることで、統計的推測に革命をもたらしました。
ブートストラップの背後にある洞察は美しくシンプルです。標本は母集団の最良の推定値です。標本から（復元抽出で）リサンプリングすることで、母集団から繰り返しサンプリングした場合に何が起こるかをシミュレートします。このアプローチは標準偏差にとって特に価値があります。従来の信頼区間の公式は正規性を仮定しますが、実際にはこの仮定が成り立たないことが多いからです。
ブートストラップは現代のデータサイエンスで不可欠になっています。任意の統計量（中央値、相関、回帰係数、ニューラルネットワークの重みなど）に対して機能し、データの基礎となる分布について仮定を置きません。
なぜ標準偏差にブートストラップを使うのか？
標準偏差の従来の信頼区間は、データが正規分布に従うことを仮定しています。この仮定が破られると（よくあることです）、これらの区間は大幅に不正確になる可能性があります。ブートストラップは分布に依存しない代替手段を提供します。
従来の方法が失敗する場合
標準偏差にブートストラップを使う主な利点：
- 分布の仮定不要： 正規分布、歪んだ分布、裾の重い分布のいずれでも同等に機能
- 小標本での性能： n < 30ではパラメトリック法よりも精度が高いことが多い
- 複雑な統計量にも対応： トリム標準偏差、MAD、カスタム変動指標にも同じアプローチが使える
- 視覚的な洞察： ブートストラップ分布が数字だけでなく、何が起こっているかを見せてくれる
ブートストラップの手順
ブートストラップのアルゴリズムは驚くほど簡単です。n個の観測値からなる元の標本から：
ブートストラップ標本を抽出する
統計量を計算する
何千回も繰り返す
分布を分析する
なぜ復元抽出なのか？
何回のブートストラップ標本が必要か？ B = 1,000は大まかな推定や仮説検定に十分な場合が多いです。信頼区間にはB = 10,000で安定したパーセンタイルが得られます。論文品質のBCa区間にはB = 15,000以上が推奨されます。
ブートストラップ信頼区間の方法
ブートストラップ標本から信頼区間を構築する方法がいくつかあり、それぞれにトレードオフがあります。
1. パーセンタイル法（最も簡単）
最も直感的なアプローチです。ブートストラップ分布のパーセンタイルをそのまま使用します。
パーセンタイル信頼区間
10,000個のブートストラップ標本の場合、250番目と9,750番目の順序値です。シンプルですが、ブートストラップ分布が歪んでいる場合はバイアスが生じることがあります。
2. 基本（ピボタル）ブートストラップ
標本統計量とブートストラップ統計量の関係を利用します。
基本ブートストラップ信頼区間
θ̂は元の標本SDです。パーセンタイル区間を標本推定値の周りに「反転」させたものです。
3. BCa法（バイアス補正・加速法）
精度における最高水準。BCaはブートストラップ分布のバイアスと加速（パラメータ値に応じて標準誤差がどう変化するか）の両方を補正します。計算はより複雑ですが、2次の精度を持つ区間を提供します。
| 方法 | 長所 | 短所 |
|---|---|---|
| パーセンタイル法 | シンプルで直感的 | 歪んだデータではバイアスが生じうる |
| 基本法 | 対称な区間 | 負の値が出る可能性がある |
| BCa法 | 最も正確、変換に対して整合的 | 計算コストが高い |
計算例：非正規データ
応答時間（ミリ秒）の15個の測定値を考えます：245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312。このデータは右に歪んでいます（非常に遅い応答がいくつかある）。
標本SDを計算する
ブートストラップ標本を生成する
ブートストラップSDを計算する
パーセンタイルを求める
95%信頼区間を構成する
ブートストラップ信頼区間は非対称（高い側がより広い）で、データの右歪みの性質を反映しています。カイ2乗信頼区間はこの非対称性を捉えません。
Python実装
複数の信頼区間法による完全なブートストラップ実装：
import numpy as np
from scipy import stats
def bootstrap_sd_ci(data, n_bootstrap=10000, ci=0.95, method='percentile'):
"""
Bootstrap confidence interval for standard deviation.
Parameters:
-----------
data : array-like - Original sample
n_bootstrap : int - Number of bootstrap samples
ci : float - Confidence level (e.g., 0.95)
method : str - 'percentile', 'basic', or 'bca'
Returns:
--------
tuple : (lower_bound, upper_bound, bootstrap_sds)
"""
data = np.array(data)
n = len(data)
original_sd = np.std(data, ddof=1)
# Generate bootstrap samples and calculate SDs
bootstrap_sds = np.array([
np.std(np.random.choice(data, size=n, replace=True), ddof=1)
for _ in range(n_bootstrap)
])
alpha = 1 - ci
if method == 'percentile':
lower = np.percentile(bootstrap_sds, 100 * alpha/2)
upper = np.percentile(bootstrap_sds, 100 * (1 - alpha/2))
elif method == 'basic':
lower = 2*original_sd - np.percentile(bootstrap_sds, 100*(1-alpha/2))
upper = 2*original_sd - np.percentile(bootstrap_sds, 100*alpha/2)
elif method == 'bca':
# Bias correction
prop_less = np.mean(bootstrap_sds < original_sd)
z0 = stats.norm.ppf(prop_less)
# Acceleration (jackknife estimate)
jackknife_sds = np.array([
np.std(np.delete(data, i), ddof=1) for i in range(n)
])
jack_mean = jackknife_sds.mean()
a = np.sum((jack_mean - jackknife_sds)**3) / \
(6 * np.sum((jack_mean - jackknife_sds)**2)**1.5)
# Adjusted percentiles
z_alpha = stats.norm.ppf([alpha/2, 1-alpha/2])
adj_percentiles = stats.norm.cdf(
z0 + (z0 + z_alpha) / (1 - a*(z0 + z_alpha))
) * 100
lower = np.percentile(bootstrap_sds, adj_percentiles[0])
upper = np.percentile(bootstrap_sds, adj_percentiles[1])
return lower, upper, bootstrap_sds
# Example usage
response_times = [245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312]
for method in ['percentile', 'basic', 'bca']:
lower, upper, _ = bootstrap_sd_ci(response_times, method=method)
print(f"{method.upper():12s} 95% CI: [{lower:.1f}, {upper:.1f}]")