편향보정 가속 붓스트랩 (BCa)¶

from scipy.stats import norm
z0, a = norm.ppf(0.62), 0.042
for z in (-1.959964, 1.959964):
    v = z0 + (z0 + z) / (1 - a * (z0 + z))
    print(round(v, 4), round(norm.cdf(v), 4))
# -1.2415  0.1072
#  2.8091  0.9975

-1.2415 0.1072
2.8091 0.9975

import numpy as np

print(f"z0 = {z0:.6f},  공식이 깨지는 a = 1/(z0+1.96) = {1 / (z0 + 1.959964):.6f}")

rng = np.random.default_rng(2)
x = rng.exponential(1.0, 15)
x = x / x.std(ddof=1) * np.sqrt(8.4)      # 표본분산을 8.4 로 맞춘다
print(f"s^2 = {x.var(ddof=1):.4f}")

B = 10_000
qs = {p: [] for p in (2.5, 10.72, 97.5, 99.75)}
for s in range(20):
    r = np.random.default_rng(100 + s)
    bs = np.var(x[r.integers(0, 15, (B, 15))], axis=1, ddof=1)
    for p in qs:
        qs[p].append(np.percentile(bs, p))
print(f"\n{'백분위점':>9}{'평균':>10}{'SD':>9}{'꼬리에 남는 복제값':>20}")
for p in (2.5, 10.72, 97.5, 99.75):
    v = np.array(qs[p])
    tail = min(p, 100 - p) / 100 * B
    print(f"{p:>9}{v.mean():>10.4f}{v.std(ddof=1):>9.4f}{tail:>20.0f}")

z0 = 0.305481,  공식이 깨지는 a = 1/(z0+1.96) = 0.441414
s^2 = 8.4000

     백분위점        평균       SD          꼬리에 남는 복제값
      2.5    2.7828   0.0376                 250
    10.72    4.1762   0.0390                1072
     97.5   13.9051   0.0886                 250
    99.75   16.5086   0.2320                  25

import numpy as np
from scipy.stats import norm
for z0 in (0.0, 0.1, 0.3):
    for a in (0.0, 0.05, 0.15):
        zl = -1.959964
        v = z0 + (z0 + zl) / (1 - a * (z0 + zl))
        print(f"z0={z0:.1f} a={a:.2f} -> alpha1={norm.cdf(v):.4f}")

z0=0.0 a=0.00 -> alpha1=0.0250
z0=0.0 a=0.05 -> alpha1=0.0371
z0=0.0 a=0.15 -> alpha1=0.0649
z0=0.1 a=0.00 -> alpha1=0.0392
z0=0.1 a=0.05 -> alpha1=0.0546
z0=0.1 a=0.15 -> alpha1=0.0878
z0=0.3 a=0.00 -> alpha1=0.0869
z0=0.3 a=0.05 -> alpha1=0.1088
z0=0.3 a=0.15 -> alpha1=0.1517

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

def ci_pct_bca(x, stat, B=800, alpha=0.05):
    n = len(x); th = stat(x)
    idx = rng.integers(0, n, (B, n))
    bs = np.array([stat(x[i]) for i in idx])
    pct = np.percentile(bs, [100*alpha/2, 100*(1-alpha/2)])
    z0 = stats.norm.ppf(np.clip((bs < th).mean(), 1e-6, 1-1e-6))
    jk = np.array([stat(np.delete(x, i)) for i in range(n)])
    d = jk.mean() - jk
    den = ((d**2).sum())**1.5
    a = (d**3).sum() / (6*den) if den > 0 else 0.0
    zl, zu = stats.norm.ppf(alpha/2), stats.norm.ppf(1-alpha/2)
    a1 = stats.norm.cdf(z0 + (z0+zl)/(1 - a*(z0+zl)))
    a2 = stats.norm.cdf(z0 + (z0+zu)/(1 - a*(z0+zu)))
    return pct, np.percentile(bs, [100*a1, 100*a2])

for n in (20, 50, 200):
    M = 800; cp = cb = 0
    for _ in range(M):
        x = rng.exponential(1, n)
        pct, bca = ci_pct_bca(x, lambda v: v.var(ddof=1))
        cp += pct[0] <= 1 <= pct[1]
        cb += bca[0] <= 1 <= bca[1]
    print(n, round(cp/M, 3), round(cb/M, 3))

20 0.686 0.756
50 0.82 0.852
200 0.892 0.904

import numpy as np
rng = np.random.default_rng(0)
x = rng.normal(0, 1, 41)
jk = np.array([np.median(np.delete(x, i)) for i in range(41)])
print(len(np.unique(np.round(jk, 10))))    # 3

3

rng2 = np.random.default_rng(0)
def acc(gen, stat, n=41, M=2000):
    out = []
    for _ in range(M):
        z = gen(n)
        jk = np.array([stat(np.delete(z, i)) for i in range(n)])
        d = jk.mean() - jk
        den = ((d**2).sum())**1.5
        out.append((d**3).sum() / (6*den) if den > 0 else 0.0)
    return np.array(out)

a_med = acc(lambda n: rng2.exponential(1, n), np.median)
a_mean = acc(lambda n: rng2.exponential(1, n), np.mean)
print("median: %.5f +- %.5f" % (a_med.mean(), a_med.std()))
print("mean  : %.5f +- %.5f" % (a_mean.mean(), a_mean.std()))

median: 0.00001 +- 0.00089
mean  : 0.04130 +- 0.01584

import numpy as np
from scipy import stats
rng = np.random.default_rng(6)
n = 40
x = rng.exponential(1, n)

# 두 호출이 같은 재표본을 쓰도록 인덱스를 미리 뽑아 공유한다
idx = rng.integers(0, n, (20000, n))

def bca(x, stat, idx, alpha=0.05):
    n = len(x); th = stat(x)
    bs = np.array([stat(x[i]) for i in idx])
    z0 = stats.norm.ppf((bs < th).mean())
    jk = np.array([stat(np.delete(x, i)) for i in range(n)])
    d = jk.mean() - jk
    a = (d**3).sum() / (6 * ((d**2).sum())**1.5)
    zl, zu = stats.norm.ppf(alpha/2), stats.norm.ppf(1-alpha/2)
    a1 = stats.norm.cdf(z0 + (z0+zl)/(1 - a*(z0+zl)))
    a2 = stats.norm.cdf(z0 + (z0+zu)/(1 - a*(z0+zu)))
    return np.percentile(bs, [100*a1, 100*a2]), z0, a

ci1, z01, a1_ = bca(x, np.mean, idx)                    # 평균
ci2, z02, a2_ = bca(x, lambda v: np.log(v.mean()), idx) # 로그 평균
print(np.round(ci1, 6), round(z01, 4), round(a1_, 4))
print(np.round(np.exp(ci2), 6), round(z02, 4), round(a2_, 4))

[0.624251 1.161511] 0.0392 0.0421
[0.624594 1.164064] 0.0392 0.044

# 위와 같은 자료·같은 난수열로 평균의 붓스트랩 복제값을 만든다
rng = np.random.default_rng(6)
n = 40
x = rng.exponential(1, n)
th = x.mean()
idx = rng.integers(0, n, (20000, n))
bs = x[idx].mean(axis=1)

# 같은 붓스트랩 복제값으로 기본 구간을 두 척도에서 만든다
lo, hi = np.percentile(bs, [2.5, 97.5])
print(np.round([2*th - hi, 2*th - lo], 5))            # 원척도
llo, lhi = np.percentile(np.log(bs), [2.5, 97.5])
lth = np.log(th)
print(np.round(np.exp([2*lth - lhi, 2*lth - llo]), 5)) # 로그척도 후 되돌림

[0.56742 1.08876]
[0.63586 1.18821]