붓스트랩 신뢰구간: 신뢰수준의 시각적 해석¶

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.utils import resample

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

# 소득 자료를 흉내 낸다. 실제 대출 소득 자료를 읽어 써도 된다
loans_income = np.random.exponential(scale=50000, size=5000) + 20000

# 모집단에서 크기 20 짜리 표본 하나를 뽑는다
original_sample = resample(loans_income, n_samples=20, replace=False)
original_mean = original_sample.mean()

print(f"Original sample size: {len(original_sample)}")   # 20
print(f"Original sample mean: ${original_mean:,.0f}")    # $67,846

# 붓스트랩: 그 표본에서 500번 재표집한다
bootstrap_means = []
for _ in range(500):
    bootstrap_sample = resample(original_sample)  # with replacement
    bootstrap_means.append(bootstrap_sample.mean())

bootstrap_means = pd.Series(bootstrap_means)

# 신뢰구간 계산
ci_90_lower, ci_90_upper = bootstrap_means.quantile([0.05, 0.95])
ci_95_lower, ci_95_upper = bootstrap_means.quantile([0.025, 0.975])

print("90% CI: [${:,.0f}, ${:,.0f}]".format(ci_90_lower, ci_90_upper))
# 90% CI: [$49,742, $90,491]
print("95% CI: [${:,.0f}, ${:,.0f}]".format(ci_95_lower, ci_95_upper))
# 95% CI: [$47,026, $95,545]
print(f"Mean of bootstrap means: ${bootstrap_means.mean():,.0f}")   # $68,444

Original sample size: 20
Original sample mean: $67,846
90% CI: [$49,742, $90,491]
95% CI: [$47,026, $95,545]
Mean of bootstrap means: $68,444

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 첫째 그림: 90% 신뢰구간
ax1.hist(bootstrap_means, bins=30, color='steelblue', edgecolor='black', alpha=0.7)
ax1.axvline(ci_90_lower, color='darkred', linestyle='--', linewidth=2.5, label='90% CI limits')
ax1.axvline(ci_90_upper, color='darkred', linestyle='--', linewidth=2.5)
ax1.axvspan(ci_90_lower, ci_90_upper, alpha=0.2, color='green', label='90% CI')

ci_90_mid = (ci_90_lower + ci_90_upper) / 2
ax1.text(ci_90_mid, 35, f'90% CI\n[${ci_90_lower:,.0f}, ${ci_90_upper:,.0f}]',
         ha='center', va='center', fontsize=10,
         bbox=dict(boxstyle='round', facecolor='white', edgecolor='darkred', linewidth=1.5))
ax1.axvline(original_mean, color='black', linestyle='-', linewidth=2,
            label=f'Sample mean: ${original_mean:,.0f}')

ax1.set_xlabel('Bootstrap Sample Mean ($)', fontsize=11)
ax1.set_ylabel('Frequency', fontsize=11)
ax1.set_title('90% Bootstrap Confidence Interval', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9)
ax1.spines[['top', 'right']].set_visible(False)
ax1.grid(True, alpha=0.3, axis='y')

# 둘째 그림: 95% 신뢰구간
ax2.hist(bootstrap_means, bins=30, color='steelblue', edgecolor='black', alpha=0.7)
ax2.axvline(ci_95_lower, color='darkblue', linestyle='--', linewidth=2.5, label='95% CI limits')
ax2.axvline(ci_95_upper, color='darkblue', linestyle='--', linewidth=2.5)
ax2.axvspan(ci_95_lower, ci_95_upper, alpha=0.2, color='orange', label='95% CI')

ci_95_mid = (ci_95_lower + ci_95_upper) / 2
ax2.text(ci_95_mid, 35, f'95% CI\n[${ci_95_lower:,.0f}, ${ci_95_upper:,.0f}]',
         ha='center', va='center', fontsize=10,
         bbox=dict(boxstyle='round', facecolor='white', edgecolor='darkblue', linewidth=1.5))
ax2.axvline(original_mean, color='black', linestyle='-', linewidth=2,
            label=f'Sample mean: ${original_mean:,.0f}')

ax2.set_xlabel('Bootstrap Sample Mean ($)', fontsize=11)
ax2.set_ylabel('Frequency', fontsize=11)
ax2.set_title('95% Bootstrap Confidence Interval', fontsize=12, fontweight='bold')
ax2.legend(loc='upper left', fontsize=9)
ax2.spines[['top', 'right']].set_visible(False)
ax2.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

n = len(original_sample)
se_lim = original_sample.std(ddof=0) / np.sqrt(n)
print(f"sigma-hat = {original_sample.std(ddof=0):,.0f},  극한 SE = {se_lim:,.0f}"
      f"   (B=500 모의 SD = {bootstrap_means.std():,.0f})")
for lab, z, lo, hi in (("90%", 1.644854, ci_90_lower, ci_90_upper),
                       ("95%", 1.959964, ci_95_lower, ci_95_upper)):
    print(f"{lab}: 폭 {hi - lo:,.0f}   예측 2 z SE = {2 * z * se_lim:,.0f}")
print(f"폭의 비 = {(ci_95_upper-ci_95_lower)/(ci_90_upper-ci_90_lower):.4f}"
      f"   (정규이론 {1.959964/1.644854:.4f})")

print(f"\n{'B':>7}{'하한 평균':>12}{'하한 SD':>10}{'상한 평균':>12}{'상한 SD':>10}{'폭의 SD':>10}")
for B in (500, 2000, 10000):
    los, his = [], []
    for s in range(30):
        r = np.random.default_rng(s)
        bm = original_sample[r.integers(0, n, (B, n))].mean(axis=1)
        lo, hi = np.percentile(bm, [2.5, 97.5])
        los.append(lo); his.append(hi)
    los, his = np.array(los), np.array(his)
    print(f"{B:>7}{los.mean():>12,.0f}{los.std(ddof=1):>10,.0f}"
          f"{his.mean():>12,.0f}{his.std(ddof=1):>10,.0f}{(his-los).std(ddof=1):>10,.0f}")

sigma-hat = 56,357,  극한 SE = 12,602   (B=500 모의 SD = 12,261)
90%: 폭 40,749   예측 2 z SE = 41,457
95%: 폭 48,519   예측 2 z SE = 49,399
폭의 비 = 1.1907   (정규이론 1.1916)

      B       하한 평균     하한 SD       상한 평균     상한 SD     폭의 SD
    500      46,568       906      94,860     1,654     1,832
   2000      46,448       539      95,074       995     1,069
  10000      46,366       192      95,099       414       459

# "95% 신뢰"는 한 번 만든 구간에 대한 확률이 아니라, 같은 절차를 되풀이할
# 때 참값을 담는 비율에 대한 약속이다. 그 비율을 직접 세어 본다.
np.random.seed(42)

true_pop = np.random.exponential(scale=50000, size=10000) + 20000
true_mean = true_pop.mean()

# 치우친 모집단에서 n=20 만 뽑으므로 실제 포함확률이 95%에 못 미친다.
# 붓스트랩이 만능이 아니라는 것을 보여 주는 대목이다.
n_simulations = 2000
ci_covers = []

for sim in range(n_simulations):
    sample = np.random.choice(true_pop, size=20, replace=False)
    boot_means = np.array([np.mean(np.random.choice(sample, size=len(sample)))
                           for _ in range(500)])
    ci_lower, ci_upper = np.percentile(boot_means, [2.5, 97.5])
    ci_covers.append(ci_lower <= true_mean <= ci_upper)

print(f"Coverage across {n_simulations} simulations: {100*np.mean(ci_covers):.1f}%")
# Coverage across 2000 simulations: 90.6%

Coverage across 2000 simulations: 90.9%

rng = np.random.default_rng(7)
for nn in (20, 50, 100, 200):
    cov, M = 0, 2000
    for _ in range(M):
        s = rng.choice(true_pop, size=nn, replace=False)
        bm = s[rng.integers(0, nn, (500, nn))].mean(axis=1)
        lo, hi = np.percentile(bm, [2.5, 97.5])
        cov += lo <= true_mean <= hi
    print(f"n={nn:>4}: 포함확률 {cov/M:.4f}   부족분 {0.95-cov/M:.4f}"
          f"   (몬테카를로 오차 {np.sqrt((cov/M)*(1-cov/M)/M):.4f})")

n=  20: 포함확률 0.8980   부족분 0.0520   (몬테카를로 오차 0.0068)
n=  50: 포함확률 0.9270   부족분 0.0230   (몬테카를로 오차 0.0058)
n= 100: 포함확률 0.9405   부족분 0.0095   (몬테카를로 오차 0.0053)
n= 200: 포함확률 0.9420   부족분 0.0080   (몬테카를로 오차 0.0052)

# 백분위수법 — 가장 간단하고 앞에서 쓴 방법이다.
ci_percentile = (bootstrap_means.quantile(0.025), bootstrap_means.quantile(0.975))

# BCa — 편향과 가속을 보정한다. scipy 가 이미 구현해 두었으므로
# 직접 짤 필요가 없다. vectorized=True 는 통계량 함수가 축 인자를
# 받는다는 뜻이고, 그래야 속도가 난다.
from scipy.stats import bootstrap

def statistic(x, axis=-1):
    return np.mean(x, axis=axis)

result = bootstrap((original_sample,), statistic, n_resamples=5000,
                   method='bca', vectorized=True)
ci_bca = result.confidence_interval

from scipy import stats

n = len(original_sample)
g1 = stats.skew(original_sample)
a = g1 / (6 * np.sqrt(n))
bm = original_sample[np.random.default_rng(0).integers(0, n, (20000, n))].mean(axis=1)
z0 = stats.norm.ppf(np.mean(bm < original_sample.mean()))
za, z1 = stats.norm.ppf(0.025), stats.norm.ppf(0.975)
p_lo = stats.norm.cdf(z0 + (z0 + za) / (1 - a * (z0 + za)))
p_hi = stats.norm.cdf(z0 + (z0 + z1) / (1 - a * (z0 + z1)))

print(f"표본평균 {original_sample.mean():,.0f},  왜도 {g1:.4f}")
print(f"a = g1/(6 sqrt n) = {a:.6f},   z0 (B=20000) = {z0:.6f}")
print(f"조정된 백분위점 {100 * p_lo:.2f}% / {100 * p_hi:.2f}%  (보정 없으면 2.5% / 97.5%)")
print(f"백분위수 [{ci_percentile[0]:,.0f}, {ci_percentile[1]:,.0f}]"
      f"   폭 {ci_percentile[1] - ci_percentile[0]:,.0f}")
print(f"BCa      [{ci_bca.low:,.0f}, {ci_bca.high:,.0f}]"
      f"   폭 {ci_bca.high - ci_bca.low:,.0f}")

표본평균 67,846,  왜도 2.1166
a = g1/(6 sqrt n) = 0.078880,   z0 (B=20000) = 0.076527
조정된 백분위점 5.90% / 99.38%  (보정 없으면 2.5% / 97.5%)
백분위수 [47,026, 95,545]   폭 48,519
BCa      [50,525, 105,939]   폭 55,414

import numpy as np
from scipy import stats
rng = np.random.default_rng(3)
n, B = 50, 50000
x = rng.normal(100, 15, n)
bs = x[rng.integers(0, n, (B, n))].mean(axis=1)

for lvl in (0.90, 0.95, 0.99):
    a = (1 - lvl) / 2
    lo, hi = np.percentile(bs, [100*a, 100*(1-a)])
    print(f"{lvl:.0%}: [{lo:.3f}, {hi:.3f}]  폭 {hi-lo:.3f}")

90%: [96.570, 104.157]  폭 7.587
95%: [95.837, 104.895]  폭 9.058
99%: [94.350, 106.417]  폭 12.067

import numpy as np
rng = np.random.default_rng(42)
pop = rng.exponential(50000, 200000) + 20000
tm, tmed = pop.mean(), np.median(pop)

def coverage(n, stat, target, M=1200, B=400):
    c = 0
    for _ in range(M):
        s = rng.choice(pop, size=n, replace=False)
        bs = stat(s[rng.integers(0, n, (B, n))], axis=1)
        lo, hi = np.percentile(bs, [2.5, 97.5])
        c += lo <= target <= hi
    return round(c / M, 3)

for n in (20, 50, 100, 400):
    print(n, coverage(n, np.mean, tm), coverage(n, np.median, tmed))

20 0.892 0.93
50 0.933 0.941
100 0.932 0.944
400 0.945 0.945

import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
n = 20
x = rng.normal(100, 15, n)

# 빈도주의 95% t 구간
se = x.std(ddof=1) / np.sqrt(n)
print(np.round(stats.t.interval(0.95, n-1, x.mean(), se), 3))

# Bayes: 무정보 사전분포 pi(mu, sigma^2) ∝ 1/sigma^2 아래의 사후 신용구간
# 사후 mu | x ~ t_{n-1}(xbar, s/sqrt(n))  -- 수치적으로 동일하다
print(np.round(stats.t.interval(0.95, n-1, x.mean(), se), 3))

[ 91.124 103.381]
[ 91.124 103.381]

import numpy as np
from scipy.stats import bootstrap
rng = np.random.default_rng(3)
n = 20
x = rng.exponential(50000, n) + 20000

# 직접 구현한 백분위수 구간
B = 20000
bs = x[rng.integers(0, n, (B, n))].mean(axis=1)
print("percentile:", np.round(np.percentile(bs, [2.5, 97.5]), 0))

# scipy 의 BCa
res = bootstrap((x,), np.mean, n_resamples=20000, method='bca',
                random_state=0)
print("BCa       :", np.round(
    [res.confidence_interval.low, res.confidence_interval.high], 0))

# scipy 의 percentile / basic 도 확인
for m in ('percentile', 'basic'):
    r = bootstrap((x,), np.mean, n_resamples=20000, method=m, random_state=0)
    print(f"{m:11s}:", np.round([r.confidence_interval.low,
                                 r.confidence_interval.high], 0))

percentile: [45433. 82402.]
BCa       : [47963. 86984.]
percentile : [45796. 82555.]
basic      : [42773. 79533.]