대안으로서의 붓스트랩¶

표본평균은 \(\bar{x} = 6.37\%\), 표본표준편차는 \(s = 9.20\), 표본왜도는 \(0.257\)이다. 모집단 분포를 가정하지 않고 모평균의 95% 신뢰구간을 만들기 위해 붓스트랩을 적용한다.

1단계. 크기 12인 붓스트랩 표본을 복원추출로 \(B = 5000\)개 뽑는다.

2단계. 각 붓스트랩 표본의 평균을 계산하여 붓스트랩 분포 \(\{\hat{\theta}^{*1}, \ldots, \hat{\theta}^{*5000}\}\)을 만든다.

3단계. 붓스트랩 평균들을 정렬하여 2.5번째와 97.5번째 백분위수를 뽑는다.

아래 코드(난수 씨앗 42)를 실행하면 \(\hat{\theta}^*_{(0.025)} = 1.48\), \(\hat{\theta}^*_{(0.975)} = 11.33\)이 나온다. 모평균 수익률의 95% 백분위수 붓스트랩 신뢰구간은 \([1.48\%,\; 11.33\%]\)이다.

비교를 위해 표준 \(t\) 구간은

곧 \([0.52\%,\; 12.21\%]\)이다.

(1) 해석적으로. 자료 \(\mathbf{x} = (x_1, \ldots, x_n)\)을 고정하고 생각한다. 붓스트랩 추출은 경험분포 \(\hat F_n\)에서의 독립동일분포 추출이고, \(\hat F_n\)은 각 \(x_i\)에 확률 \(1/n\)을 준다. 그러므로 붓스트랩 한 점 \(X^*\)의 평균과 분산은

이다. 여기서 \(s_{\hat F}^2 = \frac{1}{n}\sum (x_i - \bar x)^2\)으로, \(n-1\)이 아니라 \(n\)으로 나눈 것이다. \(X_1^*, \ldots, X_n^*\)이 독립이므로

곧 평균의 붓스트랩 표준오차는 \(s_{\hat F}/\sqrt{n}\)이다. 근사가 아니라 정확한 등식이다. \(B\)를 쓰는 모의실험은 이 값을 추정하는 수단일 뿐이고, \(B \to \infty\)에서 바로 이 값으로 간다.

이 표본에서는 \(s_{\hat F} = 8.8095\), \(\sqrt{n} = 3.4641\)이므로

이다. 그러므로 답은 \(s_{\hat F}/\sqrt{n} = 2.54\) 쪽이다. \(s/\sqrt{n} = 9.2012/3.4641 = 2.6562\)가 아니다. 두 값의 비는 \(\sqrt{(n-1)/n} = \sqrt{11/12} = 0.9574\)로, \(n\)이 작을수록 벌어진다.

(2) 수치적으로. 먼저 쪽에 실린 코드를 그대로 돌린다.

"""표본평균의 붓스트랩 신뢰구간.

치우친 자료에 비모수 붓스트랩을 적용하고, 정규성을 가정한 t 구간과
견준다. 붓스트랩은 분포를 가정하지 않고 자료 자체를 모집단으로 삼아
재표본을 뽑으므로, 정규성이 의심스러울 때 쓸 만한 대안이 된다.
"""

import numpy as np

# ===================================================================
# 자료 — 수익률 12개. 표본이 작고 치우쳐 있어 정규성을 믿기 어렵다.
# ===================================================================
data = np.array([-8.2, 3.1, 15.7, 2.4, -1.5, 22.3,
                  6.8, -3.4, 11.2, 1.9, 18.6, 7.5])

# ===================================================================
# 붓스트랩 — 자료에서 복원추출로 같은 크기의 재표본을 5000번 만들고
# 그때마다 평균을 기록한다. 그 5000개의 흩어짐이 곧 표본평균의
# 표집분포에 대한 추정이다.
# ===================================================================
rng = np.random.default_rng(42)
B = 5000
n = len(data)

bootstrap_means = np.array([
    np.mean(rng.choice(data, size=n, replace=True))
    for _ in range(B)
])

# ===================================================================
# 백분위 신뢰구간 — 재표본 평균들의 2.5·97.5 백분위점을 그대로 쓴다.
# t 구간과 달리 좌우가 대칭일 까닭이 없어, 치우친 자료에서는 한쪽이
# 더 길게 나온다. 그것이 오히려 자료를 제대로 담은 결과다.
# ===================================================================
ci_lower, ci_upper = np.percentile(bootstrap_means, [2.5, 97.5])

if __name__ == "__main__":
    print(f"Sample mean: {np.mean(data):.2f}%")
    print(f"Bootstrap SE: {np.std(bootstrap_means, ddof=0):.2f}%")
    print(f"95% Percentile CI: [{ci_lower:.2f}%, {ci_upper:.2f}%]")

출력:

Sample mean: 6.37%
Bootstrap SE: 2.50%
95% Percentile CI: [1.48%, 11.33%]

이제 (1)의 이론값과 맞춰 보고, 구간 폭을 분해한다.

import numpy as np
from scipy import stats

data = np.array([-8.2, 3.1, 15.7, 2.4, -1.5, 22.3,
                  6.8, -3.4, 11.2, 1.9, 18.6, 7.5])
n = len(data)
s_hat = data.std(ddof=0)      # 1/n 로 나눈 것. 경험분포의 표준편차다.
s = data.std(ddof=1)          # 1/(n-1) 로 나눈 보통의 표본표준편차

print(f"이론  s_hat/sqrt(n) = {s_hat:.4f} / {np.sqrt(n):.4f} = {s_hat/np.sqrt(n):.4f}")
print(f"대조  s/sqrt(n)     = {s:.4f} / {np.sqrt(n):.4f} = {s/np.sqrt(n):.4f}")

# B 를 키우면 몬테카를로 추정이 이론값으로 간다.
for B in [5_000, 50_000, 500_000]:
    rng = np.random.default_rng(42)
    idx = rng.integers(0, n, size=(B, n))
    boot = data[idx].mean(axis=1)
    print(f"B = {B:7d}  붓스트랩 SE = {boot.std(ddof=0):.4f}")

# 구간 폭을 세 단계로 분해한다.
t = stats.t.ppf(0.975, n - 1)
z = stats.norm.ppf(0.975)
print(f"t 구간 폭            = 2 x {t:.4f} x {s/np.sqrt(n):.4f} = {2*t*s/np.sqrt(n):.3f}")
print(f"t 를 z 로 바꾸면     = 2 x {z:.4f} x {s/np.sqrt(n):.4f} = {2*z*s/np.sqrt(n):.3f}")
print(f"s 를 s_hat 로 바꾸면 = 2 x {z:.4f} x {s_hat/np.sqrt(n):.4f} = {2*z*s_hat/np.sqrt(n):.3f}")

출력:

이론  s_hat/sqrt(n) = 8.8095 / 3.4641 = 2.5431
대조  s/sqrt(n)     = 9.2012 / 3.4641 = 2.6562
B =    5000  붓스트랩 SE = 2.4994
B =   50000  붓스트랩 SE = 2.5431
B =  500000  붓스트랩 SE = 2.5396
t 구간 폭            = 2 x 2.2010 x 2.6562 = 11.692
t 를 z 로 바꾸면     = 2 x 1.9600 x 2.6562 = 10.412
s 를 s_hat 로 바꾸면 = 2 x 1.9600 x 2.5431 = 9.969

(1)의 확인. 유도한 값은 \(2.5431\)이고 \(B = 5000\)의 모의실험은 \(2.4994\)를 주었다. 차이가 \(0.044\)인데, 표준편차 추정의 몬테카를로 오차가 대략 \(\mathrm{SE}/\sqrt{2B} = 2.543/100 = 0.025\)이므로 \(1.7\)배쯤 벗어난 셈이다. 우연으로 설명되는 크기다. \(B\)를 열 배로 키우면 \(2.5431\), 백 배로 키우면 \(2.5396\)(오차 척도 \(0.0025\))이 되어 이론값 주위를 맴돈다. \(B\)를 아무리 키워도 \(2.6562\) 쪽으로는 가지 않는다. 쪽에 적힌 \(2.50\%\)는 \(s/\sqrt{n}\)의 추정값이 아니라 \(s_{\hat F}/\sqrt{n}\)의 추정값이다.

여기서 눈여겨볼 것은 \(B\)가 지우는 오차와 \(n\)이 지우는 오차가 다른 것이라는 점이다. \(B\)를 키우면 몬테카를로 잡음만 사라지고 \(2.5431\)이라는 값 자체는 꼼짝하지 않는다. 그 값은 자료 12개가 정해 놓은 것이다. 붓스트랩 복제를 백만 번 한다고 표본이 커지지는 않는다.

(2) 폭의 차이는 치우침에서 오지 않는다. 분해하면 이렇다.

세 조각의 몫이 \(1.280 : 0.443 : 0.127\)이다. 가장 큰 조각은 \(t\) 분위수를 잃은 것이다. 백분위수 붓스트랩은 \(n-1 = 11\)이라는 자유도 보정을 아예 하지 않으므로, 꼬리가 두꺼운 \(t_{11}\) 대신 정규 분위수에 해당하는 폭을 쓴다. 둘째 조각이 (1)에서 본 \(\sqrt{11/12}\) 축소다. 남은 \(0.127\)은 붓스트랩 분포가 정규가 아닌 데서 오는 몫과 \(B = 5000\)의 몬테카를로 잡음이 섞인 것이다.

그러므로 "치우침을 반영했기 때문"이라는 설명은 옳지 않다. 실제로 치우침은 거의 없다. 이 표본의 왜도는 \(g_1 = 0.2569\)(scipy.stats.skew의 기본값 bias=True)로 가볍고, 붓스트랩 평균 5000개의 왜도는 \(g_1 = 0.0719\)로 사실상 대칭이다. 구간의 좌우도 중심 \(\bar x = 6.367\)에서 왼쪽 \(4.883\), 오른쪽 \(4.958\)로 \(1.5\%\)밖에 차이 나지 않는다. 치우침이 폭을 \(1.85\)나 줄일 수는 없다.

쪽 앞의 경고 상자가 말하는 그대로다. 백분위수 붓스트랩이 작은 표본에서 좁아지는 것은 자료의 모양을 잘 담아서가 아니라 자유도 보정과 분산 보정을 둘 다 빼먹어서다. \(n\)이 작을 때 BCa나 붓스트랩-\(t\)를 권하는 이유가 여기 있다.

1. 크기 \(n\)인 원래 표본에서 (각각 크기 \(n\)인) 붓스트랩 표본 \(B\)개를 복원추출로 뽑는다.
2. 각 붓스트랩 표본 \(b = 1, \dots, B\)에 대해 \(\bar{X}^*_b\)(붓스트랩 표본평균)을 계산한다.
3. \(\bar{X}^*\)의 붓스트랩 분포에서 신뢰구간을 만든다.
4. 백분위수 방법: \(B\)개 붓스트랩 평균의 \(\alpha/2\)와 \(1-\alpha/2\) 분위수를 쓴다.
5. 붓스트랩-t 방법: 붓스트랩 t 통계량을 계산하고 그 분위수를 쓴다.
6. 흔한 선택은 \(B = 10{,}000\)번의 붓스트랩 복제이다.

정규성 가정이 필요 없다. \(\bar{X}^*\)의 붓스트랩 분포가 \(\bar{X}\)의 표집분포를 흉내 내므로 치우쳤거나 꼬리가 두꺼운 자료에서도 타당한 추론이 가능하다.

1. 작은 \(n\)에서의 극단 분위수: 붓스트랩은 원래 표본의 범위를 벗어나는 값을 만들어 낼 수 없다. 작은 \(n\)에서 극단 분위수(예: 99번째 백분위수)를 추정하면 꼬리의 변동을 과소평가한다.

2. 종속 자료: 표준 비모수 붓스트랩은 관측값이 독립동일분포라고 가정한다. 자료가 종속이면(시계열, 공간 자료) 개별 관측값을 재표본추출하는 것이 종속 구조를 파괴하여 타당하지 않은 결과를 낸다. 블록 붓스트랩이나 정상 붓스트랩이 필요하다.

추가로: (3) \(n\)이 매우 작을 때(경험분포가 모집단을 잘 근사하지 못한다), (4) 불규칙한 통계량(예: 표본최댓값, 붓스트랩 분포가 제대로 수렴하지 않는다).

중간 정도의 \(n\)을 갖는 오른쪽 치우침 자료에서는 붓스트랩 백분위수 구간의 포함확률이 더 나을 것으로 기대된다. 표집분포의 대칭성을 가정하지 않기 때문이다.

\(t\) 구간 \(\bar{X} \pm t_{\alpha/2} s/\sqrt{n}\)은 \(\bar{X}\)를 중심으로 대칭이지만, \(\bar{X}\)의 참 표집분포는 (자료에서 물려받아) 치우쳐 있다. 그래서 \(t\) 구간은 한쪽을 과다 포함하고 다른 쪽을 과소 포함한다.

붓스트랩 백분위수 구간은 이 비대칭을 자연스럽게 포착한다. \(\bar{X}^*\)의 붓스트랩 분포가 치우쳐 있으면 구간도 비대칭이 되어 참 표집분포에 더 잘 맞는다. 치우친 자료에서 최선의 성능을 원한다면 편향보정 가속(BCa) 붓스트랩 구간이 권장된다.

import numpy as np

rng = np.random.default_rng(42)
data = rng.exponential(scale=5, size=50)  # skewed data

B = 10000
boot_medians = np.array([
    np.median(rng.choice(data, size=len(data), replace=True))
    for _ in range(B)
])

ci_lower = np.percentile(boot_medians, 2.5)
ci_upper = np.percentile(boot_medians, 97.5)
print(f"Sample median: {np.median(data):.3f}")
print(f"95% Bootstrap CI: ({ci_lower:.3f}, {ci_upper:.3f})")

출력:

Sample median: 3.243
95% Bootstrap CI: (1.831, 5.420)

중앙값의 표집분포에는 간단한 닫힌 형태의 공식이 없지만, 붓스트랩은 평균일 때와 똑같은 절차로 구간을 준다. 백분위수 방법을 쓴 것이다. 이 구간은 정규성 가정 없이 타당하며, 반복문 안에서 계산하는 통계량만 바꾸면 어떤 통계량(평균, 중앙값, 상관 등)에도 쓸 수 있다.