붓스트랩 재표집 방법¶

import numpy as np
import pandas as pd

# 난수 씨앗 고정
np.random.seed(seed=1)

# 소득 자료를 흉내 낸다. 중앙값의 강건함을 보이기에 알맞은 모양이다
loans_income = pd.Series(np.random.exponential(scale=50000, size=5000) + 20000)

# 원래 표본의 중앙값
original_median = loans_income.median()
print(f"Original sample median: ${original_median:,.0f}")
# Original sample median: $54,971

# 붓스트랩: 1000번 재표집한다
bootstrap_medians = []
for nrepeat in range(1000):
    # 원래 표본과 같은 크기로 복원추출한다
    bootstrap_sample = loans_income.sample(frac=1, replace=True)
    bootstrap_medians.append(bootstrap_sample.median())

bootstrap_medians = pd.Series(bootstrap_medians)

# 재표본마다 통계량을 구한다
bootstrap_mean = bootstrap_medians.mean()
bootstrap_std = bootstrap_medians.std()
bias = bootstrap_mean - original_median

print(f"  Mean of bootstrap distribution:  ${bootstrap_mean:,.0f}")   # $55,008
print(f"  Standard error of median:        ${bootstrap_std:,.0f}")    # $756
print(f"  Bias of median estimator:        ${bias:,.0f}")             # $37

Original sample median: $54,971
  Mean of bootstrap distribution:  $55,008
  Standard error of median:        $756
  Bias of median estimator:        $37

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 5))
ax.hist(bootstrap_medians, bins=40, color='steelblue', edgecolor='black', alpha=0.7)
ax.axvline(original_median, color='darkred', linestyle='--', linewidth=2.5,
           label=f'Original median: ${original_median:,.0f}')
ax.axvline(bootstrap_mean, color='green', linestyle='--', linewidth=2.5,
           label=f'Bootstrap mean: ${bootstrap_mean:,.0f}')
ax.set_xlabel('Median Income ($)', fontsize=11)
ax.set_ylabel('Frequency', fontsize=11)
ax.set_title('Bootstrap Distribution of the Sample Median', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.spines[['top', 'right']].set_visible(False)
ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()

theta, shift, n = 50000.0, 20000.0, len(loans_income)
m_pop = shift + theta * np.log(2)
f_m = 1 / (2 * theta)
print(f"모집단 중앙값 = {m_pop:,.2f},  f(m) = {f_m:.1e}")
print(f"중앙값의 점근 SE = {1 / (2 * np.sqrt(n) * f_m):,.3f}")
print(f"평균의   SE      = {theta / np.sqrt(n):,.3f}")
print(f"붓스트랩이 준 값 = {bootstrap_std:,.3f}"
      f"   (B=1000 의 몬테카를로 요동 {bootstrap_std / np.sqrt(2 * 1000):.1f})")

# 표본을 바꾸어 가며 붓스트랩 표준오차가 얼마나 흔들리는지 잰다.
rng = np.random.default_rng(7)
R, Bk = 20, 500
se_med, se_mean = [], []
for _ in range(R):
    x = rng.exponential(theta, n) + shift
    rs = x[rng.integers(0, n, (Bk, n))]
    se_med.append(np.median(rs, axis=1).std(ddof=1))
    se_mean.append(rs.mean(axis=1).std(ddof=1))
se_med, se_mean = np.array(se_med), np.array(se_mean)
print(f"\n표본 {R} 개, 각 B = {Bk}")
print(f"  중앙값 SE: 평균 {se_med.mean():.1f},  SD {se_med.std(ddof=1):.1f}"
      f"  (상대 {se_med.std(ddof=1) / se_med.mean():.3f}),"
      f"  범위 {se_med.min():.1f}~{se_med.max():.1f}")
print(f"  평균   SE: 평균 {se_mean.mean():.1f},  SD {se_mean.std(ddof=1):.1f}"
      f"  (상대 {se_mean.std(ddof=1) / se_mean.mean():.3f}),"
      f"  범위 {se_mean.min():.1f}~{se_mean.max():.1f}")
print(f"  몬테카를로 몫만이면 SD = {707.107 / np.sqrt(2 * Bk):.1f}")

모집단 중앙값 = 54,657.36,  f(m) = 1.0e-05
중앙값의 점근 SE = 707.107
평균의   SE      = 707.107
붓스트랩이 준 값 = 755.791   (B=1000 의 몬테카를로 요동 16.9)

표본 20 개, 각 B = 500
  중앙값 SE: 평균 684.9,  SD 60.4  (상대 0.088),  범위 598.8~857.9
  평균   SE: 평균 703.0,  SD 18.9  (상대 0.027),  범위 674.4~742.5
  몬테카를로 몫만이면 SD = 22.4

# 붓스트랩 반복 횟수 B 를 늘리면 추정이 안정된다. 다만 B 는 붓스트랩
# 자체의 몬테카를로 오차만 줄일 뿐, 표본크기 n 이 주는 한계를 넘지는 못한다.
# 표준오차 추정에는 수백 번이면 족하고, 신뢰구간에는 수천 번이 필요하다.
np.random.seed(1)
income = loans_income.values

for B in [100, 500, 1000, 5000]:
    boot_medians = np.array([
        np.median(np.random.choice(income, size=len(income), replace=True))
        for _ in range(B)
    ])
    print(f"B = {B:5d}: SE = ${boot_medians.std():8,.1f}")

B =   100: SE = $   696.8
B =   500: SE = $   760.6
B =  1000: SE = $   741.7
B =  5000: SE = $   759.5

rng = np.random.default_rng(3)
n = len(income)
print(f"{'B':>6}{'반복':>6}{'SE 평균':>12}{'SE 의 SD':>12}{'예측 SE/sqrt(2B)':>18}")
for B, R in [(100, 100), (500, 40), (1000, 30)]:
    ses = []
    for _ in range(R):
        rs = income[rng.integers(0, n, (B, n))]
        ses.append(np.median(rs, axis=1).std(ddof=0))
    ses = np.array(ses)
    print(f"{B:>6}{R:>6}{ses.mean():>12.1f}{ses.std(ddof=1):>12.1f}"
          f"{ses.mean() / np.sqrt(2 * B):>18.1f}")

     B    반복       SE 평균     SE 의 SD    예측 SE/sqrt(2B)
   100   100       746.5        57.6              52.8
   500    40       756.2        20.5              23.9
  1000    30       756.8        17.5              16.9

def block_bootstrap(data, block_size, n_bootstrap, rng=None):
    """이동블록 붓스트랩. 서로 독립이 아닌 자료에 쓴다.

    보통의 붓스트랩은 관측값을 하나씩 뽑으므로 시간 순서의 구조가 모두
    부서진다. 시계열처럼 이웃한 값끼리 얽혀 있는 자료에서는 표준오차를
    크게 낮잡게 된다. 그래서 낱값이 아니라 길이 block_size 짜리 덩어리를
    통째로 뽑아, 덩어리 안의 상관은 그대로 남긴다.
    """
    rng = np.random.default_rng() if rng is None else rng
    n = len(data)
    n_blocks = int(np.ceil(n / block_size))
    samples = []
    for _ in range(n_bootstrap):
        starts = rng.integers(0, n - block_size + 1, n_blocks)
        sample = np.concatenate([data[s:s + block_size] for s in starts])[:n]
        samples.append(sample)
    return np.array(samples)

rho, n, B = 0.7, 500, 1000
rng = np.random.default_rng(0)

# AR(1) 을 주변분산 1 로 만든다.
x = np.empty(n)
x[0] = rng.normal(0, 1)
eps = rng.normal(0, np.sqrt(1 - rho ** 2), n)
for t in range(1, n):
    x[t] = rho * x[t - 1] + eps[t]

k = np.arange(1, n)
var_true = (n + 2 * np.sum((n - k) * rho ** k)) / n ** 2
print(f"표본 sd = {x.std(ddof=1):.4f},  표본 1차 자기상관 = "
      f"{np.corrcoef(x[:-1], x[1:])[0, 1]:.4f}")
print(f"참 SD(xbar)      = {np.sqrt(var_true):.5f}   (점근 "
      f"{np.sqrt((1 + rho) / (1 - rho) / n):.5f})")
print(f"독립이라면       = {np.sqrt(1 / n):.5f}   -> 부풀림 "
      f"{np.sqrt(var_true * n):.4f} 배")

naive = x[rng.integers(0, n, (B, n))].mean(axis=1)
print(f"순진한 붓스트랩 SE = {naive.std(ddof=1):.5f}")

print(f"\n{'블록 길이':>9}{'블록붓스트랩 SE':>17}{'예측 sqrt(c/n)*s':>19}{'참값 대비':>11}")
for L in (1, 5, 10, 25, 50, 100):
    kk = np.arange(1, L)
    c = 1 + 2 * np.sum((1 - kk / L) * rho ** kk) if L > 1 else 1.0
    se = block_bootstrap(x, L, B, rng=np.random.default_rng(1)).mean(axis=1).std(ddof=1)
    print(f"{L:>9}{se:>17.5f}{np.sqrt(c / n) * x.std(ddof=0):>19.5f}"
          f"{se / np.sqrt(var_true):>11.3f}")

표본 sd = 0.9872,  표본 1차 자기상관 = 0.6782
참 SD(xbar)      = 0.10617   (점근 0.10646)
독립이라면       = 0.04472   -> 부풀림 2.3739 배
순진한 붓스트랩 SE = 0.04391

    블록 길이        블록붓스트랩 SE     예측 sqrt(c/n)*s      참값 대비
        1          0.04396            0.04411      0.414
        5          0.07466            0.07739      0.703
       10          0.08745            0.08990      0.824
       25          0.09242            0.09906      0.871
       50          0.09254            0.10207      0.872
      100          0.08801            0.10354      0.829

import numpy as np

rng = np.random.default_rng(0)
# 치우친 자료. 백분위수법이 어긋나기 쉬운 경우라 백분위수-t 를 쓴다.
data = rng.exponential(2.0, 40)
n, B = len(data), 2000

original_statistic = data.mean()
original_se = data.std(ddof=1) / np.sqrt(n)

# 각 붓스트랩 표본에서 통계량과 그 표준오차를 함께 계산한다
idx = rng.integers(0, n, (B, n))
resamples = data[idx]
bootstrap_statistics = resamples.mean(axis=1)
bootstrap_ses = resamples.std(axis=1, ddof=1) / np.sqrt(n)

# 백분위수-t 는 통계량 자체가 아니라 스튜던트화한 값의 분포를 쓴다.
# 이러면 치우침이 보정되어 포함확률이 명목수준에 더 가까워진다.
# 대신 붓스트랩 표본마다 표준오차를 또 계산해야 해 비용이 든다.
bootstrap_t_stats = (bootstrap_statistics - original_statistic) / bootstrap_ses
# 위아래가 뒤집혀 들어간다. t 분포의 위쪽 분위점이 구간의 아래끝을 만든다.
ci_lower = original_statistic - np.percentile(bootstrap_t_stats, 97.5) * original_se
ci_upper = original_statistic - np.percentile(bootstrap_t_stats, 2.5) * original_se

print(f"theta_hat = {original_statistic:.4f}")
print(f"bootstrap-t CI = ({ci_lower:.4f}, {ci_upper:.4f})")
print("percentile  CI = ({:.4f}, {:.4f})".format(
    *np.percentile(bootstrap_statistics, [2.5, 97.5])))

theta_hat = 2.3157
bootstrap-t CI = (1.6976, 3.3636)
percentile  CI = (1.6498, 3.0813)

from scipy import stats

# 지수분포에서는 2n*xbar/theta ~ chi^2_{2n} 이라 정확 구간이 있다.
lo_e = 2 * n * original_statistic / stats.chi2.ppf(0.975, 2 * n)
up_e = 2 * n * original_statistic / stats.chi2.ppf(0.025, 2 * n)
pc = np.percentile(bootstrap_statistics, [2.5, 97.5])
print(f"정확 (카이제곱) CI = ({lo_e:.4f}, {up_e:.4f})")
print(f"표본평균의 왜도 (이론 2/sqrt(n)) = {2 / np.sqrt(n):.4f}")
for name, (a, b) in [("백분위수", pc), ("백분위수-t", (ci_lower, ci_upper)),
                     ("정확", (lo_e, up_e))]:
    print(f"  {name:>10}: 왼쪽 {original_statistic - a:.4f}"
          f"  오른쪽 {b - original_statistic:.4f}"
          f"  비 {(b - original_statistic) / (original_statistic - a):.3f}"
          f"  폭 {b - a:.4f}")

# 포함확률 모의실험
M, Bc, theta = 1000, 1000, 2.0
rng2 = np.random.default_rng(123)
cp = ct = ce = 0
lp = rp = lt = rt = 0
for _ in range(M):
    d = rng2.exponential(theta, n)
    th = d.mean(); se = d.std(ddof=1) / np.sqrt(n)
    rs = d[rng2.integers(0, n, (Bc, n))]
    m = rs.mean(axis=1); s = rs.std(axis=1, ddof=1) / np.sqrt(n)
    a, b = np.percentile(m, [2.5, 97.5])
    cp += (a <= theta <= b); lp += theta < a; rp += theta > b
    tt = (m - th) / s
    lo = th - np.percentile(tt, 97.5) * se
    up = th - np.percentile(tt, 2.5) * se
    ct += (lo <= theta <= up); lt += theta < lo; rt += theta > up
    le = 2 * n * th / stats.chi2.ppf(0.975, 2 * n)
    ue = 2 * n * th / stats.chi2.ppf(0.025, 2 * n)
    ce += (le <= theta <= ue)
print(f"\nM = {M}, B = {Bc}, n = {n}, theta = {theta}"
      f"   (포함률의 몬테카를로 오차 {np.sqrt(0.95 * 0.05 / M):.4f})")
print(f"백분위수   포함률 {cp / M:.3f}  (왼쪽 밖 {lp / M:.3f}, 오른쪽 밖 {rp / M:.3f})")
print(f"백분위수-t 포함률 {ct / M:.3f}  (왼쪽 밖 {lt / M:.3f}, 오른쪽 밖 {rt / M:.3f})")
print(f"정확       포함률 {ce / M:.3f}")

정확 (카이제곱) CI = (1.7374, 3.2414)
표본평균의 왜도 (이론 2/sqrt(n)) = 0.3162
        백분위수: 왼쪽 0.6659  오른쪽 0.7656  비 1.150  폭 1.4315
      백분위수-t: 왼쪽 0.6180  오른쪽 1.0480  비 1.696  폭 1.6660
          정확: 왼쪽 0.5783  오른쪽 0.9257  비 1.601  폭 1.5040

M = 1000, B = 1000, n = 40, theta = 2.0   (포함률의 몬테카를로 오차 0.0069)
백분위수   포함률 0.916  (왼쪽 밖 0.021, 오른쪽 밖 0.063)
백분위수-t 포함률 0.934  (왼쪽 밖 0.030, 오른쪽 밖 0.036)
정확       포함률 0.948

import numpy as np, pandas as pd
np.random.seed(1)
x = pd.Series(np.random.exponential(scale=50000, size=5000) + 20000)
rng = np.random.default_rng(0)

idx = rng.integers(0, 5000, (5000, 5000))
boot_mean = x.values[idx].mean(axis=1)
boot_med  = np.median(x.values[idx], axis=1)

print("평균  붓스트랩 SE:", round(boot_mean.std(ddof=1), 1))
print("이론값 s/sqrt(n)  :", round(x.std(ddof=1) / np.sqrt(5000), 1))
print("중앙값 붓스트랩 SE:", round(boot_med.std(ddof=1), 1))

평균  붓스트랩 SE: 681.9
이론값 s/sqrt(n)  : 690.9
중앙값 붓스트랩 SE: 755.0

import numpy as np
rng = np.random.default_rng(0)
n, phi, M = 500, 0.8, 2000

def ar1(n):
    e = rng.normal(0, 1, n + 200)
    x = np.zeros(n + 200)
    for t in range(1, n + 200):
        x[t] = phi * x[t-1] + e[t]
    return x[200:]

# 참 SE: 시계열을 여러 번 생성
true_se = np.array([ar1(n).mean() for _ in range(M)]).std(ddof=1)

x = ar1(n)
B = 4000
iid_se = x[rng.integers(0, n, (B, n))].mean(axis=1).std(ddof=1)

def block_se(x, L, B):
    n = len(x); nb = int(np.ceil(n / L))
    starts = rng.integers(0, n - L + 1, (B, nb))
    means = np.array([np.concatenate([x[s:s+L] for s in row])[:n].mean()
                      for row in starts])
    return means.std(ddof=1)

print("참 SE       :", round(true_se, 4))
print("표준 붓스트랩:", round(iid_se, 4))
for L in (5, 20, 50):
    print(f"블록 L={L:2d}  :", round(block_se(x, L, B), 4))

참 SE       : 0.2229
표준 붓스트랩: 0.0742
블록 L= 5  : 0.1421
블록 L=20  : 0.1996
블록 L=50  : 0.1882

import numpy as np
from scipy import stats
rng = np.random.default_rng(4)
n, B, M = 200, 800, 600

for p, target in [(50, 0.0), (90, stats.norm.ppf(0.90)),
                  (99, stats.norm.ppf(0.99))]:
    cov = 0
    for _ in range(M):
        x = rng.normal(0, 1, n)
        q = np.percentile(x[rng.integers(0, n, (B, n))], p, axis=1)
        lo, hi = np.percentile(q, [2.5, 97.5])
        cov += lo <= target <= hi
    print(p, round(cov / M, 3))

50 0.942
90 0.942
99 0.857

x = rng.normal(0, 1, 200)
q = np.percentile(x[rng.integers(0, 200, (5000, 200))], 99, axis=1)
print(len(np.unique(np.round(q, 6))))   # 60  ← 5000개 복제값이 60가지 값만 갖는다

60