Bartlett 검정¶

import numpy as np
import scipy.optimize as opt
import scipy.stats as stats

# (1) 에서 유도한 닫힌 꼴: nu = n-1, C = 1 + 1/(2 nu).
seed, size = 1, 100
nu = size - 1
C = 1 + 1 / (2 * nu)
x = stats.norm(loc=0, scale=1).rvs(size, random_state=seed)

print(f"nu = {nu},  보정항 C = 1 + 1/(2nu) = {C:.7f}")
print(f"\n{'sigma_y':>8}{'bartlett':>12}{'closed':>12}{'p':>9}")
for scale in [1.00, 1.05, 1.10, 1.15, 1.20]:
    y = stats.norm(loc=1, scale=scale).rvs(size, random_state=seed)
    stat, pval = stats.bartlett(x, y)
    closed = 2 * nu * np.log((1 + scale ** 2) / (2 * scale)) / C
    print(f"{scale:>8.2f}{stat:>12.8f}{closed:>12.8f}{pval:>9.4f}")

# p = 0.05 경계. 닫힌 꼴을 거꾸로 푼다 — 격자 탐색이 아니다.
crit = stats.chi2(df=1).ppf(0.95)
g = lambda s: 2 * nu * np.log((1 + s ** 2) / (2 * s)) / C - crit
root = opt.brentq(g, 1.0001, 5.0)
print(f"\n임계값 chi2_0.95(1) = {crit:.4f}")
print(f"chi2_B = 임계값이 되는 sigma_y = {root:.4f}")

# 같은 경계를 등분산 F 검정에서 구하면 1/sqrt(F_0.025(99,99)) 다.
print(f"F 검정의 같은 경계             = {1 / np.sqrt(stats.f(nu, nu).ppf(0.025)):.4f}")

nu = 99,  보정항 C = 1 + 1/(2nu) = 1.0050505

 sigma_y    bartlett      closed        p
    1.00 -0.00000000  0.00000000   1.0000
    1.05  0.23439030  0.23439030   0.6283
    1.10  0.89344836  0.89344836   0.3445
    1.15  1.91785747  1.91785747   0.1661
    1.20  3.25635555  3.25635555   0.0711

임계값 chi2_0.95(1) = 3.8415
chi2_B = 임계값이 되는 sigma_y = 1.2191
F 검정의 같은 경계             = 1.2191

import numpy as np
from scipy import stats

rng = np.random.default_rng(1234)
B = 6_000
dists = {
    "정규": lambda n: rng.normal(0, 1, n),
    "t(5) 두꺼운 꼬리": lambda n: rng.standard_t(5, n),
    "균등 (가벼운 꼬리)": lambda n: rng.uniform(-1, 1, n),
    "지수 (치우침 2)": lambda n: rng.exponential(1, n),
    "로그정규 (치우침 6)": lambda n: np.exp(rng.normal(0, 1, n)),
}
print("k=3, 집단당 n=20, 분산이 실제로 모두 같음, 명목 0.05")
print(f"{'분포':>22s} {'바틀렛':>8s} {'레빈(평균)':>11s} {'브라운-포사이드':>14s}")
for lab, f in dists.items():
    a = b = c = 0
    for _ in range(B):
        gs = [f(20) for _ in range(3)]
        a += stats.bartlett(*gs).pvalue < 0.05
        b += stats.levene(*gs, center="mean").pvalue < 0.05
        c += stats.levene(*gs, center="median").pvalue < 0.05
    print(f"{lab:>22s} {a / B:8.4f} {b / B:11.4f} {c / B:14.4f}")

k=3, 집단당 n=20, 분산이 실제로 모두 같음, 명목 0.05
                    분포      바틀렛      레빈(평균)       브라운-포사이드
                    정규   0.0508      0.0547         0.0367
           t(5) 두꺼운 꼬리   0.2153      0.0632         0.0418
           균등 (가벼운 꼬리)   0.0027      0.0617         0.0293
            지수 (치우침 2)   0.3795      0.1817         0.0462
          로그정규 (치우침 6)   0.6742      0.2510         0.0370

import numpy as np
from scipy import stats

rng = np.random.default_rng(1234)
B = 6_000
print("k=3, n=20, 정규, σ=(1, 1, r), 명목 0.05")
print(f"{'r':>5s} {'바틀렛':>8s} {'레빈(평균)':>11s} {'브라운-포사이드':>14s}")
for r in [1.5, 2.0, 3.0]:
    a = b = c = 0
    for _ in range(B):
        gs = [rng.normal(0, 1, 20), rng.normal(0, 1, 20), rng.normal(0, r, 20)]
        a += stats.bartlett(*gs).pvalue < 0.05
        b += stats.levene(*gs, center="mean").pvalue < 0.05
        c += stats.levene(*gs, center="median").pvalue < 0.05
    print(f"{r:5.1f} {a / B:8.4f} {b / B:11.4f} {c / B:14.4f}")

k=3, n=20, 정규, σ=(1, 1, r), 명목 0.05
    r      바틀렛      레빈(평균)       브라운-포사이드
  1.5   0.4300      0.3942         0.3338
  2.0   0.8777      0.8167         0.7665
  3.0   0.9983      0.9943         0.9903

import numpy as np
from scipy import stats

def bartlett_manual(groups, correct=True):
    """바틀렛 검정. correct=False 면 보정인자 C 를 쓰지 않는다."""
    nu = np.array([len(g) - 1 for g in groups], float)
    s2 = np.array([g.var(ddof=1) for g in groups])
    N = nu.sum()
    sp2 = (nu * s2).sum() / N
    num = N * np.log(sp2) - (nu * np.log(s2)).sum()
    C = (1 + ((1 / nu).sum() - 1 / N) / (3 * (len(groups) - 1))
         if correct else 1.0)
    chi = num / C
    return chi, C, stats.chi2.sf(chi, len(groups) - 1)

rng = np.random.default_rng(31)
gs = [rng.normal(0, 1, 12), rng.normal(0, 1.6, 9), rng.normal(0, 2.2, 15)]
c, C, p = bartlett_manual(gs)
r = stats.bartlett(*gs)
print(f"자작 구현: chi² = {c:.6f}, p = {p:.6f}  (보정인자 C = {C:.6f})")
print(f"scipy   : chi² = {r.statistic:.6f}, p = {r.pvalue:.6f}")
c0, _, p0 = bartlett_manual(gs, correct=False)
print(f"보정 없이 : chi² = {c0:.6f}, p = {p0:.6f}")

print("\n보정인자 C 가 표본 크기에 따라 얼마나 커지는가 (k=3, 균형)")
print(f"{'집단당 n':>9s} {'C':>9s}")
for n in [3, 5, 10, 20, 50, 100]:
    nu = np.array([n - 1.0] * 3)
    N = nu.sum()
    print(f"{n:9d} {1 + ((1 / nu).sum() - 1 / N) / 6:9.5f}")

print("\n보정의 효과: 등분산인데 기각하는 비율 (k=3, 정규, B=8000)")
rng = np.random.default_rng(77)
B = 8_000
print(f"{'집단당 n':>9s} {'보정 있음':>10s} {'보정 없음':>10s}")
for n in [3, 5, 10, 20]:
    a = b = 0
    for _ in range(B):
        g = [rng.normal(0, 1, n) for _ in range(3)]
        a += bartlett_manual(g)[2] < 0.05
        b += bartlett_manual(g, correct=False)[2] < 0.05
    print(f"{n:9d} {a / B:10.4f} {b / B:10.4f}")

자작 구현: chi² = 5.444673, p = 0.065721  (보정인자 C = 1.042839)
scipy   : chi² = 5.444673, p = 0.065721
보정 없이 : chi² = 5.677918, p = 0.058487

보정인자 C 가 표본 크기에 따라 얼마나 커지는가 (k=3, 균형)
    집단당 n         C
        3   1.22222
        5   1.11111
       10   1.04938
       20   1.02339
       50   1.00907
      100   1.00449

보정의 효과: 등분산인데 기각하는 비율 (k=3, 정규, B=8000)
    집단당 n      보정 있음      보정 없음
        3     0.0459     0.0848
        5     0.0517     0.0685
       10     0.0494     0.0578
       20     0.0501     0.0532

import numpy as np
from scipy import stats

rng = np.random.default_rng(55)
print("k=2 에서 바틀렛과 F 검정 (분산비) 비교")
print(f"{'':4s} {'n1':>4s} {'n2':>4s} {'s1²/s2²':>9s} "
      f"{'바틀렛 p':>10s} {'F 검정 p':>10s}")
for i in range(5):
    n1, n2 = int(rng.integers(8, 30)), int(rng.integers(8, 30))
    x = rng.normal(0, 1, n1)
    y = rng.normal(0, rng.uniform(1, 2.2), n2)
    F = x.var(ddof=1) / y.var(ddof=1)
    pf = 2 * min(stats.f.cdf(F, n1 - 1, n2 - 1),
                 stats.f.sf(F, n1 - 1, n2 - 1))
    print(f"{i:4d} {n1:4d} {n2:4d} {F:9.4f} "
          f"{stats.bartlett(x, y).pvalue:10.4f} {pf:10.4f}")

B = 20_000
print("\n제1종 오류율 (n1=n2=15, B=20000)")
for lab, seed, f in [("정규", 56, lambda r, n: r.normal(0, 1, n)),
                     ("t(5)", 57, lambda r, n: r.standard_t(5, n))]:
    rng2 = np.random.default_rng(seed)
    a = b = 0
    for _ in range(B):
        x, y = f(rng2, 15), f(rng2, 15)
        a += stats.bartlett(x, y).pvalue < 0.05
        F = x.var(ddof=1) / y.var(ddof=1)
        b += 2 * min(stats.f.cdf(F, 14, 14), stats.f.sf(F, 14, 14)) < 0.05
    print(f"  {lab:6s} 바틀렛 {a / B:.4f}   F 검정 {b / B:.4f}")

k=2 에서 바틀렛과 F 검정 (분산비) 비교
       n1   n2   s1²/s2²      바틀렛 p     F 검정 p
   0   29   26    0.5143     0.0914     0.0900
   1   28   24    0.6797     0.3398     0.3344
   2   29   14    0.3577     0.0258     0.0220
   3    9   14    0.2171     0.0330     0.0369
   4   20   22    0.7977     0.6192     0.6246

제1종 오류율 (n1=n2=15, B=20000)
  정규     바틀렛 0.0515   F 검정 0.0516
  t(5)   바틀렛 0.1484   F 검정 0.1485

등분산 검정이 필요한가?
    │
    ├─ 분산분석의 가정 점검 목적
    │     └─→ 하지 말 것. 처음부터 웰치를 쓴다
    │         (사전검정의 역설 — 가정 페이지 연습문제 7)
    │
    └─ 분산 자체가 연구 질문
          │
          ├─ 정규성이 보장됨 ──→ 바틀렛 (또는 k=2 면 F 검정)
          │
          └─ 그 외 ──→ 브라운-포사이드
                          또는 부트스트랩 구간

집단별 표준편차: 1.02, 1.11, 2.94  (비 2.9배)
브라운-포사이드 검정: W = ..., p = ...
→ 등분산을 가정하지 않고 웰치 분산분석을 사용