ブートストラップ(bootstrap)とは、ブーツ(boot)のストラップ(strap)、つまり靴ひものことです。英語に “pull oneself up by one’s bootstraps” という成句があります。自分の靴ひもを引っ張って自分自身を持ち上げるという意味で、有名な『ほら吹き男爵の冒険』の主人公ミュンヒハウゼン男爵が自分の髪の毛を引っ張って沼から脱出したという逸話の別バージョンのようです。
ここでいうブートストラップは、1979年に米国の統計学者エフロン(Bradley Efron, 1938--)が考えた方法です。大きさ $n$ の一つのサンプルから、復元抽出で大きさ $n$ のサンプルをたくさん作り出し、それでサンプル誤差などを推測します。母集団が正規分布から大きく外れていても使える方法です。
正規分布から大きく外れる分布のデータを用意します。よく例として使われるのが対数正規分布です。これは $\log X \sim N(\mu, \sigma^2)$ に従う $X$ の分布で、Python では例えば
import numpy as np rng = np.random.default_rng() n = 100 x = rng.lognormal(mean=0, sigma=1, size=n)
で生成できます。平均値は $\exp(\mu + \sigma^2/2)$ になります。右裾が長く伸びた分布です。
一方で、教育統計などで使うデータは、テストの点数にせよLikertデータにせよ、裾は有限です。左右対称に近ければサンプル平均は $t$ 分布が良い近似になりますが、ここでは左右対称でない極端な例として、次のような ◺ 型の分布を考えます。
import numpy as np rng = np.random.default_rng(20251220) n = 100 x = np.abs(rng.random(n) - rng.random(n))
密度関数は $-2x+2$ ($0 \leq x \leq 1$) です。$x$ の平均は 1/3、分散は 1/18 になるはずです。
まずは $x$ の平均値の分布を $t$ 分布で近似し、95%信頼区間を求めてみましょう。
from scipy import stats low, high = stats.ttest_1samp(x, 0).confidence_interval() print(low, high)
0.268701576150326 0.36325204150481255
平均値の95%信頼区間は $[0.269, 0.363]$ です。
この x からブートストラップサンプルを一つ生成するには、rng.choice(x, n, replace=True) または同じことですが x[rng.integers(0, n, size=n)] とします。後者のほうがやや速いはずです。このサンプルから平均を求めることを B 回繰り返し、その結果を配列 boot に代入します。リストでよければ
B = 10000 boot = [np.mean(x[rng.integers(0, n, size=n)]) for _ in range(B)]
NumPy配列にするには右辺を np.array() で囲んでもいいのですが、より効率的な次のような書き方もできます。
B = 10000 boot = np.fromiter((np.mean(x[rng.integers(0, n, size=n)]) for _ in range(B)), float, B)
これで
low, high = np.quantile(boot, [0.025, 0.975]) print(low, high)
0.26856405025665725 0.3623282666882384
ベクトル化すればもっと速くなりますが、メモリは消費します:
indices = rng.integers(0, n, size=(B, n)) boot = x[indices].mean(axis=1) low, high = np.quantile(boot, [0.025, 0.975]) print(low, high)
大きなデータでは次のようなバッチ処理が安全です:
batch_size = 1000
boot = np.empty(B)
for start in range(0, B, batch_size):
stop = min(start + batch_size, B)
indices = rng.integers(0, n, size=(stop - start, n))
boot[start:stop] = x[indices].mean(axis=1)
low, high = np.quantile(boot, [0.025, 0.975])
print(low, high)
上の方法で求めた平均値の95%信頼区間は $[0.269, 0.362]$ です。この方法をパーセンタイル信頼区間(percentile confidence interval, percentile CI)といいます。
実は、これより少し良い次の方法があります。
low, high = np.quantile(boot, [0.025, 0.975]) m = np.mean(boot) print(2 * m - high, 2 * m - low)
0.2693755468507872 0.36313976328236836
平均値の95%信頼区間は $[0.269, 0.363]$ です。この方法をベーシック信頼区間(basic CI)または反転パーセンタイル信頼区間(reverse percentile CI)といいます。
さらに改良したものとして、BCa (Bias-Corrected and accelerated) と呼ばれる方法や、Studentized bootstrap (bootstrap-$t$) と呼ばれる方法があります。次はBCaの例です:
import numpy as np
from scipy.stats import bootstrap
# サンプルデータ
rng = np.random.default_rng(20251220)
n = 100
x = np.abs(rng.random(n) - rng.random(n))
result = bootstrap(
(x,),
np.mean,
confidence_level=0.95,
n_resamples=10000,
method="BCa",
vectorized=True,
batch=1000,
rng=rng,
)
ci = result.confidence_interval
print(ci.low, ci.high)
0.2700347835500123 0.3636602466490255
一方で、このような非対称な分布でも、最初の $t$ 分布による方法が十分頑健であることもわかります。