일표본 평균에 대한 붓스트랩 검정¶

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)

for n in (10, 20, 50, 200):
    M, B = 2000, 999
    s_t = s_stud = s_unstud = 0
    for _ in range(M):
        x = rng.exponential(1, n)
        s_t += stats.ttest_1samp(x, 1).pvalue < 0.05

        xs = x - x.mean() + 1.0          # 중심화
        xb = xs[rng.integers(0, n, (B, n))]

        tobs = (x.mean() - 1) / (x.std(ddof=1) / np.sqrt(n))
        tb = (xb.mean(axis=1) - 1) / (xb.std(axis=1, ddof=1) / np.sqrt(n))
        s_stud += (np.abs(tb) >= abs(tobs)).mean() < 0.05

        uobs = x.mean() - 1
        ub = xb.mean(axis=1) - 1
        s_unstud += (np.abs(ub) >= abs(uobs)).mean() < 0.05
    print(n, round(s_t/M, 3), round(s_stud/M, 3), round(s_unstud/M, 3))

10 0.104 0.074 0.143
20 0.07 0.056 0.091
50 0.069 0.06 0.074
200 0.053 0.048 0.054

import numpy as np
rng = np.random.default_rng(5)
n, B, M = 30, 999, 2000
mu0 = 1.0

p_right, p_wrong = [], []
for _ in range(M):
    x = rng.exponential(1, n)          # H0 가 참이다
    tobs = x.mean() - mu0

    # 올바른 구현: 중심화 후 재표집
    xs = x - x.mean() + mu0
    tb = xs[rng.integers(0, n, (B, n))].mean(axis=1) - mu0
    p_right.append((np.abs(tb) >= abs(tobs)).mean())

    # 잘못된 구현: 중심화 없이 재표집
    bb = x[rng.integers(0, n, (B, n))].mean(axis=1) - mu0
    p_wrong.append((np.abs(bb) >= abs(tobs)).mean())

p_right, p_wrong = np.array(p_right), np.array(p_wrong)
print("올바름: 크기 %.3f, p 중앙값 %.3f" % ((p_right < .05).mean(),
                                            np.median(p_right)))
print("잘못됨: 크기 %.3f, p 중앙값 %.3f" % ((p_wrong < .05).mean(),
                                            np.median(p_wrong)))

올바름: 크기 0.079, p 중앙값 0.453
잘못됨: 크기 0.000, p 중앙값 0.571

import numpy as np
rng = np.random.default_rng(17)
n, B, M = 25, 4000, 1500
mu0 = 1.0
agree = {k: 0 for k in ("unstud~basic", "stud~boott", "stud~pct",
                        "unstud~pct", "stud~basic")}

for _ in range(M):
    x = rng.exponential(1, n)
    th = x.mean(); se = x.std(ddof=1) / np.sqrt(n)
    idx = rng.integers(0, n, (B, n))

    # 두 검정 (중심화된 자료에서)
    xsb = (x - th + mu0)[idx]
    p_un = (np.abs(xsb.mean(axis=1) - mu0) >= abs(th - mu0)).mean()
    tbs = (xsb.mean(axis=1) - mu0) / (xsb.std(axis=1, ddof=1) / np.sqrt(n))
    p_st = (np.abs(tbs) >= abs((th - mu0)/se)).mean()

    # 세 구간 (원자료에서)
    xb = x[idx]; bs = xb.mean(axis=1)
    lo, hi = np.percentile(bs, [2.5, 97.5])
    basic = (2*th - hi, 2*th - lo)
    tt = (bs - th) / (xb.std(axis=1, ddof=1) / np.sqrt(n))
    tl, tu = np.percentile(tt, [2.5, 97.5])
    bt = (th - tu*se, th - tl*se)

    r_un, r_st = p_un < 0.05, p_st < 0.05
    o_b = not (basic[0] <= mu0 <= basic[1])
    o_t = not (bt[0] <= mu0 <= bt[1])
    o_p = not (lo <= mu0 <= hi)

    agree["unstud~basic"] += r_un == o_b
    agree["stud~boott"]   += r_st == o_t
    agree["stud~pct"]     += r_st == o_p
    agree["unstud~pct"]   += r_un == o_p
    agree["stud~basic"]   += r_st == o_b

for k, v in agree.items():
    print(k, round(v/M, 4))

unstud~basic 0.9853
stud~boott 0.972
stud~pct 0.9707
unstud~pct 0.982
stud~basic 0.954

import numpy as np
from scipy import stats
rng = np.random.default_rng(13)
n, B, M = 20, 999, 2500

def pvals(shift, one_sided=False):
    pt, pb = [], []
    alt = 'greater' if one_sided else 'two-sided'
    for _ in range(M):
        x = rng.exponential(1, n) * shift      # 참 평균 = shift
        pt.append(stats.ttest_1samp(x, 1, alternative=alt).pvalue)
        xs = x - x.mean() + 1.0
        xb = xs[rng.integers(0, n, (B, n))]
        tobs = (x.mean() - 1) / (x.std(ddof=1) / np.sqrt(n))
        tb = (xb.mean(axis=1) - 1) / (xb.std(axis=1, ddof=1) / np.sqrt(n))
        pb.append((tb >= tobs).mean() if one_sided
                  else (np.abs(tb) >= abs(tobs)).mean())
    return np.array(pt), np.array(pb)

p0t, p0b = pvals(1.0)
ct, cb = np.percentile(p0t, 5), np.percentile(p0b, 5)   # 크기 보정 임계값
for shift in (2.0, 3.0, 5.0):
    p1t, p1b = pvals(shift)
    print(shift, round((p1t < ct).mean(), 3), round((p1b < cb).mean(), 3))

2.0 0.334 0.318
3.0 0.788 0.652
5.0 0.962 0.82