붓스트랩 신뢰구간 방법의 비교¶

동기¶

앞의 절들에서 백분위수, BCa, 붓스트랩-\(t\)의 세 가지 붓스트랩 신뢰구간 방법을 소개했다. 각각 이론적 성질, 계산 요구사항, 실무적 강점이 다르다. 그중에서 고르려면 이 맞교환을 이해해야 한다. 이 절은 실무자의 판단을 돕기 위한 체계적 비교를 제공한다.

요약표¶

포함확률 정확도¶

가장 중요한 구별은 정확도의 차수이다. 명목 \(100(1-\alpha)\%\) 구간의 실제 포함확률은 다음을 만족한다.

• 백분위수: \(P(\theta \in \text{CI}) = 1 - \alpha + O(n^{-1/2})\)
• BCa: \(P(\theta \in \text{CI}) = 1 - \alpha + O(n^{-1})\)
• 붓스트랩-\(t\): \(P(\theta \in \text{CI}) = 1 - \alpha + O(n^{-1})\)

\(n = 20\)인 95% 구간에서 1차 오차는 수 퍼센트포인트(실제 포함확률 90--93%)일 수 있는 반면, 2차 오차는 보통 1퍼센트포인트 미만(실제 94--96%)이다.

변환 불변성¶

백분위수 구간과 BCa 구간은 변환 불변이다. \([L, U]\)가 \(\theta\)의 구간이면 임의의 단조증가 \(m\)에 대해 \([m(L), m(U)]\)가 \(\phi = m(\theta)\)의 구간이다. 붓스트랩-\(t\) 구간은 이 성질을 갖지 않는다.

변환 불변성이 중요한 경우는 다음과 같다.

• 모수에 자연스러운 경계가 있을 때(예: \(\sigma^2 > 0\), \(0 < p < 1\))
• 추론에 "옳은" 척도가 무엇인지 분명하지 않을 때
• 분석자가 변환된 척도로 결과를 보고할 수 있을 때(예: 확률 대신 로그 오즈)

계산 요구사항¶

백분위수: \(\hat{\theta}\)의 붓스트랩 복제값 \(B\)개만 있으면 된다. 가장 싸다.

BCa: 붓스트랩 복제값 \(B\)개에 더해 가속 \(\hat{a}\)를 위한 \(n\)번의 잭나이프 계산이 필요하다. 총 \(B + n\)번.

붓스트랩-\(t\): 각 붓스트랩 표본마다 \(\hat{\theta}^{*(b)}\)와 \(\hat{\text{se}}^{*(b)}\)를 모두 계산해야 한다. \(\hat{\text{se}}\)의 닫힌 공식이 있으면 비용이 백분위수와 비슷하다. 각 복제마다 잭나이프나 중첩 붓스트랩으로 표준오차를 추정해야 하면 비용이 \(B \times n\)(잭나이프) 또는 \(B \times B_2\)(중첩 붓스트랩)이다.

각 방법을 언제 쓰는가¶

백분위수법¶

• 빠르고 근사적인 구간으로 충분할 때
• 표본크기가 클 때(\(n > 100\))
• 통계량의 표본분포가 거의 대칭임이 알려져 있을 때
• 계산자원이나 구현 복잡도에 제약이 있을 때

BCa 방법¶

• 2차 정확도가 필요할 때(중간 \(n\), 치우친 통계량)
• 변환 불변성이 필요할 때
• 통계량에 간단한 표준오차 공식이 없을 때
• 잭나이프가 계산상 가능할 때(즉 \(n\)이 지나치게 크지 않을 때)

붓스트랩-t 방법¶

• 닫힌 형태의 표준오차 공식이 있을 때(예: 평균)
• 추축량 \((\hat{\theta} - \theta)/\hat{\text{se}}\)의 분포가 \(\theta\)에 거의 무관할 때
• 변환 불변성이 필요 없을 때
• \(\hat{\text{se}}^*\)의 분포가 안정적일 때(\(0\)에 가까운 값이 없을 때)

실무 권고¶

입문 통계학의 대부분 응용에서 다음 결정 규칙이 합리적이다.

1. 기본 선택: BCa --- 2차 정확도, 변환 불변, 폭넓게 적용 가능
2. BCa를 쓸 수 없을 때: 백분위수 --- 단순하고 로버스트하며 큰 \(n\)에 충분
3. 표본평균의 경우: 붓스트랩-\(t\) --- 추가 비용이 거의 없이 2차 정확도
4. 계산이 매우 비싼 통계량: 백분위수 --- BCa의 잭나이프 부담을 피한다

모의실험 보기 1: 잘 작동하는 경우¶

\(\chi^2_3\) 분포(오른쪽으로 치우쳐 있고 참 평균이 \(3\))에서 \(n = 100\)인 표본의 평균을 추정한다. 2000회 반복, \(B = 2000\).

\(n = 100\)이고 통계량이 평균이므로 다섯 방법이 사실상 구별되지 않는다. 포함확률이 \(0.934\)--\(0.947\), 폭이 \(0.946\)--\(0.986\) 범위에 모여 있다.

이 상황에서는 가장 싸고 간단한 백분위수법을 쓰면 된다. BCa의 잭나이프 \(100\)번이나 붓스트랩-\(t\)의 추가 계산이 아무 이득도 주지 않는다.

모의실험 보기 2: 무너지는 경우¶

이제 어려운 문제로 바꾼다. \(\text{Exp}(1)\) 자료에서 \(n = 15\)로 모분산 \(\sigma^2 = 1\)을 추정한다. \(s^2\)의 표본분포는 오른쪽으로 심하게 치우쳐 있다.

3000회 반복, \(B = 1500\):

세 방법 모두 명목 \(0.95\)에 크게 못 미친다. 어떤 붓스트랩 방법도 이 문제를 풀지 못한다.

두 모의실험을 하나의 평면에 올려 놓으면 방법 선택이 언제 문제가 되는지가 한눈에 들어온다. 가로축에 평균 폭을, 세로축에 실제 포함확률을 두고 각 방법을 점 하나로 찍었다.

왼쪽이 쉬운 문제다. 다섯 점이 포함확률 \(0.934\)--\(0.947\), 폭 \(0.946\)--\(0.986\)의 좁은 구역에 모여 있다. 가로축의 눈금이 \(0.94\)에서 \(0.99\)까지라는 점에 주의하라. 점들이 흩어져 보이지만 실제 차이는 폭에서 \(4\%\), 포함확률에서 \(1.3\)%p에 지나지 않으며, 2000회 반복의 몬테카를로 표준오차 \(0.005\)의 두세 배 수준이다. 이 구역에서는 가장 싼 백분위수법을 쓰고 계산을 아끼는 것이 옳다.

오른쪽이 어려운 문제다. 가로축이 로그 눈금으로 바뀌었다는 사실 자체가 상황을 말해 준다. 백분위수 \((1.59,\ 0.677)\)에서 BCa \((1.88,\ 0.736)\)로 가는 이동은 폭을 \(18\%\) 늘려 포함확률을 \(6\)%p 사는 정직한 거래다. 그러나 붓스트랩-\(t\) \((14.15,\ 0.878)\)는 다르다. 포함확률이 가장 높지만 폭이 \(14\)이고, 참값이 \(1\)인 문제에서 폭 \(14\)짜리 구간은 아무것도 말해 주지 않는다. 세로축만 보면 승자이고, 두 축을 함께 보면 탈락이다.

그리고 오른쪽 그림에서 정작 눈에 띄는 것은 세 점이 모두 붉은 명목선 \(0.95\) 아래에 한참 떨어져 있다는 사실이다. 방법을 바꾸어 점을 옮기는 일은 이 그림 안에서의 이동일 뿐, 그림 전체를 위로 올리지는 못한다. 연습문제 2에서 확인하듯 통계량의 어려움이 표본크기보다 결정적이며, \(\text{Exp}(1)\)의 분산을 \(n = 15\)로 추정한다는 설정 자체가 명목수준을 포기하게 만든다.

이런 자료에서는 어떻게 해야 하는가.

1. 변환한다. \(\log s^2\)의 척도에서 구간을 만들고 되돌리면 치우침이 크게 줄어든다.
2. 모수적 방법을 쓴다. 지수분포임을 안다면 \(\sigma^2 = 1/\lambda^2\)의 정확 구간을 쓸 수 있다.
3. 표본을 늘린다. \(n = 15\)로 4차 적률에 의존하는 양을 추정하는 것 자체가 무리이다.

그 밖의 방법¶

문헌에 등장하는 다른 두 방법이 있다.

정규구간 (\(\hat{\theta} \pm z_{\alpha/2} \cdot \widehat{\text{SE}}_{\text{boot}}\)): 1차 정확도이며 표본분포의 대칭성을 가정한다. 백분위수법보다 간단하지만 변환 불변이 아니고 모수 경계를 존중하지 않는다.

기본(추축) 구간 (\([2\hat{\theta} - \hat{\theta}^*_{(1-\alpha/2)}, \; 2\hat{\theta} - \hat{\theta}^*_{(\alpha/2)}]\)): 붓스트랩 분포의 위치 편향을 보정하지만 여전히 1차 정확도이다. 변환 불변이 아니다.

이들은 덜 권장되지만 특정 맥락에서 유용할 수 있다.

연습문제¶

import numpy as np
from scipy import stats
rng = np.random.default_rng(7)
n, B = 100, 10000
x = rng.chisquare(3, n)
th = x.mean(); se = x.std(ddof=1) / np.sqrt(n)

idx = rng.integers(0, n, (B, n)); xb = x[idx]
bs = xb.mean(axis=1)

lo_p, hi_p = np.percentile(bs, [2.5, 97.5])
ci = {
    "normal":  (th - 1.96*bs.std(ddof=1), th + 1.96*bs.std(ddof=1)),
    "pct":     (lo_p, hi_p),
    "basic":   (2*th - hi_p, 2*th - lo_p),
    "t-int":   (th - stats.t.ppf(0.975, n-1)*se, th + stats.t.ppf(0.975, n-1)*se),
}
for k, (lo, hi) in ci.items():
    print(f"{k:7s} [{lo:.4f}, {hi:.4f}]  중심에서 {th-lo:.4f} / {hi-th:.4f}")

normal  [2.1442, 2.9287]  중심에서 0.3923 / 0.3923
pct     [2.1522, 2.9388]  중심에서 0.3842 / 0.4023
basic   [2.1341, 2.9207]  중심에서 0.4023 / 0.3842
t-int   [2.1373, 2.9356]  중심에서 0.3991 / 0.3991

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

def coverage(gen, stat, target, n, M=800, B=600):
    cp = cb = 0
    for _ in range(M):
        x = gen(n); th = stat(x)
        idx = rng.integers(0, n, (B, n))
        bs = np.array([stat(x[i]) for i in idx])
        lo, hi = np.percentile(bs, [2.5, 97.5])
        cp += lo <= target <= hi
        z0 = stats.norm.ppf(np.clip((bs < th).mean(), 1e-6, 1-1e-6))
        jk = np.array([stat(np.delete(x, i)) for i in range(n)])
        d = jk.mean() - jk
        den = ((d**2).sum())**1.5
        a = (d**3).sum()/(6*den) if den > 0 else 0.0
        zl, zu = -1.959964, 1.959964
        a1 = stats.norm.cdf(z0 + (z0+zl)/(1 - a*(z0+zl)))
        a2 = stats.norm.cdf(z0 + (z0+zu)/(1 - a*(z0+zu)))
        lo, hi = np.percentile(bs, [100*a1, 100*a2])
        cb += lo <= target <= hi
    return cp/M, cb/M

ex = lambda n: rng.exponential(1, n)
for n in (15, 50, 200):
    print("mean", n, coverage(ex, np.mean, 1.0, n))
    print("var ", n, coverage(ex, lambda v: v.var(ddof=1), 1.0, n))

mean 15 (0.90375, 0.91125)
var  15 (0.6625, 0.7075)
mean 50 (0.93125, 0.9375)
var  50 (0.81, 0.85)
mean 200 (0.9475, 0.94625)
var  200 (0.87875, 0.9075)

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

def run(n, M=1200, B=1200, trim=None):
    c, w = 0, []
    for _ in range(M):
        x = rng.exponential(1, n)
        xb = x[rng.integers(0, n, (B, n))]
        cc = xb - xb.mean(axis=1, keepdims=True)
        m2 = (cc**2).mean(axis=1); m4 = (cc**4).mean(axis=1)
        bs = xb.var(axis=1, ddof=1)
        seb = np.sqrt(np.maximum((m4 - m2**2)/n, 1e-12))
        co = x - x.mean(); th = x.var(ddof=1)
        se = np.sqrt(max(((co**4).mean() - ((co**2).mean())**2)/n, 1e-12))
        t = (bs - th) / seb
        if trim is not None:
            t = t[np.abs(t) < np.percentile(np.abs(t), 100 - trim)]
        tl, tu = np.percentile(t, [2.5, 97.5])
        lo, hi = th - tu*se, th - tl*se
        c += lo <= 1 <= hi; w.append(hi - lo)
    return round(c/M, 3), round(np.mean(w), 2), round(np.median(w), 2)

for n in (15, 50, 200):
    print(n, run(n))
print("trim 1%", run(15, trim=1))
print("trim 5%", run(15, trim=5))

15 (0.895, 14.83, 4.14)
50 (0.911, 2.74, 1.67)
200 (0.933, 0.96, 0.83)
trim 1% (0.84, 9.34, 2.93)
trim 5% (0.798, 8.25, 2.11)

정리하며¶

백분위수, BCa, 붓스트랩-\(t\) 방법은 정교함이 커지는 위계를 이룬다. 백분위수법은 가장 단순하지만 1차 정확도에 그친다. BCa와 붓스트랩-\(t\)는 서로 다른 방식으로 2차 정확도를 달성한다. BCa는 편향과 가속 보정으로 분위수 수준을 조정하고, 붓스트랩-\(t\)는 통계량을 스튜던트화하여 근사 추축량을 만든다. 일반적인 용도로는 BCa가 권장 기본값이고, 표준오차 공식을 바로 쓸 수 있으면 붓스트랩-\(t\)가 낫다.