백분위수법¶

백분위수 붓스트랩 절차:

1. 중앙값의 붓스트랩 복제값 \(B = 10{,}000\)개를 생성한다.
2. 복제값을 정렬한다.
3. 95% 신뢰구간은 \([\hat{\theta}^*_{[250]}, \hat{\theta}^*_{[9750]}]\)이다.

정렬된 붓스트랩 중앙값이 \(\hat{\theta}^*_{[250]} = 11.8\), \(\hat{\theta}^*_{[9750]} = 16.5\)였다고 하자. 95% 백분위수 구간은 \([11.8, 16.5]\)이다.

이 구간이 관측된 중앙값 \(14.3\)을 중심으로 비대칭임에 주목하라. 하한까지의 거리 \(2.5\)가 상한까지의 거리 \(2.2\)와 다르며, 이는 표본분포의 치우침을 반영한 것이다.

import numpy as np
rng = np.random.default_rng(0)
n = 40
x = rng.lognormal(0, 1, n)
B = 20000
idx = rng.integers(0, n, (B, n))
m = x[idx].mean(axis=1)          # 원 척도 붓스트랩 복제값
lm = np.log(m)                    # 로그 척도 복제값

# 백분위수 구간
print(np.round(np.percentile(m, [2.5, 97.5]), 4))
print(np.round(np.exp(np.percentile(lm, [2.5, 97.5])), 4))

# 정규근사 구간
th = x.mean()
print(np.round([th - 1.96*m.std(ddof=1), th + 1.96*m.std(ddof=1)], 4))
lth = np.log(th); lse = lm.std(ddof=1)
print(np.round(np.exp([lth - 1.96*lse, lth + 1.96*lse]), 4))

출력:

[0.9639 1.591 ]
[0.9639 1.591 ]
[0.9468 1.5776]
[0.9816 1.6231]

백분위수 구간이 소수점 이하 전부 일치한다. 우연이 아니라 정의에서 따라 나오는 항등식이다. \(\log\)가 단조증가이므로 \(\log(m)\)의 \(q\)번째 분위수는 \(m\)의 \(q\)번째 분위수의 로그와 정확히 같다.

정규근사 구간은 척도에 따라 다르다. 원 척도에서 \([0.947, 1.578]\), 로그 척도를 거치면 \([0.982, 1.623]\)이다. 하한이 4% 차이 난다.

어느 쪽이 옳은가? 로그 척도를 거친 쪽이 낫다. 로그정규 평균의 표본분포가 오른쪽으로 치우쳐 있어 로그를 취하면 정규에 가까워지기 때문이다. 그러나 이를 알려면 사용자가 올바른 변환을 미리 알아야 한다.

백분위수법의 가치가 여기에 있다. 어떤 변환이 정규화 변환인지 몰라도, 그 변환을 적용한 결과와 같은 답을 준다.

import numpy as np
rng = np.random.default_rng(3)
for n in (25, 100, 400):
    M, B = 1500, 1000
    cn = cp = cb = 0
    for _ in range(M):
        x = rng.exponential(1, n)
        th = x.var(ddof=1)
        idx = rng.integers(0, n, (B, n))
        bs = x[idx].var(axis=1, ddof=1)
        se = bs.std(ddof=1)
        cn += th - 1.96*se <= 1.0 <= th + 1.96*se
        lo, hi = np.percentile(bs, [2.5, 97.5])
        cp += lo <= 1.0 <= hi
        cb += 2*th - hi <= 1.0 <= 2*th - lo
    print(n, round(cn/M, 3), round(cp/M, 3), round(cb/M, 3))

출력:

25 0.743 0.744 0.726
100 0.845 0.851 0.831
400 0.914 0.916 0.907

세 방법 모두 명목수준에 크게 못 미친다. \(n = 25\)에서 \(0.74\), \(n = 400\)에서도 \(0.92\)에 그친다.

포함확률 오차가 \(n^{-1/2}\)에 비례하여 줄어드는 것을 확인할 수 있다. \(0.95 - 0.744 = 0.206\), \(0.95 - 0.851 = 0.099\), \(0.95 - 0.916 = 0.034\)이고, \(n\)이 16배가 될 때 오차가 대략 4배 줄었다(\(0.206/0.099 = 2.08 \approx \sqrt{4}\), \(0.099/0.034 = 2.9 \approx \sqrt{4}\)). 이것이 "1차 정확도, \(O(n^{-1/2})\)"의 구체적 의미이다.

왜 분산이 이렇게 어려운가. \(s^2\)의 표본분포는 모집단의 4차 적률에 지배되는데, 지수분포의 초과첨도는 \(6\)이다. 붓스트랩 재표집이 원표본의 4차 적률을 그대로 물려받으므로, \(n\)이 작으면 그 추정 자체가 매우 불안정하다.

어떻게 개선하는가. BCa나 붓스트랩-\(t\)를 쓰면 2차 정확도 \(O(n^{-1})\)를 얻는다. 또는 로그변환처럼 분산을 안정화하는 변환 위에서 구간을 만든 뒤 되돌리는 방법도 있다.

자료를 고정하고 붓스트랩 절차만 반복한다. 이때 관측되는 변동은 순전히 몬테카를로 오차이다.

import numpy as np
rng = np.random.default_rng(5)
n = 50
x = rng.exponential(1, n)          # 자료 고정

for B in (500, 2000, 10000, 50000):
    los, his = [], []
    for _ in range(200):
        idx = rng.integers(0, n, (B, n))
        bs = x[idx].mean(axis=1)
        lo, hi = np.percentile(bs, [2.5, 97.5])
        los.append(lo); his.append(hi)
    print(B, round(np.std(los), 5), round(np.std(his), 5),
          round(np.mean(his) - np.mean(los), 5))

출력:

500 0.01 0.01093 0.35653
2000 0.00514 0.00597 0.35877
10000 0.00197 0.00235 0.35929
50000 0.0011 0.00112 0.35964

끝점의 몬테카를로 표준편차가 \(\sqrt{B}\)에 반비례하여 줄어든다. \(B\)가 100배가 되면 변동이 약 10배 줄어든다(\(0.0100 \to 0.0011\)).

실무적 판단 기준. 구간 폭이 \(0.359\)인데 \(B = 500\)에서 끝점의 변동이 \(0.010\)이다. 폭의 2.8%로, 무시할 만하다고 보기에는 애매하다. \(B = 10{,}000\)이면 \(0.0020\)으로 폭의 0.55%가 되어 안전하다.

경험칙: 끝점의 몬테카를로 표준오차가 구간 폭의 1% 아래가 되도록 \(B\)를 정한다. 대체로 \(B = 5{,}000\)--\(10{,}000\)이면 충분하다.

한 가지 함정. \(B\)를 늘려도 자료 자체의 변동은 전혀 줄지 않는다. 위 표의 평균 구간 폭이 \(B\)에 따라 거의 변하지 않는 것에 주목하라. \(B\)는 붓스트랩 근사의 잡음만 줄일 뿐, 표본이 하나뿐이라는 근본적 한계는 그대로이다.

\(\hat{\theta} = \bar{X}^2\)은

이므로 \(1/n\)만큼 위로 편향되어 있다.

import numpy as np
rng = np.random.default_rng(8)
n, B, M = 20, 2000, 3000
mu, true = 0.5, 0.25

cov_p = cov_b = 0
for _ in range(M):
    x = rng.normal(mu, 1, n)
    th = x.mean() ** 2
    idx = rng.integers(0, n, (B, n))
    bs = x[idx].mean(axis=1) ** 2
    lo, hi = np.percentile(bs, [2.5, 97.5])
    cov_p += lo <= true <= hi
    cov_b += 2*th - hi <= true <= 2*th - lo
print(round(cov_p/M, 3), round(cov_b/M, 3))

출력:

0.932 0.732

\(\mu = 0.5\), \(n = 20\), 참값 \(\theta = 0.25\):

기본 구간이 훨씬 나쁘다. 편향을 보정하도록 설계된 방법이 오히려 무너진다.

왜 그런가. 두 힘이 반대로 작용한다.

1. 편향. 붓스트랩 복제값 \(\bar{X}^{*2}\)의 평균이 \(\hat{\theta}\)보다 크다. 실제로 확인하면

   x = rng.normal(0.5, 1, 20)
   idx = rng.integers(0, 20, (20000, 20))
   bs = x[idx].mean(axis=1) ** 2
   print(round(x.mean()**2, 4), round(bs.mean(), 4))   # 0.7267  0.8094
   

   출력:

   0.7267 0.8094
   

   차이 \(0.083\)은 정확히 \(\hat\sigma^2/n\)이다(\(\bar{X}^{*2}\)의 기댓값이 \(\hat\theta + \text{Var}^*(\bar{X}^*)\)이므로). 이 편향만 보면 기본 구간의 반사가 옳은 방향이다.

2. 치우침. \(\bar{X}^2\)은 제곱이므로 \(0\)이라는 바닥이 있고 붓스트랩 분포가 오른쪽으로 심하게 치우친다. 참 표본분포도 같은 방향으로 치우쳐 있다. 기본 구간의 반사가 이 치우침을 왼쪽으로 뒤집어 버린다.

이 자료에서는 두 번째 효과가 압도적이라 기본 구간이 크게 실패한다. 붓스트랩 방법 연습문제 1에서 지수분포 중앙값에 대해 본 것과 정확히 같은 현상이다.

그렇다면 편향은 어떻게 다루는가. 반사가 아니라 명시적인 편향보정을 쓴다.

• 복제값에서 \(\widehat{\text{Bias}} = \bar{\hat\theta}^* - \hat\theta\)를 빼고 백분위수 구간을 만든다.
• 또는 BCa를 쓴다. \(\hat z_0\)이 편향을, \(\hat a\)가 치우침을 각각 따로 조정하므로 두 효과가 서로 상쇄되지 않는다.

경계 모수의 함정. \(\mu = 0\)이면 상황이 더 나빠진다. \(\theta = 0\)이 모수공간의 경계에 있어 \(\hat\theta \ge 0\)이 언제나 성립하고, 어떤 붓스트랩 구간도 \(0\)을 잘 덮지 못한다. 경계 모수는 붓스트랩의 알려진 실패 지점이며 BCa로도 고쳐지지 않는다.