← 目次 ｜ ← 前: 18 ｜ 次: 20 →

19. リサンプリング法

テーマ: ブートストラップ信頼区間と並べ替え検定（A/Bテストの評価）

学習点: 再標本抽出, ノンパラメトリック推論, p値の意味, 正規近似との比較

依存: NumPy, SciPy ／ 難易度: 中級

実行方法

uv run 19_bootstrap_test.py

スクリプト冒頭の PEP 723 メタデータ(# /// script)により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

EC サイトの施策を対照群 A（180 人）と処置群 B（175 人）で比較する A/B テストを題材に、リサンプリング法による統計的推論を実演するプログラムです。購入金額のような右に大きく歪んだデータ（対数正規分布で生成、歪度 2.5〜3.2）では正規分布の仮定が疑わしいため、分布の形を仮定しないブートストラップ信頼区間と並べ替え検定を使います。

さらに Welch の t 検定・Mann-Whitney の U 検定と p 値を並べて手法間の整合性を確かめ、効果量（Cohen’s d）と検出力の確認まで行う、実務の A/B テスト評価のひな型になっています。

コードの読みどころ

理論的背景

ブートストラップ（Efron 1979）は、手元の標本の経験分布を母集団の代わりに使い、そこからの復元抽出を繰り返して統計量の標本分布を近似する方法です。並べ替え検定は、帰無仮説「2 群の分布が同一」の下では観測値のラベルが交換可能であることを利用し、分布仮定なしに厳密な p 値を与えます。どちらも計算力で数学的仮定を置き換える現代的な推論手法です。

実行結果の見方

平均のブートストラップ標準誤差 194.4（A）が中心極限定理による通常の SE 193.7 とほぼ一致する点は、手法の妥当性の検算になっています。平均差 +832.8 円（+31.9%）に対し、並べ替え検定 p=0.0174、Welch の t 検定 p=0.0182、Mann-Whitney p=0.0205 と 3 手法がそろって 5% 水準で有意になり、結論が特定の仮定に依存していないことがわかります。

ただし Cohen’s d は 0.253 と効果量としては小さめです。最後の検出力の表は、同じ効果でも n=100/群では 12% しか検出できず、n=1000/群でも 64% にとどまることを示しており、「有意でない＝効果がない」と読んではいけない理由（第 2 種の過誤）を数字で確認できます。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "numpy",
#     "scipy",
# ]
# ///
"""19: リサンプリング法 ---------------------------------------------------
テーマ: ブートストラップ信頼区間と並べ替え検定（A/Bテストの評価）
学習点: 再標本抽出, ノンパラメトリック推論, p値の意味, 正規近似との比較
根拠: ブートストラップ(Efron 1979)は経験分布から復元抽出を繰り返して
      統計量の標本分布を近似する。並べ替え検定は帰無仮説「2群の分布が同一」
      の下でラベルが交換可能であることを利用し、分布仮定なしに厳密な
      p 値を与える。
"""
import numpy as np
from scipy import stats

rng = np.random.default_rng(5)


def bootstrap_ci(data, stat_fn, n_boot=20_000, alpha=0.05):
    """パーセンタイル法によるブートストラップ信頼区間。"""
    data = np.asarray(data)
    idx = rng.integers(0, len(data), size=(n_boot, len(data)))
    stats_ = np.apply_along_axis(stat_fn, 1, data[idx])
    lo, hi = np.quantile(stats_, [alpha / 2, 1 - alpha / 2])
    return stats_.mean(), stats_.std(ddof=1), lo, hi


def permutation_test(a, b, n_perm=20_000):
    """2群の平均差に対する並べ替え検定（両側）。"""
    a, b = np.asarray(a), np.asarray(b)
    obs = a.mean() - b.mean()
    pool = np.concatenate([a, b])
    na = len(a)
    diffs = np.empty(n_perm)
    for i in range(n_perm):
        p = rng.permutation(pool)
        diffs[i] = p[:na].mean() - p[na:].mean()
    pval = (np.abs(diffs) >= abs(obs)).mean()
    return obs, pval, diffs


def main() -> None:
    # A/Bテスト: 施策後の一人当たり購入金額（右に歪んだ分布）
    a = rng.lognormal(mean=7.6, sigma=0.9, size=180)     # 対照群
    b = rng.lognormal(mean=7.72, sigma=0.9, size=175)    # 処置群

    print(f"対照群 A: n={len(a)}, 平均 {a.mean():,.1f}円, "
          f"中央値 {np.median(a):,.1f}円, 標準偏差 {a.std(ddof=1):,.1f}")
    print(f"処置群 B: n={len(b)}, 平均 {b.mean():,.1f}円, "
          f"中央値 {np.median(b):,.1f}円, 標準偏差 {b.std(ddof=1):,.1f}")
    print(f"歪度: A={stats.skew(a):.2f}, B={stats.skew(b):.2f} "
          "-> 正規性は疑わしい\n")

    print("[平均のブートストラップ信頼区間]")
    for name, d in [("A", a), ("B", b)]:
        m, se, lo, hi = bootstrap_ci(d, np.mean)
        print(f"  {name}: 平均 {d.mean():>9,.1f}  BS標準誤差 {se:>8,.1f}  "
              f"95%CI [{lo:>9,.1f}, {hi:>9,.1f}]")
    print("  参考: 中心極限定理による通常のSE = "
          f"{a.std(ddof=1)/np.sqrt(len(a)):,.1f} (A)\n")

    print("[中央値のブートストラップCI（解析解が難しい統計量）]")
    m, se, lo, hi = bootstrap_ci(b, np.median)
    print(f"  B の中央値 {np.median(b):,.1f}  95%CI [{lo:,.1f}, {hi:,.1f}]\n")

    print("[平均差の検定 — 3つの手法の比較]")
    obs, pval, diffs = permutation_test(a, b, 20_000)
    t_res = stats.ttest_ind(b, a, equal_var=False)       # Welch の t 検定
    u_res = stats.mannwhitneyu(b, a, alternative="two-sided")
    print(f"  観測された平均差 (B-A) = {b.mean()-a.mean():+,.1f}円 "
          f"({b.mean()/a.mean()-1:+.1%})")
    print(f"  並べ替え検定     : p = {pval:.4f}")
    print(f"  Welchのt検定     : t = {t_res.statistic:.3f}, "
          f"p = {t_res.pvalue:.4f}")
    print(f"  Mann-WhitneyのU  : p = {u_res.pvalue:.4f}")
    print("\n  ※ p値は「効果がない世界でこれ以上に極端な差が出る確率」であり、"
          "\n    「効果がない確率」ではない。効果量と信頼区間を併記すべき。")

    print(f"\n[効果量] Cohen's d = "
          f"{(b.mean()-a.mean()) / np.sqrt((a.var(ddof=1)+b.var(ddof=1))/2):.3f}")

    print("\n[検出力の確認] 同じ効果量で標本サイズを変えたときの有意判定率")
    true_lift = 0.12
    for n_each in [50, 100, 200, 500, 1000]:
        hits = 0
        for _ in range(300):
            x = rng.lognormal(7.6, 0.9, n_each)
            y = rng.lognormal(7.6 + np.log(1 + true_lift) - 0.0, 0.9, n_each)
            hits += stats.ttest_ind(y, x, equal_var=False).pvalue < 0.05
        print(f"  n={n_each:>5}/群: 有意となった割合 = {hits/300:>6.1%}")
    print("  ※ 標本が小さいと真の効果があっても検出できない（第2種の過誤）。")


if __name__ == "__main__":
    main()

実行結果

対照群 A: n=180, 平均 2,606.6円, 中央値 1,685.9円, 標準偏差 2,598.9
処置群 B: n=175, 平均 3,439.4円, 中央値 2,145.5円, 標準偏差 3,867.0
歪度: A=2.50, B=3.18 -> 正規性は疑わしい

[平均のブートストラップ信頼区間]
  A: 平均   2,606.6  BS標準誤差    194.4  95%CI [  2,244.3,   3,004.0]
  B: 平均   3,439.4  BS標準誤差    288.3  95%CI [  2,910.2,   4,036.8]
  参考: 中心極限定理による通常のSE = 193.7 (A)

[中央値のブートストラップCI（解析解が難しい統計量）]
  B の中央値 2,145.5  95%CI [1,884.9, 2,409.5]

[平均差の検定 — 3つの手法の比較]
  観測された平均差 (B-A) = +832.8円 (+31.9%)
  並べ替え検定     : p = 0.0174
  Welchのt検定     : t = 2.375, p = 0.0182
  Mann-WhitneyのU  : p = 0.0205

  ※ p値は「効果がない世界でこれ以上に極端な差が出る確率」であり、
    「効果がない確率」ではない。効果量と信頼区間を併記すべき。

[効果量] Cohen's d = 0.253

[検出力の確認] 同じ効果量で標本サイズを変えたときの有意判定率
  n=   50/群: 有意となった割合 =   7.7%
  n=  100/群: 有意となった割合 =  12.0%
  n=  200/群: 有意となった割合 =  21.7%
  n=  500/群: 有意となった割合 =  39.7%
  n= 1000/群: 有意となった割合 =  64.0%
  ※ 標本が小さいと真の効果があっても検出できない（第2種の過誤）。

← 目次 ｜ ← 前: 18 ｜ 次: 20 →