재표집 (신발 판매 A/B 검정) (코드)¶

import numpy as np

# 판촉 전후 12주치 신발 판매량. 평균이 2.25 올랐는데, 이것이 판촉의
# 효과인지 그저 주마다의 들쭉날쭉인지를 아래에서 가린다.
BEFORE = np.array([23, 21, 19, 24, 35, 17, 18, 24, 33, 27, 21, 23])
AFTER  = np.array([31, 28, 19, 24, 32, 27, 16, 28, 29, 26, 25, 27])

print(f"Before mean: {BEFORE.mean():.2f}")   # 23.75
print(f"After  mean: {AFTER.mean():.2f}")    # 26.00
print(f"Difference:  {AFTER.mean() - BEFORE.mean():.2f}")   # 2.25

Before mean: 23.75
After  mean: 26.00
Difference:  2.25

n = len(BEFORE)
s1, s2 = BEFORE.std(ddof=1), AFTER.std(ddof=1)
r = np.corrcoef(BEFORE, AFTER)[0, 1]
d = AFTER - BEFORE
print(f"s_전 = {s1:.6f},  s_후 = {s2:.6f},  r = {r:.6f}")
print(f"s_d^2 직접 = {d.var(ddof=1):.6f},  공식 = {s1 ** 2 + s2 ** 2 - 2 * r * s1 * s2:.6f}")

se_paired = d.std(ddof=1) / np.sqrt(n)
se_unpaired = np.sqrt((s1 ** 2 + s2 ** 2) / n)
print(f"대응 SE = {se_paired:.6f},  비대응 SE = {se_unpaired:.6f}")
print(f"비 = {se_paired / se_unpaired:.6f},"
      f"  공식 = {np.sqrt(1 - 2 * r * s1 * s2 / (s1 ** 2 + s2 ** 2)):.6f}")

# 재표집 귀무분포의 표준편차로 견주어도 같은 이야기다.
z = np.concatenate([BEFORE, AFTER])
print(f"섞기 귀무 SD       = {z.std(ddof=1) * np.sqrt(2 / n):.6f}")
print(f"부호뒤집기 귀무 SD = {np.sqrt((d ** 2).sum()) / n:.6f}")

s_전 = 5.561638,  s_후 = 4.612237,  r = 0.609567
s_d^2 직접 = 20.931818,  공식 = 20.931818
대응 SE = 1.320726,  비대응 SE = 2.085756
비 = 0.633212,  공식 = 0.633212
섞기 귀무 SD       = 2.093165
부호뒤집기 귀무 SD = 1.421560

def permutation_test(before, after, n_perm=199_999, rng=None):
    """평균이 올랐는지에 대한 단측 순열검정.

    "올랐는가"만 묻는 것이므로 오른쪽 꼬리만 본다. 양측으로 하면 p-값이
    두 배가 되는데, 어느 쪽으로 할지는 자료를 보기 전에 정해야 한다.
    """
    rng = rng or np.random.default_rng(42)
    observed_diff = after.mean() - before.mean()
    combined = np.concatenate([before, after])
    n_before = len(before)

    P = np.array([rng.permutation(combined) for _ in range(n_perm)])
    perm_diffs = P[:, n_before:].mean(axis=1) - P[:, :n_before].mean(axis=1)
    p_value = ((perm_diffs >= observed_diff - 1e-12).sum() + 1) / (n_perm + 1)
    return observed_diff, perm_diffs, p_value

from collections import defaultdict
from scipy import stats

obs_diff, perm_diffs, p_value = permutation_test(BEFORE, AFTER)
print(f"관측 차이 = {obs_diff:.2f},  순열 단측 p (B=199999) = {p_value:.5f}")
print(f"순열 귀무 SD 모의 = {perm_diffs.std(ddof=1):.6f}")

zi = np.concatenate([BEFORE, AFTER]).astype(int)
total = int(zi.sum())
# k 개를 골랐을 때 합이 s 인 가짓수를 값 하나씩 갱신해 가며 센다.
dp = [defaultdict(int) for _ in range(n + 1)]
dp[0][0] = 1
for v in zi:
    for k in range(n - 1, -1, -1):
        for s, c in list(dp[k].items()):
            dp[k + 1][s + v] += c
counts = dp[n]
n_subsets = sum(counts.values())
hit = sum(c for s, c in counts.items() if s >= AFTER.sum())
lo_ = sum(c for s, c in counts.items() if s <= total - AFTER.sum())

print(f"24개 중 12개를 고르는 방법 = {n_subsets:,}")
print(f"관측된 '후'의 합 = {AFTER.sum()},  전체 합 = {total}")
print(f"합이 {AFTER.sum()} 이상인 배정 = {hit:,}  ->  정확 단측 p = {hit / n_subsets:.6f}")
print(f"정확 양측 p = {(hit + lo_) / n_subsets:.6f}")
print(f"Welch t 검정 (양측) = {stats.ttest_ind(AFTER, BEFORE, equal_var=False).pvalue:.6f}")

관측 차이 = 2.25,  순열 단측 p (B=199999) = 0.15329
순열 귀무 SD 모의 = 2.094198
24개 중 12개를 고르는 방법 = 2,704,156
관측된 '후'의 합 = 312,  전체 합 = 597
합이 312 이상인 배정 = 417,216  ->  정확 단측 p = 0.154287
정확 양측 p = 0.308574
Welch t 검정 (양측) = 0.292782

def bootstrap_ci(before, after, n_boot=100_000, ci=95, rng=None):
    """평균 차이에 대한 백분위수 붓스트랩 신뢰구간.

    검정은 "효과가 있다/없다"까지만 말한다. 판촉에 돈을 쓸지 정하려면
    효과가 얼마나 되는지를 알아야 하고, 그것은 구간이 말해 준다.
    """
    rng = rng or np.random.default_rng(42)
    b = before[rng.integers(0, len(before), (n_boot, len(before)))].mean(axis=1)
    a = after[rng.integers(0, len(after), (n_boot, len(after)))].mean(axis=1)
    diffs = a - b
    lo = (100 - ci) / 2
    return diffs, np.percentile(diffs, [lo, 100 - lo])

sig1, sig2 = BEFORE.std(ddof=0), AFTER.std(ddof=0)
print(f"붓스트랩 SE 의 극한 = {np.sqrt((sig1 ** 2 + sig2 ** 2) / n):.6f}")
print(f"고전 SE             = {se_unpaired:.6f}")
print(f"비 = sqrt((n-1)/n) = {np.sqrt((n - 1) / n):.6f}")

ups, los, ses = [], [], []
for s in range(10):
    diffs_s, ci_s = bootstrap_ci(BEFORE, AFTER, rng=np.random.default_rng(s))
    los.append(ci_s[0]); ups.append(ci_s[1]); ses.append(diffs_s.std(ddof=1))
print(f"\n씨앗 10 개에서")
print(f"  붓스트랩 SE : {min(ses):.4f} ~ {max(ses):.4f}")
print(f"  95% 하한    : {sorted(set(round(v, 4) for v in los))}")
print(f"  95% 상한    : {sorted(set(round(v, 4) for v in ups))}   (격자 간격 {1 / n:.4f})")

diffs, ci95 = bootstrap_ci(BEFORE, AFTER)
_, ci90 = bootstrap_ci(BEFORE, AFTER, ci=90)
print(f"\n기본 씨앗: 90% = [{ci90[0]:.4f}, {ci90[1]:.4f}],"
      f"  95% = [{ci95[0]:.4f}, {ci95[1]:.4f}],  SE = {diffs.std(ddof=1):.6f}")

붓스트랩 SE 의 극한 = 1.996959
고전 SE             = 2.085756
비 = sqrt((n-1)/n) = 0.957427

씨앗 10 개에서
  붓스트랩 SE : 1.9943 ~ 2.0050
  95% 하한    : [-1.8333]
  95% 상한    : [5.9167, 6.0]   (격자 간격 0.0833)

기본 씨앗: 90% = [-1.1667, 5.4167],  95% = [-1.8333, 6.0000],  SE = 1.997389

def permutation_test_two_sided(before, after, n_perm=199_999, rng=None):
    rng = rng or np.random.default_rng(42)
    obs = after.mean() - before.mean()
    combined = np.concatenate([before, after])
    n_b = len(before)
    P = np.array([rng.permutation(combined) for _ in range(n_perm)])
    d = P[:, n_b:].mean(axis=1) - P[:, :n_b].mean(axis=1)
    return obs, ((np.abs(d) >= abs(obs) - 1e-12).sum() + 1) / (n_perm + 1)

import numpy as np, itertools
from scipy import stats

d = AFTER - BEFORE
print(d)          # [ 8  7  0  0 -3 10 -2  4 -4 -1  4  4]
print(d.mean(), d.std(ddof=1))     # 2.25  4.575

S = np.array(list(itertools.product([1, -1], repeat=12)))
m = (S * d).mean(axis=1)
p_one = (m >= d.mean() - 1e-12).sum() / 4096
p_two = (np.abs(m) >= abs(d.mean()) - 1e-12).sum() / 4096
print(p_one, p_two)                # 0.06543  0.13086

[ 8  7  0  0 -3 10 -2  4 -4 -1  4  4]
2.25 4.575130400526108
0.0654296875 0.130859375

theta_hat = AFTER.mean() - BEFORE.mean()      # 2.25
boot_diffs, _ = bootstrap_ci(BEFORE, AFTER, n_boot=100_000)
q_lo, q_hi = np.percentile(boot_diffs, [2.5, 97.5])

print(f"Percentile: [{q_lo:.2f}, {q_hi:.2f}]")
print(f"Basic:      [{2*theta_hat - q_hi:.2f}, {2*theta_hat - q_lo:.2f}]")

Percentile: [-1.83, 6.00]
Basic:      [-1.50, 6.33]

import numpy as np
rng = np.random.default_rng(42)

def power_estimate(n, true_diff=2.0, sd=5.0, M=800, B=999):
    r = 0
    for _ in range(M):
        b = rng.normal(24, sd, n); a = rng.normal(24 + true_diff, sd, n)
        obs = a.mean() - b.mean()
        z = np.concatenate([b, a])
        P = np.array([rng.permutation(z) for _ in range(B)])
        d = P[:, n:].mean(1) - P[:, :n].mean(1)
        r += ((d >= obs - 1e-12).sum() + 1)/(B+1) < 0.05
    return r / M