붓스트랩 신뢰구간 시각화 (코드)¶

import numpy as np

def simulate_income_data(n=5000, seed=3):
    """오른쪽으로 치우친 소득 자료를 만든다.

    모집단을 우리가 만들었으므로 참 평균을 알고 있다. 뒤에서 신뢰구간이
    그 값을 정말 95% 담는지 세어 볼 수 있다.
    """
    rng = np.random.default_rng(seed)
    return rng.exponential(scale=50_000, size=n) + 20_000

population = simulate_income_data()
print(population.mean(), np.median(population))    # 69526  54542

# 표본은 20개뿐이다. 치우친 모집단에서 이만큼만 뽑으면 붓스트랩 구간의
# 실제 포함확률이 95%에 못 미친다 — 아래 모의실험에서 확인한다.
rng = np.random.default_rng(303)
sample = rng.choice(population, size=20, replace=False)

69525.69401218835 54541.90941528454

c, sc, N = 20_000.0, 50_000.0, len(population)
print(f"이론: 평균 = c + 1/lam = {c + sc:,.1f},  중앙값 = c + ln2/lam = {c + sc*np.log(2):,.1f},"
      f"  표준편차 = {sc:,.1f}")
print(f"모의 모집단: 평균 {population.mean():,.2f},  중앙값 {np.median(population):,.2f}")
se = sc / np.sqrt(N)
print(f"둘의 표준오차 = sigma/sqrt(N) = 1/(2 f(m) sqrt(N)) = {se:,.1f}")
print(f"벗어난 정도: 평균 {(population.mean() - (c + sc)) / se:+.2f} SE,"
      f"  중앙값 {(np.median(population) - (c + sc * np.log(2))) / se:+.2f} SE")

이론: 평균 = c + 1/lam = 70,000.0,  중앙값 = c + ln2/lam = 54,657.4,  표준편차 = 50,000.0
모의 모집단: 평균 69,525.69,  중앙값 54,541.91
둘의 표준오차 = sigma/sqrt(N) = 1/(2 f(m) sqrt(N)) = 707.1
벗어난 정도: 평균 -0.67 SE,  중앙값 -0.16 SE

def bootstrap_sampling_distribution(sample, n_bootstrap=20_000, rng=None):
    """표본평균의 붓스트랩 표집분포. 반복문 없이 한 번에 계산한다.

    (n_bootstrap, n) 모양의 색인 배열을 만들어 한꺼번에 뽑으면 파이썬
    반복문이 사라진다. 붓스트랩처럼 같은 일을 만 번 되풀이하는 계산에서
    속도가 크게 달라진다.
    """
    rng = rng or np.random.default_rng(0)
    n = len(sample)
    return sample[rng.integers(0, n, (n_bootstrap, n))].mean(axis=1)

from scipy import stats

bd = bootstrap_sampling_distribution(sample)
print(f"표본: 평균 {sample.mean():,.2f},  sigma-hat {sample.std(ddof=0):,.2f},"
      f"  왜도 {stats.skew(sample):.4f}")
print(f"붓스트랩 평균 {bd.mean():,.2f}  (예측 = 표본평균)")
print(f"붓스트랩 SD   {bd.std(ddof=1):,.2f}  (예측 sigma-hat/sqrt(n) = "
      f"{sample.std(ddof=0) / np.sqrt(20):,.2f},"
      f"  몬테카를로 오차 {sample.std(ddof=0) / np.sqrt(20) / np.sqrt(2 * 20000):.1f})")
print(f"붓스트랩 왜도 {stats.skew(bd):.4f}  (예측 g1/sqrt(n) = "
      f"{stats.skew(sample) / np.sqrt(20):.4f},"
      f"  몬테카를로 오차 {np.sqrt(6 / 20000):.4f})")
print(f"붓스트랩 분포의 범위 [{bd.min():,.0f}, {bd.max():,.0f}]"
      f"   자료의 범위 [{sample.min():,.0f}, {sample.max():,.0f}]")

표본: 평균 78,122.71,  sigma-hat 54,427.21,  왜도 1.3169
붓스트랩 평균 78,195.53  (예측 = 표본평균)
붓스트랩 SD   12,284.09  (예측 sigma-hat/sqrt(n) = 12,170.29,  몬테카를로 오차 60.9)
붓스트랩 왜도 0.3040  (예측 g1/sqrt(n) = 0.2945,  몬테카를로 오차 0.0173)
붓스트랩 분포의 범위 [38,491, 130,423]   자료의 범위 [24,550, 223,179]

def compute_confidence_intervals(bootstrap_dist):
    """여러 신뢰수준에서의 백분위수 신뢰구간.

    신뢰수준을 올리면 구간이 넓어진다. 확신을 더 얻는 대가로 말해 주는
    범위가 흐려지는 맞바꿈이다.
    """
    return {
        '90%': np.percentile(bootstrap_dist, [5, 95]),
        '95%': np.percentile(bootstrap_dist, [2.5, 97.5]),
        '99%': np.percentile(bootstrap_dist, [0.5, 99.5]),
    }

cis = compute_confidence_intervals(bd)
se_b = bd.std(ddof=1)
for lab, z in (('90%', 1.644854), ('95%', 1.959964), ('99%', 2.575829)):
    lo, hi = cis[lab]
    print(f"{lab}: [{lo:,.0f}, {hi:,.0f}]  폭 {hi-lo:,.0f}"
          f"   예측 2 z SE = {2*z*se_b:,.0f}   (비 {(hi-lo)/(2*z*se_b):.4f})")
w = {k: v[1] - v[0] for k, v in cis.items()}
print(f"폭의 비 99/90 = {w['99%']/w['90%']:.4f}   (정규이론 {2.575829/1.644854:.4f})")

90%: [59,118, 99,527]  폭 40,409   예측 2 z SE = 40,411   (비 1.0000)
95%: [56,169, 104,028]  폭 47,859   예측 2 z SE = 48,153   (비 0.9939)
99%: [50,628, 112,855]  폭 62,226   예측 2 z SE = 63,283   (비 0.9833)
폭의 비 99/90 = 1.5399   (정규이론 1.5660)

def simulate_coverage(population, true_mean, n=20, B=1000, N=2000, rng=None):
    """백분위수 붓스트랩 구간의 실제 포함확률을 센다.

    모집단에서 표본을 새로 뽑는 일을 N 번 되풀이하며, 그때마다 만든 구간이
    참 평균을 담는지 센다. 95%로 약속한 구간의 실제 성적표다.
    """
    rng = rng or np.random.default_rng(303)
    hits = 0
    for _ in range(N):
        s = rng.choice(population, size=n, replace=False)
        bm = s[rng.integers(0, n, (B, n))].mean(axis=1)
        lo, hi = np.percentile(bm, [2.5, 97.5])
        hits += lo <= true_mean <= hi
    return hits / N

cov = simulate_coverage(population, population.mean())
print(f"실제 포함확률 = {cov:.4f}"
      f"   (N=2000 의 몬테카를로 오차 {np.sqrt(cov * (1 - cov) / 2000):.4f})")

rng2 = np.random.default_rng(303)
lo_miss = hi_miss = 0
for _ in range(2000):
    s = rng2.choice(population, size=20, replace=False)
    bm = s[rng2.integers(0, 20, (1000, 20))].mean(axis=1)
    lo, hi = np.percentile(bm, [2.5, 97.5])
    if population.mean() > hi:
        hi_miss += 1
    elif population.mean() < lo:
        lo_miss += 1
print(f"구간이 참값보다 왼쪽에 (아래로 빗나감): {hi_miss / 2000:.4f}")
print(f"구간이 참값보다 오른쪽에 (위로 빗나감): {lo_miss / 2000:.4f}")

실제 포함확률 = 0.9105   (N=2000 의 몬테카를로 오차 0.0064)
구간이 참값보다 왼쪽에 (아래로 빗나감): 0.0705
구간이 참값보다 오른쪽에 (위로 빗나감): 0.0190

import numpy as np
rng = np.random.default_rng(303)
population = simulate_income_data()

for n in (20, 50, 100, 200):
    widths = []
    for _ in range(400):
        s = rng.choice(population, n, replace=False)
        b = s[rng.integers(0, n, (2000, n))].mean(axis=1)
        lo, hi = np.percentile(b, [2.5, 97.5])
        widths.append(hi - lo)
    w = np.mean(widths)
    print(f"n={n:3d}  mean width = {w:,.0f}   width*sqrt(n) = {w*np.sqrt(n):,.0f}")

n= 20  mean width = 40,078   width*sqrt(n) = 179,233
n= 50  mean width = 26,088   width*sqrt(n) = 184,469
n=100  mean width = 18,908   width*sqrt(n) = 189,081
n=200  mean width = 13,472   width*sqrt(n) = 190,526

import numpy as np
from scipy import stats
rng = np.random.default_rng(303)
population = simulate_income_data()
true_mean = population.mean()          # 69,526
z = stats.norm.ppf(0.975)

cp = cn = cb = 0
N = 2000
for _ in range(N):
    s = rng.choice(population, 20, replace=False)
    bm = s[rng.integers(0, 20, (1000, 20))].mean(axis=1)
    lo, hi = np.percentile(bm, [2.5, 97.5]); cp += lo <= true_mean <= hi
    cn += abs(s.mean() - true_mean) < z * bm.std(ddof=1)
    # BCa
    z0 = stats.norm.ppf(np.clip((bm < s.mean()).mean(), 1e-6, 1-1e-6))
    jk = (s.sum() - s) / 19; jm = jk.mean()
    num = ((jm-jk)**3).sum(); den = 6*(((jm-jk)**2).sum())**1.5
    a = num/den if den > 0 else 0.0
    p1 = stats.norm.cdf(z0 + (z0-z)/(1 - a*(z0-z)))
    p2 = stats.norm.cdf(z0 + (z0+z)/(1 - a*(z0+z)))
    l2, h2 = np.percentile(bm, [100*p1, 100*p2]); cb += l2 <= true_mean <= h2

import numpy as np
rng = np.random.default_rng(303)
population = simulate_income_data()
sample = np.random.default_rng(3).choice(population, 20, replace=False)

idx = rng.integers(0, 20, (20_000, 20))
boot_means   = sample[idx].mean(axis=1)
boot_medians = np.median(sample[rng.integers(0, 20, (20_000, 20))], axis=1)