붓스트랩 신뢰구간 방법 (코드)¶

def bootstrap_percentile_ci(data, statistic, n_boot=10_000, alpha=0.05, rng=None):
    """백분위수 붓스트랩 신뢰구간.

    붓스트랩 분포의 2.5·97.5 백분위점을 그대로 쓴다. 가장 간단하고
    직관적이지만, 통계량이 치우쳐 있거나 편향이 있으면 어긋난다.
    """
    rng = rng or np.random.default_rng(0)
    n = len(data)
    boot_stats = np.array([
        statistic(data[rng.integers(0, n, n)])
        for _ in range(n_boot)
    ])
    lo = np.percentile(boot_stats, 100 * alpha / 2)
    hi = np.percentile(boot_stats, 100 * (1 - alpha / 2))
    return lo, hi, boot_stats

import numpy as np

rng = np.random.default_rng(11)
d = rng.exponential(1.0, 30)
B = 10_000
bs = d[rng.integers(0, 30, (B, 30))].mean(axis=1)   # 복제값 한 벌을 돌려 쓴다
th = d.mean()

pc = np.percentile(bs, [2.5, 97.5])                 # 원 척도의 백분위수 구간
pc_log = np.percentile(np.log(bs), [2.5, 97.5])     # 로그 척도에서 구한 것
ba = np.array([2 * th - pc[1], 2 * th - pc[0]])     # 원 척도의 기본 구간
ba_log = np.array([2 * np.log(th) - pc_log[1],
                   2 * np.log(th) - pc_log[0]])

print(f"theta_hat = {th:.6f}")
print(f"백분위수 (원 척도)  = [{pc[0]:.6f}, {pc[1]:.6f}]")
print(f"백분위수 (로그->역) = [{np.exp(pc_log[0]):.6f}, {np.exp(pc_log[1]):.6f}]"
      f"   같은가 {np.allclose(pc, np.exp(pc_log))}")
print(f"기본     (원 척도)  = [{ba[0]:.6f}, {ba[1]:.6f}]")
print(f"기본     (로그->역) = [{np.exp(ba_log[0]):.6f}, {np.exp(ba_log[1]):.6f}]"
      f"   같은가 {np.allclose(ba, np.exp(ba_log))}")

theta_hat = 1.087874
백분위수 (원 척도)  = [0.732291, 1.487806]
백분위수 (로그->역) = [0.732291, 1.487806]   같은가 True
기본     (원 척도)  = [0.687943, 1.443457]
기본     (로그->역) = [0.795447, 1.616120]   같은가 False

def bootstrap_basic_ci(data, statistic, boot_stats, alpha=0.05):
    """기본(역백분위수) 붓스트랩 신뢰구간.

    백분위점을 추정값 둘레로 되비춘다. 붓스트랩 분포가 오른쪽으로 치우쳐
    있으면 구간은 왼쪽으로 늘어나는데, 이는 "추정값이 참값보다 크게 나오는
    경향이 있다면 구간을 아래쪽으로 넓혀야 한다"는 셈법에서 나온다.
    백분위수법과 정반대 방향으로 움직이는 것이 처음에는 어리둥절하다.
    """
    theta_hat = statistic(data)
    lo = 2 * theta_hat - np.percentile(boot_stats, 100 * (1 - alpha / 2))
    hi = 2 * theta_hat - np.percentile(boot_stats, 100 * alpha / 2)
    return lo, hi

rng = np.random.default_rng(5)
M, Bb = 200, 1000
print(f"{'n':>5}{'기본 하한 < 0':>16}{'백분위수 하한 < 0':>20}")
for n in (10, 20, 30, 50):
    neg_b = neg_p = 0
    for _ in range(M):
        x = rng.exponential(1.0, n)
        bsv = np.var(x[rng.integers(0, n, (Bb, n))], axis=1, ddof=1)
        t = x.var(ddof=1)
        p = np.percentile(bsv, [2.5, 97.5])
        neg_b += (2 * t - p[1]) < 0
        neg_p += p[0] < 0
    print(f"{n:>5}{neg_b / M:>16.3f}{neg_p / M:>20.3f}")

    n       기본 하한 < 0         백분위수 하한 < 0
   10           0.180               0.000
   20           0.175               0.000
   30           0.100               0.000
   50           0.065               0.000

def bootstrap_bca_ci(data, statistic, boot_stats, alpha=0.05):
    """BCa(편향보정 가속) 붓스트랩 신뢰구간.

    백분위점을 두 가지로 조정한다. z0 은 붓스트랩 분포가 추정값을 중심으로
    치우친 정도(편향)를, a 는 통계량의 분산이 참값에 따라 달라지는 정도
    (가속)를 잡는다. 셋 중 가장 정확하지만 계산이 가장 무겁다.
    """
    n = len(data)
    theta_hat = statistic(data)

    # 편향보정 z0: 붓스트랩 값 중 관측된 추정값보다 작은 것의 비율을
    # 정규 분위점으로 옮긴다. 치우침이 없으면 절반이라 z0 이 0 이 된다.
    z0 = stats.norm.ppf(np.mean(boot_stats < theta_hat))

    # 가속 a: 잭나이프로 구한다. 관측값을 하나씩 빼 가며 통계량을 계산해,
    # 그 값들의 왜도에서 얻는다.
    jack = np.array([statistic(np.delete(data, i)) for i in range(n)])
    jack_mean = jack.mean()
    a_num = np.sum((jack_mean - jack) ** 3)
    a_den = 6 * np.sum((jack_mean - jack) ** 2) ** 1.5
    a = a_num / a_den if a_den != 0 else 0.0

    # 두 보정을 반영해 백분위점을 옮긴다. z0=0, a=0 이면 원래 백분위수법과
    # 정확히 같아진다.
    z_alpha = stats.norm.ppf(alpha / 2)
    z_1alpha = stats.norm.ppf(1 - alpha / 2)

    p_lo = stats.norm.cdf(z0 + (z0 + z_alpha) / (1 - a * (z0 + z_alpha)))
    p_hi = stats.norm.cdf(z0 + (z0 + z_1alpha) / (1 - a * (z0 + z_1alpha)))

    lo = np.percentile(boot_stats, 100 * p_lo)
    hi = np.percentile(boot_stats, 100 * p_hi)
    return lo, hi, z0, a

from scipy import stats

def bca_parts(data, statistic, boot_stats, alpha=0.05):
    n = len(data); th = statistic(data)
    z0 = stats.norm.ppf(np.mean(boot_stats < th))
    jack = np.array([statistic(np.delete(data, i)) for i in range(n)])
    jm = jack.mean()
    a = np.sum((jm - jack) ** 3) / (6 * np.sum((jm - jack) ** 2) ** 1.5)
    za, z1 = stats.norm.ppf(alpha / 2), stats.norm.ppf(1 - alpha / 2)
    p_lo = stats.norm.cdf(z0 + (z0 + za) / (1 - a * (z0 + za)))
    p_hi = stats.norm.cdf(z0 + (z0 + z1) / (1 - a * (z0 + z1)))
    return (np.percentile(boot_stats, 100 * p_lo),
            np.percentile(boot_stats, 100 * p_hi), z0, a, p_lo, p_hi)

lo_b, hi_b, z0, a, p_lo, p_hi = bca_parts(d, np.mean, bs)
g1 = stats.skew(d)
print(f"자료 왜도 g1 = {g1:.6f},  a 닫힌 꼴 g1/(6 sqrt n) = {g1 / (6 * np.sqrt(30)):.8f}")
print(f"코드의 a = {a:.8f}")
print(f"z0 = {z0:.6f}   (1차근사 {g1 / (6 * np.sqrt(30)):.6f},"
      f"  z0 의 몬테카를로 오차 {np.sqrt(0.25 / B) / stats.norm.pdf(0):.4f})")
print(f"조정된 백분위점: {100 * p_lo:.3f}% 와 {100 * p_hi:.3f}%  (보정 없으면 2.5% 와 97.5%)")
print(f"BCa      = [{lo_b:.6f}, {hi_b:.6f}]")
print(f"백분위수 = [{pc[0]:.6f}, {pc[1]:.6f}]")
print(f"기본     = [{ba[0]:.6f}, {ba[1]:.6f}]")

자료 왜도 g1 = 1.357052,  a 닫힌 꼴 g1/(6 sqrt n) = 0.04129377
코드의 a = 0.04129377
z0 = 0.037357   (1차근사 0.041294,  z0 의 몬테카를로 오차 0.0125)
조정된 백분위점: 4.059% 와 98.659%  (보정 없으면 2.5% 와 97.5%)
BCa      = [0.768921, 1.548148]
백분위수 = [0.732291, 1.487806]
기본     = [0.687943, 1.443457]

import numpy as np
from scipy import stats

# 포아송 자료에 세 방법을 모두 적용해 구간을 견준다. 자료가 치우쳐
# 있으므로 세 구간이 조금씩 어긋난다.
rng = np.random.default_rng(0)
data = stats.poisson.rvs(3.5, size=80, random_state=42)
print(data.mean())          # 3.425

lo_p, hi_p, boots = bootstrap_percentile_ci(data, np.mean, rng=rng)
lo_b, hi_b = bootstrap_basic_ci(data, np.mean, boots)
lo_bca, hi_bca, z0, a = bootstrap_bca_ci(data, np.mean, boots)
print(z0, a)                # -0.0266  0.0126

for name, (lo, hi) in [("백분위수", (lo_p, hi_p)), ("기본", (lo_b, hi_b)),
                       ("BCa", (lo_bca, hi_bca))]:
    print(f"{name:>5}: [{lo:.4f}, {hi:.4f}]  폭 {hi - lo:.4f}")

3.425
-0.026573386823392654 0.012573560456423716
 백분위수: [3.0375, 3.8250]  폭 0.7875
   기본: [3.0250, 3.8125]  폭 0.7875
  BCa: [3.0375, 3.8250]  폭 0.7875

n = len(data)
theta_hat = data.mean()
sig2 = data.var(ddof=0)                 # 경험분포의 분산
g1 = stats.skew(data)

print(f"격자 간격 = 1/n = {1 / n:.4f}")
print(f"붓스트랩 SE 극한 = {np.sqrt(sig2 / n):.6f}   (모의 {boots.std(ddof=1):.6f},"
      f"  몬테카를로 오차 {np.sqrt(sig2 / n) / np.sqrt(2 * 10000):.6f})")
print(f"참 SE = sqrt(3.5/80) = {np.sqrt(3.5 / n):.6f}")
lo_z = theta_hat - 1.959964 * np.sqrt(sig2 / n)
hi_z = theta_hat + 1.959964 * np.sqrt(sig2 / n)
print(f"정규근사 구간 = [{lo_z:.6f}, {hi_z:.6f}]")
print(f"백분위수 구간 = [{lo_p:.6f}, {hi_p:.6f}]"
      f"   -> 격자 눈금으로 {lo_p * n:.0f}/{n} 와 {hi_p * n:.0f}/{n}")

print(f"\n자료 왜도 g1 = {g1:.6f}")
print(f"a 의 닫힌 꼴 g1/(6 sqrt n) = {g1 / (6 * np.sqrt(n)):.8f}")
below = np.mean(boots < theta_hat)
eq = np.mean(np.isclose(boots, theta_hat))
print(f"P*(theta* < theta_hat) = {below:.4f},  P*(= ) = {eq:.4f},"
      f"  P*(>) = {np.mean(boots > theta_hat):.4f}")
print(f"z0 (코드, 강한 부등호) = {stats.norm.ppf(below):+.6f}")
print(f"z0 (동점을 반씩 세면)  = {stats.norm.ppf(below + eq / 2):+.6f}")
print(f"z0 의 1차근사 = a = {g1 / (6 * np.sqrt(n)):+.6f}"
      f"   (z0 의 몬테카를로 오차 {np.sqrt(0.25 / 10000) / stats.norm.pdf(0):.4f})")

격자 간격 = 1/n = 0.0125
붓스트랩 SE 극한 = 0.199041   (모의 0.199015,  몬테카를로 오차 0.001407)
참 SE = sqrt(3.5/80) = 0.209165
정규근사 구간 = [3.034887, 3.815113]
백분위수 구간 = [3.037500, 3.825000]   -> 격자 눈금으로 243/80 와 306/80

자료 왜도 g1 = 0.674768
a 의 닫힌 꼴 g1/(6 sqrt n) = 0.01257356
P*(theta* < theta_hat) = 0.4894,  P*(= ) = 0.0237,  P*(>) = 0.4869
z0 (코드, 강한 부등호) = -0.026573
z0 (동점을 반씩 세면)  = +0.003133
z0 의 1차근사 = a = +0.012574   (z0 의 몬테카를로 오차 0.0125)

import numpy as np
rng = np.random.default_rng(0)
data = rng.exponential(scale=1.0, size=30)
print(data.mean())      # 1.1845

lo_p, hi_p, boots = bootstrap_percentile_ci(data, np.mean, rng=rng)
lo_b, hi_b = bootstrap_basic_ci(data, np.mean, boots)
lo_bca, hi_bca, z0, a = bootstrap_bca_ci(data, np.mean, boots)

1.1844536180110892

import numpy as np
from scipy import stats
rng = np.random.default_rng(42)

def run(n=20, M=2000, B=2000, df=3):
    cp = cbs = cb = ct = 0
    zl, zu = stats.norm.ppf(0.025), stats.norm.ppf(0.975)
    for _ in range(M):
        d = rng.chisquare(df, n); th = d.mean()
        b = d[rng.integers(0, n, (B, n))].mean(1)
        lo, hi = np.percentile(b, [2.5, 97.5])
        cp  += lo <= df <= hi
        cbs += (2*th - hi) <= df <= (2*th - lo)
        z0 = stats.norm.ppf(np.clip((b < th).mean(), 1e-6, 1-1e-6))
        jk = (d.sum() - d) / (n - 1); 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+zl)/(1 - a*(z0+zl)))
        p2 = stats.norm.cdf(z0 + (z0+zu)/(1 - a*(z0+zu)))
        l2, h2 = np.percentile(b, [100*p1, 100*p2]); cb += l2 <= df <= h2
        ct += stats.ttest_1samp(d, df).pvalue >= 0.05
    return cp/M, cbs/M, cb/M, ct/M