붓스트랩 가설검정 (코드)¶

(1) 해석적으로. 평행이동은 분포의 모양을 전혀 건드리지 않는다. 중심화한 경험분포의 중심적률이

로, 둘째와 셋째는 원자료의 것과 같다. 거기서 독립으로 \(n\)개 뽑아 평균을 내므로

이다. 중심이동이 바꾸는 것은 중심 하나뿐이고, 폭과 치우침은 그대로다. 바로 그래서 이 한 줄이 "귀무가설만 참인 세계"를 만든다.

귀무분포가 정규라면 \(p \approx 2\Phi(-|z|)\), \(z = \dfrac{\bar x - \mu_0}{\hat\sigma/\sqrt n}\)이다. 다만 왜도가 \(\hat\gamma_1/\sqrt n\)만큼 남아 있으므로 오른쪽 꼬리가 정규보다 두껍고, 실제 \(p\)값은 정규근사보다 클 것이다.

(2) 수치적으로. 함수는 이렇다.

import numpy as np

def bootstrap_mean_test(data, mu_0=0, n_boot=10_000, rng=None):
    """H0: 평균 = mu_0 에 대한 붓스트랩 검정. 양측 p-값을 돌려준다.

    검정을 하려면 귀무가설이 참인 상태에서 뽑아야 한다. 그래서 자료의
    중심을 mu_0 으로 옮긴 뒤 재표집한다. 이 중심 이동 한 줄이 신뢰구간을
    만들 때와 검정을 할 때를 가르는 핵심이다.
    """
    rng = rng or np.random.default_rng(0)
    n = len(data)
    # 자료를 통째로 옮겨 평균이 정확히 mu_0 이 되게 한다. 퍼짐과 모양은
    # 그대로 남으므로, 귀무가설만 참인 세계를 흉내 낸 셈이다.
    centered = data - data.mean() + mu_0
    boot_means = centered[rng.integers(0, n, (n_boot, n))].mean(axis=1)
    obs_mean = data.mean()
    # 분자와 분모에 1 을 더한다. 이래야 p-값이 0 이 되는 일을 막을 수 있다.
    # 관측된 자료 자체도 하나의 가능한 재표본으로 세는 셈이다.
    p_value = ((np.abs(boot_means - mu_0) >= abs(obs_mean - mu_0)).sum() + 1) \
              / (n_boot + 1)
    return obs_mean, p_value, boot_means

보기 4의 자료로 돌려 본다.

from scipy import stats

rng_data = np.random.default_rng(1)
rng = np.random.default_rng(7)
data = rng_data.exponential(scale=5, size=50) + 2
obs, p, boots = bootstrap_mean_test(data, mu_0=5.0, rng=rng)

n = len(data)
sig = data.std(ddof=0)
print(f"sigma-hat = {sig:.4f},  예측 귀무 SD = {sig/np.sqrt(n):.4f}"
      f"   (모의 {boots.std(ddof=1):.4f},  몬테카를로 오차 {sig/np.sqrt(n)/np.sqrt(2*10000):.4f})")
print(f"귀무분포 평균 = {boots.mean():.4f}  (예측 mu_0 = 5)")
print(f"귀무분포 왜도 = {stats.skew(boots):.4f}  (예측 g1/sqrt(n) = {stats.skew(data)/np.sqrt(n):.4f})")
z = (obs - 5) / (sig / np.sqrt(n))
print(f"z = {z:.4f},  정규근사 양측 p = {2*stats.norm.sf(z):.4f}")
print(f"붓스트랩 p = {p:.4f},   t 검정 p = {stats.ttest_1samp(data, 5.0).pvalue:.4f}")

rng9 = np.random.default_rng(7)
bm_nc = data[rng9.integers(0, n, (10000, n))].mean(axis=1)
print(f"중심이동을 빼먹으면 p = {((np.abs(bm_nc-5)>=abs(obs-5)).sum()+1)/10001:.4f}")

출력:

sigma-hat = 7.2062,  예측 귀무 SD = 1.0191   (모의 1.0270,  몬테카를로 오차 0.0072)
귀무분포 평균 = 5.0088  (예측 mu_0 = 5)
귀무분포 왜도 = 0.3969  (예측 g1/sqrt(n) = 0.4257)
z = 3.1190,  정규근사 양측 p = 0.0018
붓스트랩 p = 0.0032,   t 검정 p = 0.0033
중심이동을 빼먹으면 p = 0.4741

세 적률이 모두 예측대로다. 귀무분포의 평균이 \(5.0088\)로 \(\mu_0 = 5\)에 붙고(평균의 몬테카를로 오차 \(0.0103\)), 표준편차 \(1.0270\)이 예측 \(1.0191\)에서 \(1.1\) 몬테카를로 오차 떨어져 있으며, 왜도 \(0.3969\)도 예측 \(0.4257\) 근처다(\(\sqrt{6/B} = 0.0245\)의 \(1.2\)배).

정규근사보다 붓스트랩 \(p\)가 크다. \(0.0018\) 대 \(0.0032\)로 거의 두 배다. (1)에서 예상한 대로 귀무분포의 오른쪽 꼬리가 정규보다 두껍기 때문이며, \(t\) 검정의 \(0.0033\)과는 오히려 잘 맞는다. \(t\) 분포도 정규보다 꼬리가 두껍다는 점에서 같은 방향의 보정을 하고 있는 셈이다.

중심이동을 빼먹으면 \(p = 0.4741\)이 된다. 옮기지 않은 붓스트랩 분포는 \(\bar x = 8.179\)에 중심을 두므로 \(\mu_0 = 5\)에서 잰 거리가 복제값들의 거리와 비슷해지고, \(p\)값이 \(0.5\) 근처에 머문다. 어떤 자료를 넣어도 기각하지 못한다. 오류 없이 조용히 실패하는 종류의 버그다.

(1) 해석적으로. 합친 \(N = m + n\)개에서 복원으로 \(N\)개를 뽑아 앞의 \(m\)개를 \(X\), 나머지를 \(Y\)라 부른다. 뽑기가 독립이므로 두 집단의 평균도 독립이고, 각각의 분산이 \(\hat\sigma_p^2/m\)과 \(\hat\sigma_p^2/n\)이다. 따라서

이다(\(\hat\sigma_p^2 = \frac1N\sum_i (z_i - \bar z)^2\), \(z\)는 합친 자료).

순열검정과의 차이는 딱 하나다. 이표본 순열검정 보기 1에서 라벨을 비복원으로 섞을 때의 분산이 \(S^2(1/m + 1/n)\)임을 보았고, 여기서 \(S^2\)은 \(N-1\)로 나눈 표본분산이다. \(\hat\sigma_p^2 = \frac{N-1}{N}S^2\)이므로

다. \(N = 80\)이면 \(0.99373\)으로 \(0.6\%\) 좁다. 복원추출이 유한모집단 보정을 잃는 대신 집단의 크기를 고정하지 않기 때문이며, \(N\)이 크면 사라지는 차이다.

(2) 수치적으로. 함수는 이렇다.

def bootstrap_two_sample(x, y, n_boot=10_000, rng=None):
    """H0: 두 평균이 같다에 대한 붓스트랩 검정.

    귀무가설이 참이면 두 집단이 같은 모집단에서 나온 것이다. 그래서 둘을
    합쳐 하나의 웅덩이로 만들고 거기서 다시 뽑는다.
    """
    rng = rng or np.random.default_rng(0)
    obs_diff = x.mean() - y.mean()
    pooled = np.concatenate([x, y])
    m, N = len(x), len(pooled)
    P = pooled[rng.integers(0, N, (n_boot, N))]
    boot_diffs = P[:, :m].mean(axis=1) - P[:, m:].mean(axis=1)
    p_value = ((np.abs(boot_diffs) >= abs(obs_diff)).sum() + 1) / (n_boot + 1)
    return obs_diff, p_value, boot_diffs

보기 1의 블록에서 이어 쓴다.

x = rng_data.normal(52, 10, 40)
y = rng_data.normal(48, 10, 40)
diff, p2, boots2 = bootstrap_two_sample(x, y, rng=rng)

pooled = np.concatenate([x, y]); N = len(pooled); m = len(x)
sp = pooled.std(ddof=0)
print(f"합친 sigma-hat = {sp:.4f}")
print(f"예측 귀무 SD = sigma-hat*sqrt(1/m+1/n) = {sp*np.sqrt(1/m+1/(N-m)):.4f}"
      f"   (모의 {boots2.std(ddof=1):.4f},  몬테카를로 오차 {sp*np.sqrt(2/m)/np.sqrt(2*10000):.4f})")
print(f"순열(비복원) 판본의 SD = S*sqrt(2/m) = {pooled.std(ddof=1)*np.sqrt(2/m):.4f}"
      f"   비 = {sp/pooled.std(ddof=1):.6f}  (sqrt((N-1)/N) = {np.sqrt((N-1)/N):.6f})")
z2 = diff / (sp*np.sqrt(2/m))
print(f"z = {z2:.4f},  정규근사 양측 p = {2*stats.norm.sf(z2):.4f}")
w = stats.ttest_ind(x, y, equal_var=False)
print(f"붓스트랩 p = {p2:.4f},   Welch t p = {w.pvalue:.4f}")

출력:

합친 sigma-hat = 9.4434
예측 귀무 SD = sigma-hat*sqrt(1/m+1/n) = 2.1116   (모의 2.1242,  몬테카를로 오차 0.0149)
순열(비복원) 판본의 SD = S*sqrt(2/m) = 2.1249   비 = 0.993730  (sqrt((N-1)/N) = 0.993730)
z = 1.9657,  정규근사 양측 p = 0.0493
붓스트랩 p = 0.0518,   Welch t p = 0.0504

유도한 두 식이 맞는다. 예측 \(2.1116\)에 대해 모의값 \(2.1242\)가 몬테카를로 오차 \(0.0149\)의 \(0.8\)배 안에 있고, 순열 판본과의 비 \(0.993730\)이 \(\sqrt{79/80}\)과 여섯 자리까지 같다.

세 \(p\)값이 \(0.05\) 양쪽에 걸쳐 있다. 정규근사 \(0.0493\)은 턱걸이로 기각하고, 붓스트랩 \(0.0518\)과 Welch \(0.0504\)는 기각하지 못한다. 셋의 차이는 \(0.0025\)에 지나지 않지만 \(\alpha = 0.05\)라는 선 위에 걸려 결론이 갈린다. 붓스트랩 \(p\)값 자체의 몬테카를로 표준오차가 \(\sqrt{0.05 \times 0.95/10^4} = 0.0022\)이므로, 이 자료에서 "유의한가"를 묻는 것은 의미가 없다. 구간을 보고하거나 자료를 더 모으는 것이 옳은 대응이다.

(1) 해석적으로. 로그정규 밀도는

이다. 표본중앙값 \(\tilde x = 31{,}508.1\)을 넣으면 \(f(\tilde x) = 1.557962\times10^{-5}\)이고

이다. 모집단 중앙값 \(e^{10.5} = 36{,}315.5\) 자리의 값 \(2{,}574.7\)과 다르다는 데 주의할 것. 표본중앙값이 참값보다 작게 나왔고, 로그정규는 그 언저리에서 밀도가 더 높으므로 표준오차가 작아진다.

편향의 부호. 붓스트랩 편향은 \(\overline{\tilde x^{*}} - \tilde x\), 곧 붓스트랩 중앙값 분포의 평균과 중앙값의 차이에 가깝다. 모집단이 오른쪽으로 치우쳐 있으면 재표본 중앙값의 분포도 오른쪽으로 치우치고, 치우친 분포에서는 평균이 중앙값보다 크다. 그러므로 편향은 양수일 것이다.

크기 판정. 편향은 절대 크기가 아니라 표준오차에 견준 크기로 본다. \(\lvert\widehat{\text{bias}}\rvert/\widehat{\operatorname{SE}} < 0.25\)이면 무시할 만하다는 것이 흔한 경험칙이다.

(2) 수치적으로. 함수는 이렇다.

def bootstrap_se_median(data, n_boot=10_000, rng=None):
    """중앙값의 표준오차와 편향을 붓스트랩으로 추정한다.

    편향은 붓스트랩 값들의 평균에서 원래 추정값을 뺀 것이다. 0 에서 멀면
    그 통계량이 참값을 체계적으로 빗나간다는 뜻이다.
    """
    rng = rng or np.random.default_rng(0)
    n = len(data)
    boot_medians = np.median(data[rng.integers(0, n, (n_boot, n))], axis=1)
    se = boot_medians.std(ddof=1)
    bias = boot_medians.mean() - np.median(data)
    return se, bias, boot_medians

보기 2의 블록에서 이어 쓴다.

income = rng_data.lognormal(mean=10.5, sigma=0.8, size=200)
se_med, bias, boot_med = bootstrap_se_median(income, rng=rng)

mu_, sg_ = 10.5, 0.8
m_s = np.median(income)
f_at = np.exp(-(np.log(m_s) - mu_) ** 2 / (2 * sg_ ** 2)) / (m_s * sg_ * np.sqrt(2 * np.pi))
print(f"표본중앙값 = {m_s:,.1f},  모집단 중앙값 = {np.exp(mu_):,.1f}")
print(f"표본중앙값 자리의 참 밀도 = {f_at:.6e}  ->  SE = {1/(2*f_at*np.sqrt(200)):,.1f}")
print(f"붓스트랩 SE = {se_med:,.1f},  편향 = {bias:,.1f},  |편향|/SE = {abs(bias)/se_med:.4f}")
print(f"붓스트랩 중앙값의 왜도 = {stats.skew(boot_med):.4f},"
      f"  고유값 수 = {len(np.unique(boot_med))}")

출력:

표본중앙값 = 31,508.1,  모집단 중앙값 = 36,315.5
표본중앙값 자리의 참 밀도 = 1.557962e-05  ->  SE = 2,269.3
붓스트랩 SE = 1,912.2,  편향 = 269.2,  |편향|/SE = 0.1408
붓스트랩 중앙값의 왜도 = 0.2382,  고유값 수 = 263

편향의 부호를 맞혔다. \(+269.2\)로 양수이고, 붓스트랩 중앙값 분포의 왜도가 \(+0.2382\)인 것이 그 까닭이다. \(\lvert\text{편향}\rvert/\widehat{\operatorname{SE}} = 0.1408\)로 경험칙 \(0.25\)를 밑돈다. 보정할 만큼 크지 않다.

표준오차는 \(1{,}912.2\)로 이론값 \(2{,}269.3\)보다 \(16\%\) 작다. \(B = 10{,}000\)의 몬테카를로 요동이 \(1912/\sqrt{2B} = 13.5\)뿐이므로 \(B\) 탓이 아니다. 까닭은 중앙값의 붓스트랩 보기 1에서 본 그대로다. 붓스트랩이 쓰는 것은 참 밀도가 아니라 표본중앙값 둘레의 관측 간격이고, 이 표본은 그 자리가 모집단보다 촘촘했다. 중앙값의 표준오차를 붓스트랩으로 구할 때 한 자리 이상을 믿어서는 안 되는 이유다.

값이 \(263\)가지뿐이라는 것도 눈여겨볼 것. \(10{,}000\)개의 복제값이 서로 다른 값을 그만큼만 갖는다. 재표본 중앙값이 원자료 순서통계량 두 개의 평균일 수밖에 없기 때문이며, 그래서 중앙값의 붓스트랩 분포는 계단 모양이 된다.

(1) 해석적으로. 몬테카를로 오차의 꼴이 양마다 다르다.

둘째 줄이 문제가 된다. \(p = 0.0518\)에 \(\pm 0.0022\)가 붙으므로 \(\alpha = 0.05\)라는 선을 사이에 두고 어느 쪽인지 말할 수 없다. 첫째 줄은 \(0.0032 \pm 0.00056\)으로 결론이 흔들리지 않는다.

고전적 비교 대상이 없는 것은 셋째다. 중앙값의 표준오차는 \(\dfrac{1}{2f(m)\sqrt n}\)인데 이 식에는 모집단 밀도 \(f(m)\)이 들어 있다. 평균의 \(\sigma/\sqrt n\)에서 \(\sigma\)를 \(s\)로 바꾸듯 \(f(m)\)을 자료에서 바로 바꿔 넣을 수가 없다. 밀도를 추정하려면 띠너비를 골라야 하고 그 선택이 답을 바꾼다. 붓스트랩은 그 선택을 하지 않고 같은 일을 해낸다. 게다가 셋째는 애초에 검정이 아니라 추정의 불확실성을 재는 일이라 \(p\)값이 없다.

(2) 수치적으로.

import numpy as np
from scipy import stats

rng_data = np.random.default_rng(1)
rng = np.random.default_rng(7)

# 1. 일표본 검정: 2 만큼 이동된 지수 자료
data = rng_data.exponential(scale=5, size=50) + 2
obs, p, boots = bootstrap_mean_test(data, mu_0=5.0, rng=rng)
print(obs, p)                    # 8.1786  0.0032

# 2. 이표본 검정: 두 정규 모집단
x = rng_data.normal(52, 10, 40)
y = rng_data.normal(48, 10, 40)
diff, p2, boots2 = bootstrap_two_sample(x, y, rng=rng)
print(diff, p2)                  # 4.151  0.0518

# 3. 중앙값의 붓스트랩 표준오차: 로그정규 소득
income = rng_data.lognormal(mean=10.5, sigma=0.8, size=200)
se_med, bias, boot_med = bootstrap_se_median(income, rng=rng)
print(np.median(income), se_med, bias)   # 31508  1912  269

출력:

8.178614426992942 0.0031996800319968005
4.150721066684774 0.051794820517948204
31508.06500959787 1912.1800372703474 269.241019752415

두 \(p\)값이 \(t\) 검정과 거의 같다.

일표본 쪽은 사실상 일치한다. \(0.0032\) 대 \(0.0033\)으로 차이가 몬테카를로 오차 \(0.0006\)의 \(0.2\)배다. \(n = 50\)에서 정규근사가 이미 잘 통하고, 보기 1에서 보았듯 붓스트랩 귀무분포의 왜도가 \(t_{49}\)의 두꺼운 꼬리와 비슷한 보정을 하기 때문이다.

이표본 쪽은 \(0.0518\) 대 \(0.0504\)로 \(\alpha = 0.05\)의 양쪽에 걸쳐 있지만 구별할 수 없다. 차이 \(0.0014\)가 몬테카를로 오차 \(0.0022\)보다 작다. "유의하다/아니다"를 가르는 대신 차이의 구간을 보고해야 할 자료다.

셋째는 비교할 것이 없다. 그것이 이 절차를 쓰는 이유다.

import numpy as np
from scipy import stats
rng = np.random.default_rng(101)
data = np.random.default_rng(3).normal(0, 1, 40)
print(data.mean(), data.std(ddof=1))       # -0.0411  1.1751

_, p0, _ = bootstrap_mean_test(data, mu_0=0.0, rng=rng)
_, p05, _ = bootstrap_mean_test(data, mu_0=0.5, rng=rng)

출력:

-0.041112513627908596 1.1751281074452444

참 평균이 \(0\)이므로 \(H_0\colon \mu = 0\)에서는 큰 \(p\)값이 나와 기각하지 못한다. \(H_0\colon \mu = 0.5\)에서는 표본평균 \(-0.041\)이 \(0.5\)에서 \(2.9\) 표준오차 떨어져 있어(\(\widehat{\text{se}} = 1.175/\sqrt{40} = 0.186\)) 강하게 기각한다.

중심화가 무엇을 하는지 정확히 보자. \(\mu_0 = 0.5\)일 때 중심화된 자료는 \(x_i - (-0.041) + 0.5 = x_i + 0.541\)이다. 이 자료의 평균은 정확히 \(0.5\)이고, 재표집된 평균들은 \(0.5\) 주위에 표준편차 \(0.186\)으로 흩어진다.

관측된 편차 \(|{-0.041} - 0.5| = 0.541\)은 이 분포에서 \(2.9\) 표준편차에 해당하므로 그만큼 극단적인 재표본이 거의 나오지 않는다.

붓스트랩 \(p\)값이 \(t\) 검정의 절반이다(\(0.0032\) 대 \(0.0059\)). 이는 붓스트랩이 정규이론의 \(t\) 보정을 하지 않기 때문이다. \(n = 40\)에서 \(t_{39}\)의 꼬리가 정규분포보다 두꺼운 만큼 차이가 난다. 꼬리로 갈수록 이 차이가 커지므로, 작은 \(p\)값을 붓스트랩으로 보고할 때는 주의해야 한다.

이표본 붓스트랩 검정에서는 합친 집합에서 \(m + n\)개를 복원추출한 뒤 크기 \(m\)과 \(n\)으로 나눈다. 따라서 어떤 관측은 여러 번 나타나고 어떤 관측은 그 반복에서 아예 빠진다.

순열검정에서는 \(m + n\)개의 라벨을 비복원으로 섞으므로 모든 관측이 각 순열된 자료에 정확히 한 번 나타난다. 합친 표본의 구성이 정확히 보존된다.

비슷한 \(p\)값을 주는 조건은 표본크기 \(m\)과 \(n\)이 어느 정도 클 때이다. \(m, n \to \infty\)에서 평균차의 붓스트랩 분포가 순열분포로 수렴한다. 작은 표본에서는 순열검정이 (자료에 조건부로) 정확한 반면 붓스트랩 검정은 근사이다.

구체적인 차이. 붓스트랩 분산과 순열 분산의 비는 대략

이다(\(N = m + n\)). 순열은 유한모집단 수정계수를 자동으로 반영하지만 복원추출은 그렇지 않기 때문이다. \(N = 10\)이면 \(11\)%, \(N = 80\)이면 \(1.3\)% 차이이다.

이 차이가 실제로 얼마나 큰지는 비교 연습문제 1에서 확인했다. \(m = n = 5\)일 때 붓스트랩 꼬리 확률 \(0.0006\)과 정확 순열 \(p\)값 \(0.0397\)이 \(66\)배 차이가 났다. 위 비율만으로는 설명되지 않는 차이이며, 작은 표본에서 붓스트랩이 \(t\) 보정을 놓치는 것이 더 큰 원인이다.

순열검정은 무작위화 메커니즘을 직접 모형화하므로 교환가능성이라는 귀무가설 아래에서 더 자연스럽기도 하다. \(\square\)

붓스트랩 편향의 정의는

이다. \(E^*\)는 붓스트랩 분포(자료의 경험분포) 아래의 기댓값이고 \(\overline{\tilde x^*} = \frac{1}{B}\sum_{b=1}^{B}\tilde x^{*(b)}\)가 이를 추정한다.

중앙값의 붓스트랩 분포가 \(\tilde x\)에 대해 대칭이면, \(\tilde x^{*(b)} = \tilde x + \delta\)를 내는 반복마다 (근사적으로) \(\tilde x - \delta\)를 내는 반복이 대응한다. 따라서 \(E^*[\tilde x^*] = \tilde x\)이고

이다. \(\square\)

실제로는 대칭이 아니다. 시연의 로그정규 소득 자료에서 붓스트랩 편향이 \(+269\)였다. 표본중앙값 \(31{,}508\)의 \(0.85\)%이다.

양의 편향이 나오는 이유는 두 가지가 겹친다.

첫째, 중앙값의 붓스트랩 분포는 이산적이다. 붓스트랩 중앙값은 원자료의 값(또는 짝수 개일 때 인접 두 값의 평균)만 취할 수 있다. \(n = 200\)이면 가능한 값이 몇 백 개뿐이다.

둘째, 모분포가 오른쪽으로 치우쳐 있다. 표본중앙값 위쪽의 관측들이 아래쪽보다 더 넓게 퍼져 있으므로, 붓스트랩 중앙값이 위로 움직일 여지가 아래로 움직일 여지보다 크다.

import numpy as np
from scipy import stats
rng = np.random.default_rng(101)
data = np.random.default_rng(12).gamma(shape=2, scale=1, size=200)
print(data.mean(), np.median(data))       # 1.9889  1.7680

boot_means = data[rng.integers(0, 200, (10_000, 200))].mean(axis=1)
se_mean_boot = boot_means.std(ddof=1)
se_med, _, _ = bootstrap_se_median(data, rng=rng)

출력:

1.9888641137690775 1.7679603646706958

붓스트랩 \(\widehat{\text{SE}}\)(평균)이 \(s/\sqrt{n}\)과 소수 넷째 자리까지 일치한다(\(0.0861\) 대 \(0.0860\)). 이는 우연이 아니다. 평균의 붓스트랩 분산은 정확히 \(\hat{\sigma}^2/n = \frac{n-1}{n}s^2/n\)이므로 두 값이 이론적으로 같아야 한다.

이론값 \(0.1000\)과는 \(14\)% 차이가 난다. 이는 붓스트랩의 오차가 아니라 이 표본의 \(s = 1.216\)이 참값 \(\sigma = \sqrt{2} = 1.414\)보다 작기 때문이다. 붓스트랩은 참 \(\sigma\)를 알 수 없으므로 표본에서 추정할 수밖에 없다.

중앙값에 공식이 없는 이유. 중앙값의 점근 표준오차는

이다. 여기서 \(f\)는 모집단 밀도, \(m\)은 모중앙값이다. \(\text{Gamma}(2,1)\)에서 \(m = 1.678\), \(f(m) = 0.313\)이므로 \(1/(2 \times 0.313 \times 14.14) = 0.1128\)이다.

문제는 이 공식이 모집단 밀도를 중앙값 한 점에서 알아야 한다는 것이다. 평균의 공식이 \(\sigma\)만 필요한 것과 대조적이다. \(f(m)\)은 밀도추정을 해야 하고, 그 자체가 대역폭 선택 등의 문제를 안고 있으며 수렴이 느리다(\(n^{-2/5}\)).

붓스트랩은 이 문제를 완전히 우회한다. np.median을 반복 계산할 뿐이며 밀도를 추정할 필요가 없다. 위 표에서 붓스트랩 \(0.1202\)가 점근값 \(0.1128\)과 \(7\)% 이내로 일치한다.

대립가설 \(\mu \neq \mu_0\) 아래에서 큰 수의 법칙에 의해 \(n \to \infty\)일 때 \(\bar x \to \mu\)이므로

이다.

중심화된 자료 \(x_i^0 = x_i - \bar x + \mu_0\)의 표본평균은 정확히 \(\mu_0\)이다. 붓스트랩 중심극한정리에 의해 중심화된 자료에서 재표집한 평균 \(\bar x^{*(b)}\)는

를 만족한다. 여기서 \(\sigma^2\)은 모분산이다. 따라서 \(|\bar x^{*(b)} - \mu_0| = O_p(n^{-1/2})\)이고 이는 \(0\)으로 수렴한다.

한편 \(|\bar x - \mu_0| \to |\mu - \mu_0| > 0\)은 고정된 양의 상수이다. \(n\)이 크면 붓스트랩 반복이 \(|\bar x^{*(b)} - \mu_0| \ge |\bar x - \mu_0|\)를 만족할 확률이 무시할 수준이 된다.

따라서 검정은 일치한다. \(\square\)

수렴 속도. 더 정확히는, \(\bar{x}^* - \mu_0\)가 근사적으로 \(N(0, \sigma^2/n)\)이므로

이고, 이는 \(n\)에 대해 지수적으로 \(0\)에 접근한다. \(|\mu - \mu_0|/\sigma = 0.5\)일 때

실용적 함의가 하나 있다. \(B\)가 유한하면 이 수렴을 관측할 수 없다. \(B = 10{,}000\)이면 붓스트랩이 낼 수 있는 최소 \(p\)값이 \(1/10001 = 0.0001\)이다. \(n = 100\)의 참 \(p\)값 \(5.7 \times 10^{-7}\)은 표현조차 되지 않는다.

작은 \(p\)값을 정량적으로 보고해야 한다면 붓스트랩 대신 점근 근사를 쓰는 것이 옳다. 붓스트랩의 강점은 \(p\)값의 정밀도가 아니라 가정으로부터의 자유이다.