붓스트랩 분산 검정¶

import numpy as np

def variance_ratio(x1, x2):
    """두 표본분산의 비."""
    return np.var(x1, ddof=1) / np.var(x2, ddof=1)

def bootstrap_varratio(x1, x2, B=2000, seed=None):
    """분산비의 붓스트랩 신뢰구간과 양측 p-값을 구한다.

    F 검정과 달리 정규성을 가정하지 않는다. 각 표본에서 따로 복원추출해
    비를 다시 구하는 일을 B 번 되풀이한다.
    """
    rng = np.random.default_rng(seed)
    x1 = np.asarray(x1, dtype=float)
    x2 = np.asarray(x2, dtype=float)
    n1, n2 = len(x1), len(x2)
    stat_obs = variance_ratio(x1, x2)

    boots = np.empty(B)
    for b in range(B):
        b1 = rng.choice(x1, size=n1, replace=True)
        b2 = rng.choice(x2, size=n2, replace=True)
        boots[b] = np.log(variance_ratio(b1, b2))

    # 비는 아래로 0, 위로 무한이라 분포가 치우친다. 로그를 씌우면 그 눈금이
    # 대칭에 가까워지므로 구간을 로그 눈금에서 만든 뒤 되돌린다.
    lo, hi = np.percentile(boots, [2.5, 97.5])
    ci = (float(np.exp(lo)), float(np.exp(hi)))

    # p-값은 붓스트랩 분포를 관측값이 아니라 귀무값 log(1)=0 과 견주어
    # 만든다. 관측값과 견주면 언제나 0.5 근처가 나와 뜻이 없다.
    p_two = 2 * min(np.mean(boots <= 0.0), np.mean(boots >= 0.0))
    p_two = float(min(p_two, 1.0))

    return float(stat_obs), ci, p_two

# === 구간 뒤집기와 p 값이 정말 같은 판정을 주는가 ===
x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
x2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)

# 함수 안의 붓스트랩 값을 다시 만들어 들여다본다 (같은 씨앗, 같은 수열).
rng = np.random.default_rng(42)
B = 10000
boots = np.empty(B)
for b in range(B):
    b1 = rng.choice(x1, size=len(x1), replace=True)
    b2 = rng.choice(x2, size=len(x2), replace=True)
    boots[b] = np.log(variance_ratio(b1, b2))

lo_cnt, hi_cnt = np.mean(boots <= 0.0), np.mean(boots >= 0.0)
ties = int((boots == 0.0).sum())
print(f"P(log theta* <= 0) = {lo_cnt:.4f}")
print(f"P(log theta* >= 0) = {hi_cnt:.4f}")
print(f"두 비율의 합        = {lo_cnt + hi_cnt:.4f}   (1 + 동점비율)")
print(f"정확히 0 인 값      = {ties} 개 / {B}   -> 동점비율 {ties / B:.4f}")
print(f"2 x min            = {2 * min(lo_cnt, hi_cnt):.4f}  (자르기 전)")

p = min(2 * min(lo_cnt, hi_cnt), 1.0)
print(f"\n   alpha     lower     upper   1 in CI   p<alpha")
for a in (0.05, 0.10, 0.15, 0.1892, 0.19, 0.25, 0.40):
    lo, hi = np.exp(np.percentile(boots, [100 * a / 2, 100 * (1 - a / 2)]))
    inci = "yes" if lo <= 1.0 <= hi else "no"
    rej = "yes" if p < a else "no"
    print(f"{a:>8.4f}{lo:>10.4f}{hi:>10.4f}{inci:>10}{rej:>10}")

P(log theta* <= 0) = 0.9118
P(log theta* >= 0) = 0.0946
두 비율의 합        = 1.0064   (1 + 동점비율)
정확히 0 인 값      = 64 개 / 10000   -> 동점비율 0.0064
2 x min            = 0.1892  (자르기 전)

   alpha     lower     upper   1 in CI   p<alpha
  0.0500    0.1345    1.6092       yes        no
  0.1000    0.1705    1.2526       yes        no
  0.1500    0.1946    1.0787       yes        no
  0.1892    0.2113    0.9967        no        no
  0.1900    0.2117    0.9952        no       yes
  0.2500    0.2381    0.8930        no       yes
  0.4000    0.2904    0.7458        no       yes

import numpy as np
from scipy.special import polygamma
from scipy.stats import f as f_dist

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
x2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)

theta_hat, ci, p = bootstrap_varratio(x1, x2, B=10000, seed=42)
print(f"Observed ratio: {theta_hat:.4f}")
print(f"95% Bootstrap CI: ({ci[0]:.4f}, {ci[1]:.4f})")
print(f"Bootstrap p-value: {p:.4f}")

# === 정규이론 F 구간과 견준다 ===
n = len(x1)
fl, fh = theta_hat / f_dist(n - 1, n - 1).ppf([0.975, 0.025])
p_f = 2 * min(f_dist(n - 1, n - 1).cdf(theta_hat), f_dist(n - 1, n - 1).sf(theta_hat))
print(f"\n정규이론 F 구간 : ({fl:.4f}, {fh:.4f})   로그 폭 {np.log(fh / fl):.4f}   p = {p_f:.4f}")
print(f"붓스트랩 구간   : ({ci[0]:.4f}, {ci[1]:.4f})   로그 폭 {np.log(ci[1] / ci[0]):.4f}   p = {p:.4f}")
print(f"붓스트랩이 {100 * (1 - np.log(ci[1] / ci[0]) / np.log(fh / fl)):.1f} % 좁다")

# === 왜 좁은가: log theta* 의 산포를 정규이론 참값과 견준다 ===
# Var(log F(d,d)) = psi'(d/2) + psi'(d/2) 가 정확한 값이다.
v_exact = 2 * polygamma(1, (n - 1) / 2)
rng = np.random.default_rng(42)
boots = np.empty(10000)
for b in range(10000):
    b1 = rng.choice(x1, size=n, replace=True)
    b2 = rng.choice(x2, size=n, replace=True)
    boots[b] = np.log(variance_ratio(b1, b2))
sd_boot = boots.std(ddof=1)
print(f"\nVar(log F(7,7)) = 2 psi'(3.5) = {v_exact:.6f}  ->  SD = {np.sqrt(v_exact):.4f}  (정확한 값)")
print(f"붓스트랩 SD(log theta*)                 = {sd_boot:.4f}")
print(f"비 = {sd_boot / np.sqrt(v_exact):.4f}   (붓스트랩이 산포를 그만큼 적게 본다)")

# === 범인은 경험 첨도다. 그리고 그것에는 n 이 정한 천장이 있다 ===
def emp_kurtosis(v):
    v = np.asarray(v, float); m = v.mean()
    return ((v - m) ** 4).mean() / (((v - m) ** 2).mean()) ** 2

k1, k2 = emp_kurtosis(x1), emp_kurtosis(x2)
bound = (n * n - 3 * n + 3) / (n - 1)
print(f"\n경험 첨도  x1 = {k1:.4f},  x2 = {k2:.4f}   (정규는 3)")
print(f"n = {n} 인 표본이 가질 수 있는 최대 첨도 = (n^2-3n+3)/(n-1) = {bound:.4f}")
z = np.array([np.sqrt(n - 1)] + [-1 / np.sqrt(n - 1)] * (n - 1))
print(f"  그 최대를 내는 꼴: 한 점 {z[0]:+.4f}, 나머지 {n - 1} 점 {z[1]:+.4f}")
print(f"  그 꼴의 첨도 = {emp_kurtosis(z):.4f}   (공식과 일치)")
print(f"\n5.3절의 델타 공식 SD(log theta) ~ sqrt((k1-1)/n + (k2-1)/n)")
print(f"  경험 첨도를 넣으면 {np.sqrt((k1 - 1) / n + (k2 - 1) / n):.4f}")
print(f"  정규값 3 을 넣으면 {np.sqrt(2 * (3 - 1) / n):.4f}")
print(f"  배율 sqrt((kbar-1)/2) = {np.sqrt(((k1 + k2) / 2 - 1) / 2):.4f}   (관측된 비 {sd_boot / np.sqrt(v_exact):.4f})")
for nn in (8, 20, 30, 100):
    print(f"  n = {nn:>3}: 첨도 천장 {(nn * nn - 3 * nn + 3) / (nn - 1):8.3f}")

Observed ratio: 0.4732
95% Bootstrap CI: (0.1345, 1.6092)
Bootstrap p-value: 0.1892

정규이론 F 구간 : (0.0947, 2.3637)   로그 폭 3.2168   p = 0.3448
붓스트랩 구간   : (0.1345, 1.6092)   로그 폭 2.4816   p = 0.1892
붓스트랩이 22.9 % 좁다

Var(log F(7,7)) = 2 psi'(3.5) = 0.660716  ->  SD = 0.8128  (정확한 값)
붓스트랩 SD(log theta*)                 = 0.6136
비 = 0.7549   (붓스트랩이 산포를 그만큼 적게 본다)

경험 첨도  x1 = 1.8985,  x2 = 1.7619   (정규는 3)
n = 8 인 표본이 가질 수 있는 최대 첨도 = (n^2-3n+3)/(n-1) = 6.1429
  그 최대를 내는 꼴: 한 점 +2.6458, 나머지 7 점 -0.3780
  그 꼴의 첨도 = 6.1429   (공식과 일치)

5.3절의 델타 공식 SD(log theta) ~ sqrt((k1-1)/n + (k2-1)/n)
  경험 첨도를 넣으면 0.4556
  정규값 3 을 넣으면 0.7071
  배율 sqrt((kbar-1)/2) = 0.6443   (관측된 비 0.7549)
  n =   8: 첨도 천장    6.143
  n =  20: 첨도 천장   18.053
  n =  30: 첨도 천장   28.034
  n = 100: 첨도 천장   98.010

import numpy as np

def boot_log_ratios(x1, x2, B, rng):
    n1, n2 = len(x1), len(x2)
    i1 = rng.integers(0, n1, (B, n1))
    i2 = rng.integers(0, n2, (B, n2))
    return np.log(x1[i1].var(axis=1, ddof=1) / x2[i2].var(axis=1, ddof=1))

rng = np.random.default_rng(0)
covers = 0
for _ in range(500):
    x1 = rng.normal(0, 1, 30)
    x2 = rng.normal(0, 1, 30)
    bs = boot_log_ratios(x1, x2, 5000, rng)
    lo, hi = np.percentile(bs, [2.5, 97.5])
    if np.exp(lo) <= 1.0 <= np.exp(hi):
        covers += 1

print(f"Coverage: {covers/500:.3f}")

Coverage: 0.920

import numpy as np
from scipy.stats import f as f_dist

def p_vs_null(x1, x2, B, rng):
    n1, n2 = len(x1), len(x2)
    i1 = rng.integers(0, n1, (B, n1))
    i2 = rng.integers(0, n2, (B, n2))
    bs = np.log(x1[i1].var(axis=1, ddof=1) / x2[i2].var(axis=1, ddof=1))
    return min(2 * min((bs <= 0).mean(), (bs >= 0).mean()), 1.0)

rng = np.random.default_rng(7)
for name, gen in [("Normal", lambda n: rng.normal(0, 1, n)),
                  ("Exponential", lambda n: rng.exponential(1, n))]:
    for n in [20, 50]:
        rej = sum(p_vs_null(gen(n), gen(n), 1000, rng) < 0.05
                  for _ in range(2000))
        print(f"{name:12s} n={n}: size = {rej/2000:.4f}")

Normal       n=20: size = 0.0725
Normal       n=50: size = 0.0610
Exponential  n=20: size = 0.0915
Exponential  n=50: size = 0.0905

import numpy as np
from scipy.stats import f as f_dist

def variance_ratio(x1, x2):
    return np.var(x1, ddof=1) / np.var(x2, ddof=1)

def bootstrap_pvalue(x1, x2, B, rng):
    """양측 붓스트랩 p-값. 관측값이 아니라 귀무값 1 과 견주어 만든다."""
    n1, n2 = len(x1), len(x2)
    i1 = rng.integers(0, n1, (B, n1))
    i2 = rng.integers(0, n2, (B, n2))
    bs = np.log(x1[i1].var(axis=1, ddof=1) / x2[i2].var(axis=1, ddof=1))
    return min(2 * min((bs <= 0).mean(), (bs >= 0).mean()), 1.0)

rng = np.random.default_rng(42)
rej_f, rej_boot = 0, 0
n_sims = 2000

for _ in range(n_sims):
    x1 = rng.exponential(1, 20)
    x2 = rng.exponential(1, 20)
    F = variance_ratio(x1, x2)
    p_f = 2 * min(f_dist.cdf(F, 19, 19), f_dist.sf(F, 19, 19))
    if p_f < 0.05:
        rej_f += 1
    if bootstrap_pvalue(x1, x2, 1000, rng) < 0.05:
        rej_boot += 1

print(f"F-test false positive rate:    {rej_f/n_sims:.4f}")
print(f"Bootstrap false positive rate: {rej_boot/n_sims:.4f}")

F-test false positive rate:    0.2695
Bootstrap false positive rate: 0.0890

import numpy as np
from scipy import stats as sp_stats

def variance_ratio(x1, x2):
    return np.var(x1, ddof=1) / np.var(x2, ddof=1)

x1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
x2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)
rng = np.random.default_rng(42)
B = 10000
n1, n2 = len(x1), len(x2)

log_obs = np.log(variance_ratio(x1, x2))
boots = np.array([np.log(variance_ratio(
    rng.choice(x1, n1, replace=True),
    rng.choice(x2, n2, replace=True))) for _ in range(B)])

# 백분위수법
lo_p, hi_p = np.percentile(boots, [2.5, 97.5])

# BCa: 편향보정 항
z0 = sp_stats.norm.ppf(np.mean(boots < log_obs))

# 가속 항. 1집단에 잭나이프를 적용해 구한다
jk = np.empty(n1)
for i in range(n1):
    x1_jk = np.delete(x1, i)
    jk[i] = np.log(np.var(x1_jk, ddof=1) / np.var(x2, ddof=1))
jk_mean = jk.mean()
a_hat = np.sum((jk_mean - jk)**3) / (6 * np.sum((jk_mean - jk)**2)**1.5)

z_lo, z_hi = sp_stats.norm.ppf(0.025), sp_stats.norm.ppf(0.975)
alpha1 = sp_stats.norm.cdf(z0 + (z0 + z_lo) / (1 - a_hat * (z0 + z_lo)))
alpha2 = sp_stats.norm.cdf(z0 + (z0 + z_hi) / (1 - a_hat * (z0 + z_hi)))
lo_bca, hi_bca = np.percentile(boots, [100*alpha1, 100*alpha2])

print(f"z0 = {z0:.4f}, a_hat = {a_hat:.4f}")
print(f"BCa percentiles: {100*alpha1:.2f}%, {100*alpha2:.2f}%")
print(f"Percentile CI (ratio): ({np.exp(lo_p):.3f}, {np.exp(hi_p):.3f})")
print(f"BCa CI (ratio):        ({np.exp(lo_bca):.3f}, {np.exp(hi_bca):.3f})")

z0 = 0.0175, a_hat = 0.0560
BCa percentiles: 4.14%, 98.75%
Percentile CI (ratio): (0.135, 1.609)
BCa CI (ratio):        (0.160, 2.117)