재표집 방법 비교 (코드)¶

import numpy as np
from scipy import stats

def bootstrap_ci_demo(data, B=10_000, alpha=0.05, rng=None):
    """평균에 대한 세 가지 붓스트랩 신뢰구간을 함께 구한다.

    정규법은 붓스트랩으로 표준오차만 얻고 구간은 정규분포로 만든다.
    백분위수법은 붓스트랩 분포의 분위점을 그대로 쓴다. 기본법은 그
    분위점을 추정값 둘레로 되비춘다. 분포가 대칭이면 셋이 거의 같다.
    """
    rng = rng or np.random.default_rng(0)
    n = len(data)
    theta_hat = data.mean()
    z = stats.norm.ppf(1 - alpha / 2)

    boot_means = data[rng.integers(0, n, (B, n))].mean(axis=1)
    se_boot = boot_means.std(ddof=1)
    lo_q, hi_q = np.percentile(boot_means, [100*alpha/2, 100*(1 - alpha/2)])

    return {
        "normal":     (theta_hat - z*se_boot, theta_hat + z*se_boot),
        "percentile": (lo_q, hi_q),
        "basic":      (2*theta_hat - hi_q, 2*theta_hat - lo_q),
    }

d = np.random.default_rng(12).exponential(1/3, 30)
r = bootstrap_ci_demo(d)
th = d.mean()
for k, (a, b) in r.items():
    print(f"{k:>11}: [{a:.6f}, {b:.6f}]  폭 {b-a:.6f}  중점 {(a+b)/2:.6f}")
print(f"theta_hat = {th:.6f}")
print(f"중점 합 (기본+백분위수) = {(sum(r['basic'])/2 + sum(r['percentile'])/2):.6f}"
      f"  = 2*theta_hat = {2*th:.6f}")
print(f"정규 구간 폭 극한 = 2 z sigma-hat/sqrt(n) = "
      f"{2*1.959964*d.std(ddof=0)/np.sqrt(30):.6f}")

     normal: [0.257486, 0.511222]  폭 0.253736  중점 0.384354
 percentile: [0.261750, 0.516883]  폭 0.255133  중점 0.389317
      basic: [0.251825, 0.506958]  폭 0.255133  중점 0.379392
theta_hat = 0.384354
중점 합 (기본+백분위수) = 0.768708  = 2*theta_hat = 0.768708
정규 구간 폭 극한 = 2 z sigma-hat/sqrt(n) = 0.255376

def permutation_test_two_sample(x, y, B=9999, stat_func=None, rng=None):
    """이표본 순열검정. 통계량 함수를 바꿔 끼울 수 있다.

    평균 차이든 중앙값 차이든 절사평균 차이든, 귀무가설 아래에서 이름표가
    무의미하다는 논리는 그대로다. 순열검정이 통계량에 매이지 않는 까닭이다.
    """
    rng = rng or np.random.default_rng(0)
    if stat_func is None:
        stat_func = lambda a, b: a.mean() - b.mean()
    t_obs = stat_func(x, y)
    pooled = np.concatenate([x, y])
    m = len(x)

    count = 0
    for _ in range(B):
        p = rng.permutation(pooled)
        count += abs(stat_func(p[:m], p[m:])) >= abs(t_obs)
    return t_obs, (count + 1) / (B + 1)

nx, ny, sx, sy = 20, 50, 1.0, 3.0
true_sd = np.sqrt(sx**2/nx + sy**2/ny)
sbar2 = (nx*sx**2 + ny*sy**2)/(nx+ny)
perm_sd = np.sqrt(sbar2*(nx+ny)/(nx+ny-1)*(1/nx+1/ny))
print(f"참 SD(xbar-ybar) = {true_sd:.4f},  순열 귀무 SD 어림 = {perm_sd:.4f},"
      f"  비 = {perm_sd/true_sd:.4f}")

rng = np.random.default_rng(21)
M, B = 1000, 499
rej_mean = rej_welch = 0
sds = []
for _ in range(M):
    x = rng.normal(5, sx, nx); y = rng.normal(5, sy, ny)
    obs = x.mean() - y.mean(); z = np.concatenate([x, y])
    P = np.array([rng.permutation(z) for _ in range(B)])
    dd = P[:, :nx].mean(1) - P[:, nx:].mean(1)
    sds.append(dd.std(ddof=1))
    rej_mean += ((np.abs(dd) >= abs(obs)).sum() + 1)/(B + 1) < 0.05
    rej_welch += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
print(f"순열 귀무 SD 의 평균(모의) = {np.mean(sds):.4f}")
print(f"제1종 오류율: 평균차 순열 {rej_mean/M:.3f},  Welch t {rej_welch/M:.3f}"
      f"   (몬테카를로 오차 {np.sqrt(0.05*0.95/M):.4f})")

참 SD(xbar-ybar) = 0.4796,  순열 귀무 SD 어림 = 0.6905,  비 = 1.4398
순열 귀무 SD 의 평균(모의) = 0.6840
제1종 오류율: 평균차 순열 0.006,  Welch t 0.044   (몬테카를로 오차 0.0069)

def permutation_test_correlation(x, y, B=9999, rng=None):
    """상관계수에 대한 순열검정.

    한쪽만 섞는다. 그러면 두 변수의 짝은 부서지되 각각의 주변분포는
    그대로 남으므로, "관계가 없다"는 상태를 정확히 흉내 낼 수 있다.
    """
    rng = rng or np.random.default_rng(0)
    r_obs = np.corrcoef(x, y)[0, 1]
    count = 0
    for _ in range(B):
        count += abs(np.corrcoef(x, rng.permutation(y))[0, 1]) >= abs(r_obs)
    return r_obs, (count + 1) / (B + 1)

rng = np.random.default_rng(31)
M, B, n = 1000, 999, 20
rej = 0
ps = []
for _ in range(M):
    x = rng.normal(0, 1, n)
    y = rng.exponential(1, n)          # 독립, 주변분포는 서로 다름
    r = np.corrcoef(x, y)[0, 1]
    perm = np.array([np.corrcoef(x, rng.permutation(y))[0, 1] for _ in range(B)])
    p = ((np.abs(perm) >= abs(r)).sum() + 1) / (B + 1)
    ps.append(p); rej += p < 0.05
ps = np.array(ps)
print(f"제1종 오류율 = {rej/M:.3f}  (몬테카를로 오차 {np.sqrt(0.05*0.95/M):.4f})")
print(f"p 값의 평균 = {ps.mean():.4f} (균등이면 0.5),"
      f"  표준편차 {ps.std(ddof=1):.4f} (균등이면 0.2887)")
print(f"SD(r*) 닫힌 꼴 1/sqrt(n-1) = {1/np.sqrt(n-1):.6f}")

제1종 오류율 = 0.037  (몬테카를로 오차 0.0069)
p 값의 평균 = 0.5104 (균등이면 0.5),  표준편차 0.2904 (균등이면 0.2887)
SD(r*) 닫힌 꼴 1/sqrt(n-1) = 0.229416

def paired_permutation_test(x, y, B=9999, rng=None):
    """차이의 부호를 뒤집는 대응 순열검정.

    대응자료에서는 이름표를 섞으면 안 된다. 짝 자체가 자료의 구조이기
    때문이다. 대신 귀무가설 아래에서 각 차이의 분포가 0 을 중심으로
    대칭이므로, 부호를 아무렇게나 뒤집어도 똑같이 그럴듯하다.
    """
    rng = rng or np.random.default_rng(0)
    d = np.asarray(x) - np.asarray(y)
    t_obs = d.mean()
    signs = rng.choice([-1, 1], size=(B, len(d)))
    t_perm = (signs * d).mean(axis=1)
    return t_obs, ((np.abs(t_perm) >= abs(t_obs)).sum() + 1) / (B + 1)

rng = np.random.default_rng(41)
M, B, n = 2000, 999, 15
for name, gen in [
        ("대칭 (정규)", lambda k: rng.normal(0, 1, k)),
        ("치우침 (지수-1)", lambda k: rng.exponential(1, k) - 1),
        ("강한 치우침 (로그정규 중심화)",
         lambda k: np.exp(rng.normal(0, 1.5, k)) - np.exp(1.5**2/2))]:
    rej = 0
    for _ in range(M):
        d = gen(n); t = d.mean()
        S = rng.choice([-1, 1], size=(B, n))
        tp = (S * d).mean(axis=1)
        rej += ((np.abs(tp) >= abs(t)).sum() + 1)/(B + 1) < 0.05
    print(f"{name:>22}: 제1종 오류율 {rej/M:.3f}"
          f"  (몬테카를로 오차 {np.sqrt(0.05*0.95/M):.4f})")

               대칭 (정규): 제1종 오류율 0.045  (몬테카를로 오차 0.0049)
            치우침 (지수-1): 제1종 오류율 0.088  (몬테카를로 오차 0.0049)
     강한 치우침 (로그정규 중심화): 제1종 오류율 0.248  (몬테카를로 오차 0.0049)

def bootstrap_vs_permutation_comparison(x, y, B=9999, rng=None):
    """붓스트랩 신뢰구간과 순열 p-값을 나란히 놓는다.

    둘은 경쟁 관계가 아니다. 순열검정은 "차이가 있는가"에, 붓스트랩은
    "차이가 얼마나 되는가"에 답한다. 보고할 때는 둘 다 싣는 편이 낫다.
    """
    rng = rng or np.random.default_rng(0)
    diff_obs = x.mean() - y.mean()

    # 각 집단을 따로 재표집한 신뢰구간
    bx = x[rng.integers(0, len(x), (B, len(x)))].mean(axis=1)
    by = y[rng.integers(0, len(y), (B, len(y)))].mean(axis=1)
    ci = np.percentile(bx - by, [2.5, 97.5])

    # 순열 p 값
    _, p_perm = permutation_test_two_sample(x, y, B=B, rng=rng)
    return diff_obs, ci, p_perm

rng = np.random.default_rng(51)
M, B = 500, 999
agree = ci_only = perm_only = 0
for _ in range(M):
    x = rng.normal(0.5, 1, 25); y = rng.normal(0, 1, 25)
    bx = x[rng.integers(0, 25, (B, 25))].mean(axis=1)
    by = y[rng.integers(0, 25, (B, 25))].mean(axis=1)
    ci = np.percentile(bx - by, [2.5, 97.5])
    excl = not (ci[0] <= 0 <= ci[1])

    z = np.concatenate([x, y]); obs = x.mean() - y.mean()
    P = np.array([rng.permutation(z) for _ in range(B)])
    dd = P[:, :25].mean(1) - P[:, 25:].mean(1)
    rej = ((np.abs(dd) >= abs(obs)).sum() + 1)/(B + 1) < 0.05

    agree += (excl == rej)
    ci_only += (excl and not rej)
    perm_only += (rej and not excl)
print(f"{M} 번 중 두 결론이 같은 횟수 = {agree} ({agree/M:.3f})")
print(f"  구간만 기각 = {ci_only} ({ci_only/M:.4f}),"
      f"  순열만 기각 = {perm_only} ({perm_only/M:.4f})")

500 번 중 두 결론이 같은 횟수 = 481 (0.962)
  구간만 기각 = 19 (0.0380),  순열만 기각 = 0 (0.0000)

import numpy as np
from scipy import stats

def cover(n, M=3000, B=2000, scale=3.0, seed=0):
    """Exp(scale) 모집단에서 세 붓스트랩 구간의 경험적 포함확률."""
    rng = np.random.default_rng(seed)
    z = stats.norm.ppf(0.975)
    hits = {"normal": 0, "percentile": 0, "basic": 0}
    for _ in range(M):
        x = rng.exponential(scale, n)
        th = x.mean()
        boot = x[rng.integers(0, n, (B, n))].mean(axis=1)
        se = boot.std(ddof=1)
        lo, hi = np.percentile(boot, [2.5, 97.5])
        hits["normal"] += th - z*se <= scale <= th + z*se
        hits["percentile"] += lo <= scale <= hi
        hits["basic"] += 2*th - hi <= scale <= 2*th - lo
    return {k: round(v / M, 4) for k, v in hits.items()}

for n in (30, 100):
    print(n, cover(n, M=3000, B=2000, scale=3.0))

30 {'normal': 0.91, 'percentile': 0.913, 'basic': 0.9007}
100 {'normal': 0.9387, 'percentile': 0.9383, 'basic': 0.934}

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

def welch_t(a, b):
    return (a.mean() - b.mean()) / np.sqrt(a.var(ddof=1)/len(a)
                                          + b.var(ddof=1)/len(b))

def rate(shift, M=1500, B=999):
    raw = wel = 0
    for _ in range(M):
        x = rng.normal(5 + shift, 1, 20); y = rng.normal(5, 3, 50)
        z = np.concatenate([x, y])
        P = np.array([rng.permutation(z) for _ in range(B)])
        A, Bm = P[:, :20], P[:, 20:]
        d = A.mean(1) - Bm.mean(1); obs = x.mean() - y.mean()
        raw += ((np.abs(d) >= abs(obs)).sum() + 1)/(B+1) < 0.05
        w = d / np.sqrt(A.var(1, ddof=1)/20 + Bm.var(1, ddof=1)/50)
        wo = welch_t(x, y)
        wel += ((np.abs(w) >= abs(wo)).sum() + 1)/(B+1) < 0.05
    return round(raw/M, 3), round(wel/M, 3)

import numpy as np
from itertools import product
from scipy import stats

before = np.array([82, 78, 91, 85, 73, 88, 79, 95, 84, 76])
after  = np.array([88, 82, 95, 89, 78, 91, 84, 98, 90, 81])
d = after - before
print(d)                    # [6 4 4 4 5 3 5 3 6 5]
print(d.mean())             # 4.5

S = np.array(list(product([-1, 1], repeat=len(d))))
t_perm = (S * d).mean(axis=1)
count = (np.abs(t_perm) >= abs(d.mean()) - 1e-12).sum()
print(count, len(t_perm), count / len(t_perm))    # 2  1024  0.001953

[6 4 4 4 5 3 5 3 6 5]
4.5
2 1024 0.001953125

import numpy as np
rng = np.random.default_rng(9)
x = rng.normal(6.2, 1, 20)      # 이동 = 1.2
y = rng.normal(5.0, 3, 50)
diff, ci, p_perm = bootstrap_vs_permutation_comparison(x, y, rng=rng)