두 분산에 대한 F-검정¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(2024)
M, n1, n2 = 20_000, 20, 20

def gen(name, size):
    if name == "정규":    return rng.standard_normal(size)
    if name == "균등":    return rng.uniform(-1, 1, size)
    if name == "t(5)":    return rng.standard_t(5, size)
    if name == "지수":    return rng.exponential(1, size)
    if name == "로그정규": return rng.lognormal(0, 1, size)

print(f"{'분포':>8s} {'F 검정':>8s} {'Levene':>8s} {'B-F':>8s} {'Bartlett':>9s}")
for name in ["정규", "균등", "t(5)", "지수", "로그정규"]:
    f = lev = bf = bar = 0
    for _ in range(M):
        x, y = gen(name, n1), gen(name, n2)     # 두 분산이 실제로 같다
        F = x.var(ddof=1) / y.var(ddof=1)
        p = 2 * min(stats.f.cdf(F, n1 - 1, n2 - 1),
                    stats.f.sf(F, n1 - 1, n2 - 1))
        f   += p < 0.05
        lev += stats.levene(x, y, center='mean').pvalue < 0.05
        bf  += stats.levene(x, y, center='median').pvalue < 0.05
        bar += stats.bartlett(x, y).pvalue < 0.05
    print(f"{name:>8s} {f / M:8.4f} {lev / M:8.4f} "
          f"{bf / M:8.4f} {bar / M:9.4f}")

      분포     F 검정   Levene      B-F  Bartlett
      정규   0.0491   0.0539   0.0386    0.0491
      균등   0.0046   0.0556   0.0344    0.0046
    t(5)   0.1590   0.0559   0.0395    0.1590
      지수   0.2677   0.1438   0.0465    0.2677
    로그정규   0.4813   0.1905   0.0391    0.4813

import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
M, n1, n2 = 50_000, 15, 15
a = b = c = 0
for _ in range(M):
    x, y = rng.standard_normal(n1), rng.standard_normal(n2)
    s1, s2 = x.var(ddof=1), y.var(ddof=1)

    if s1 >= s2:                        # 큰 쪽을 분자로
        F, d1, d2 = s1 / s2, n1 - 1, n2 - 1
    else:
        F, d1, d2 = s2 / s1, n2 - 1, n1 - 1

    a += F > stats.f.ppf(0.95, d1, d2)    # ① 상단 0.05  ← 흔한 실수
    b += F > stats.f.ppf(0.975, d1, d2)   # ② 상단 0.025 ← 올바름

    F_fixed = s1 / s2                     # ③ 분자를 고정한 정식 양측
    p = 2 * min(stats.f.cdf(F_fixed, n1 - 1, n2 - 1),
                stats.f.sf(F_fixed, n1 - 1, n2 - 1))
    c += p < 0.05

print(f"① 큰 쪽 분자 + 상단 0.05   기각률 {a / M:.4f}")
print(f"② 큰 쪽 분자 + 상단 0.025  기각률 {b / M:.4f}")
print(f"③ 분자 고정 + 정식 양측    기각률 {c / M:.4f}")

① 큰 쪽 분자 + 상단 0.05   기각률 0.1008
② 큰 쪽 분자 + 상단 0.025  기각률 0.0508
③ 분자 고정 + 정식 양측    기각률 0.0508

import numpy as np
from scipy import stats

s1sq, n1 = 25.0, 16
s2sq, n2 = 10.0, 21
F = s1sq / s2sq
d1, d2 = n1 - 1, n2 - 1
lo = F / stats.f.ppf(0.975, d1, d2)
hi = F / stats.f.ppf(0.025, d1, d2)
print(f"F = {F}")
print(f"σ1²/σ2² 의 95% 구간  ({lo:.4f}, {hi:.4f})")
print(f"σ1/σ2   의 95% 구간  ({np.sqrt(lo):.4f}, {np.sqrt(hi):.4f})")

F = 2.5
σ1²/σ2² 의 95% 구간  (0.9716, 6.8898)
σ1/σ2   의 95% 구간  (0.9857, 2.6248)

rng = np.random.default_rng(5)
M = 20_000
for name, gen in [("정규", lambda n: rng.standard_normal(n)),
                  ("t(5)", lambda n: rng.standard_t(5, n)),
                  ("균등", lambda n: rng.uniform(-1, 1, n))]:
    cov = 0
    for _ in range(M):
        x, y = gen(16), gen(21)          # 참 비 = 1
        f = x.var(ddof=1) / y.var(ddof=1)
        l = f / stats.f.ppf(0.975, 15, 20)
        h = f / stats.f.ppf(0.025, 15, 20)
        cov += (l <= 1 <= h)
    print(f"{name:>5s}  포함률 {cov / M:.4f}")

   정규  포함률 0.9497
 t(5)  포함률 0.8462
   균등  포함률 0.9940

from scipy import stats

def power_F(n1, n2, ratio, alpha=0.05):
    """참 분산비가 ratio 일 때 양측 F 검정의 검정력."""
    d1, d2 = n1 - 1, n2 - 1
    lo = stats.f.ppf(alpha / 2, d1, d2)
    hi = stats.f.ppf(1 - alpha / 2, d1, d2)
    return stats.f.cdf(lo / ratio, d1, d2) + stats.f.sf(hi / ratio, d1, d2)

print(f"{'군당 n':>7s} {'비=1.5':>8s} {'비=2':>8s} {'비=3':>8s} {'비=4':>8s}")
for n in [10, 20, 30, 50, 100, 200]:
    row = " ".join(f"{power_F(n, n, r):8.4f}" for r in [1.5, 2, 3, 4])
    print(f"{n:7d} {row}")

print()
for r in [1.5, 2, 3, 4]:
    for n in range(4, 5_000):
        if power_F(n, n, r) >= 0.80:
            print(f"분산비 {r} 탐지에 필요한 군당 표본 {n}개 (검정력 80%)")
            break

   군당 n    비=1.5      비=2      비=3      비=4
     10   0.0854   0.1584   0.3347   0.4964
     20   0.1352   0.3083   0.6441   0.8375
     30   0.1862   0.4478   0.8285   0.9559
     50   0.2880   0.6703   0.9673   0.9976
    100   0.5183   0.9293   0.9997   1.0000
    200   0.8140   0.9982   1.0000   1.0000

분산비 1.5 탐지에 필요한 군당 표본 193개 (검정력 80%)
분산비 2 탐지에 필요한 군당 표본 68개 (검정력 80%)
분산비 3 탐지에 필요한 군당 표본 28개 (검정력 80%)
분산비 4 탐지에 필요한 군당 표본 19개 (검정력 80%)

import numpy as np
from scipy import stats

rng = np.random.default_rng(31)
M = 20_000

print(f"{'n1,n2':>9s} {'σ2/σ1':>7s} {'합동 t':>8s} {'Welch':>8s} {'2단계':>8s}")
for n1, n2 in [(10, 10), (10, 30), (30, 10)]:
    for r in [1.0, 2.0, 4.0]:
        a = b = c = 0
        for _ in range(M):
            x = rng.standard_normal(n1)          # 두 평균이 같다 (H0 참)
            y = rng.standard_normal(n2) * r
            p_pool = stats.ttest_ind(x, y, equal_var=True).pvalue
            p_welch = stats.ttest_ind(x, y, equal_var=False).pvalue
            s1, s2 = x.var(ddof=1), y.var(ddof=1)
            F = s1 / s2
            p_F = 2 * min(stats.f.cdf(F, n1 - 1, n2 - 1),
                          stats.f.sf(F, n1 - 1, n2 - 1))
            p_two = p_welch if p_F < 0.05 else p_pool   # 2단계 절차
            a += p_pool < 0.05
            b += p_welch < 0.05
            c += p_two < 0.05
        print(f"{n1:4d},{n2:<4d} {r:7.1f} {a / M:8.4f} "
              f"{b / M:8.4f} {c / M:8.4f}")

    n1,n2   σ2/σ1     합동 t    Welch      2단계
  10,10       1.0   0.0527   0.0517   0.0524
  10,10       2.0   0.0530   0.0490   0.0502
  10,10       4.0   0.0609   0.0503   0.0507
  10,30       1.0   0.0488   0.0495   0.0503
  10,30       2.0   0.0089   0.0491   0.0392
  10,30       4.0   0.0029   0.0476   0.0476
  30,10       1.0   0.0512   0.0525   0.0530
  30,10       2.0   0.1499   0.0529   0.0717
  30,10       4.0   0.2402   0.0476   0.0483

분산을 비교하려는 목적은?
    │
    ├─ 평균 검정의 가정 확인 ──────→ 하지 마라. 웰치를 쓴다.
    │                                (연습문제 9)
    │
    ├─ 분산 자체가 연구 질문
    │      │
    │      ├─ 자료가 정규에 가까운가
    │      │      ├─ 예 ────────→ F 검정 + 분산비 신뢰구간
    │      │      └─ 아니오 ────→ 브라운·포사이드 + 붓스트랩 구간
    │      │
    │      └─ 표본이 충분한가 (분산비 2에 군당 68개)
    │             └─ 아니면 "결정 불가"로 보고
    │
    └─ 이상점·꼬리가 진짜 관심 ───→ 분위수 비교, 로버스트 산포 측도