등분산에 대한 F-검정¶

import numpy as np
import scipy.stats as stats

def f_test(data_0, data_1):
    """두 분산이 같은지 F 검정한다.

    두 표본분산의 비가 F 분포를 따른다는 사실을 쓴다. 카이제곱과 마찬가지로
    좌우가 대칭이 아니므로 작은 쪽 꼬리를 두 배 해 양측 p-값을 만든다.
    이 검정은 정규성에 매우 민감하다 — 그래서 실제로는 Levene 을 더 쓴다.
    """
    statistic = data_0.var(ddof=1) / data_1.var(ddof=1)
    df1 = data_0.shape[0] - 1
    df2 = data_1.shape[0] - 1
    p_value = 2 * min(
        stats.f(df1, df2).cdf(statistic),
        stats.f(df1, df2).sf(statistic)
    )
    return statistic, p_value

# 표본 크기를 일부러 다르게 두어 자유도가 뒤바뀌는 것까지 보이게 한다.
rng = np.random.default_rng(7)
a = rng.normal(scale=1.0, size=20)
b = rng.normal(scale=1.7, size=35)

Fab, pab = f_test(a, b)
Fba, pba = f_test(b, a)
print(f"f_test(a, b):  F = {Fab:.6f}   p = {pab:.6f}   (자유도 19, 34)")
print(f"f_test(b, a):  F = {Fba:.6f}   p = {pba:.6f}   (자유도 34, 19)")
print(f"F 는 역수인가 ?      {np.isclose(Fab * Fba, 1.0)}")
print(f"p 는 같은가 ?        {np.isclose(pab, pba)}")

# 임계값의 역수 관계. F_{1-a/2}(d1,d2) = 1 / F_{a/2}(d2,d1).
d1, d2 = 19, 34
hi = stats.f(d1, d2).ppf(0.975)
lo_swapped = stats.f(d2, d1).ppf(0.025)
print(f"\nF_0.975(19,34)        = {hi:.6f}")
print(f"1 / F_0.025(34,19)    = {1 / lo_swapped:.6f}")

# 양측 p-값이 1 을 넘지 않는 까닭은 d1 = d2 일 때 가장 또렷하다.
# F(d,d) 는 역수에 대해 대칭이므로 중앙값이 정확히 1 이다.
print(f"\nP(F(99,99) <= 1)      = {stats.f(99, 99).cdf(1):.6f}")

f_test(a, b):  F = 0.237062   p = 0.001581   (자유도 19, 34)
f_test(b, a):  F = 4.218313   p = 0.001581   (자유도 34, 19)
F 는 역수인가 ?      True
p 는 같은가 ?        True

F_0.975(19,34)        = 2.153050
1 / F_0.025(34,19)    = 2.153050

P(F(99,99) <= 1)      = 0.500000

# Bartlett 검정과 같은 설정으로 돌려 두 검정의 p-값을 견주어 볼 수 있다.
seed, size = 1, 100
x = stats.norm(loc=0, scale=1).rvs(size, random_state=seed)

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 = f_test(x, y)
    print(f"sigma_y={scale:.2f}: F={stat:.2f}, p={pval:.3f}")

sigma_y=1.00: F=1.00, p=1.000
sigma_y=1.05: F=0.91, p=0.628
sigma_y=1.10: F=0.83, p=0.345
sigma_y=1.15: F=0.76, p=0.166
sigma_y=1.20: F=0.69, p=0.071

import numpy as np
from scipy import stats

def F_p(x, y):
    F = x.var(ddof=1) / y.var(ddof=1)
    d1, d2 = len(x) - 1, len(y) - 1
    return 2 * min(stats.f.cdf(F, d1, d2), stats.f.sf(F, d1, d2))

rng = np.random.default_rng(606)
B = 15_000
dists = {
    "정규": lambda n: rng.normal(0, 1, n),
    "균등": lambda n: rng.uniform(-np.sqrt(3), np.sqrt(3), n),
    "t(8)": lambda n: rng.standard_t(8, n),
    "t(5)": lambda n: rng.standard_t(5, n),
    "지수": lambda n: rng.exponential(1, n),
    "로그정규": lambda n: np.exp(rng.normal(0, 1, n)),
}
print("두 집단 모두 같은 분포에서 (n=25 씩), 명목 0.05")
print(f"{'분포':>20s} {'F 검정(분산)':>12s} {'웰치 t(평균)':>12s}")
for lab, f in dists.items():
    a = b = 0
    for _ in range(B):
        x, y = f(25), f(25)
        a += F_p(x, y) < 0.05
        b += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
    print(f"{lab:>20s} {a / B:12.4f} {b / B:12.4f}")

두 집단 모두 같은 분포에서 (n=25 씩), 명목 0.05
                  분포     F 검정(분산)     웰치 t(평균)
                  정규       0.0529       0.0505
                  균등       0.0043       0.0493
                t(8)       0.1115       0.0481
                t(5)       0.1653       0.0481
                  지수       0.2770       0.0477
                로그정규       0.5049       0.0355

import numpy as np
from scipy import stats

def ratio_ci(x, y, a=0.05):
    d1, d2 = len(x) - 1, len(y) - 1
    R = x.var(ddof=1) / y.var(ddof=1)
    return R / stats.f.ppf(1 - a / 2, d1, d2), R / stats.f.ppf(a / 2, d1, d2)

rng = np.random.default_rng(707)
B = 10_000
print("참 분산비 = 1 을 덮는 비율 (명목 95%)")
print(f"{'분포':>12s} {'n=25':>8s} {'n=100':>8s}")
for lab, f in {"정규": lambda n: rng.normal(0, 1, n),
               "t(8)": lambda n: rng.standard_t(8, n),
               "t(5)": lambda n: rng.standard_t(5, n),
               "지수": lambda n: rng.exponential(1, n),
               "균등": lambda n: rng.uniform(-1, 1, n)}.items():
    row = []
    for n in [25, 100]:
        c = 0
        for _ in range(B):
            lo, hi = ratio_ci(f(n), f(n))
            c += lo <= 1.0 <= hi
        row.append(c / B)
    print(f"{lab:>12s} {row[0]:8.4f} {row[1]:8.4f}")

참 분산비 = 1 을 덮는 비율 (명목 95%)
          분포     n=25    n=100
          정규   0.9503   0.9495
        t(8)   0.8971   0.8786
        t(5)   0.8345   0.7848
          지수   0.7272   0.6852
          균등   0.9955   0.9969

import numpy as np
from scipy import stats

rng = np.random.default_rng(707)

def perm_p(x, y, rng, R=600):
    """로그 분산비의 순열검정. 중앙값으로 중심을 맞춘 뒤 섞는다."""
    obs = abs(np.log(x.var(ddof=1) / y.var(ddof=1)))
    z = np.concatenate([x - np.median(x), y - np.median(y)])
    n = len(x)
    cnt = 0
    for _ in range(R):
        p = rng.permutation(z)
        cnt += abs(np.log(p[:n].var(ddof=1) / p[n:].var(ddof=1))) >= obs
    return (cnt + 1) / (R + 1)

B = 4_000
print("두 집단의 분산이 실제로 같음, n=25 씩, 명목 0.05")
print(f"{'분포':>12s} {'F 검정':>8s} {'브라운-포사이드':>14s} {'순열(로그 분산비)':>16s}")
for lab, f in {"정규": lambda n: rng.normal(0, 1, n),
               "t(5)": lambda n: rng.standard_t(5, n),
               "지수": lambda n: rng.exponential(1, n),
               "로그정규": lambda n: np.exp(rng.normal(0, 1, n))}.items():
    a = b = c = 0
    for _ in range(B):
        x, y = f(25), f(25)
        F = x.var(ddof=1) / y.var(ddof=1)
        a += 2 * min(stats.f.cdf(F, 24, 24), stats.f.sf(F, 24, 24)) < 0.05
        b += stats.levene(x, y, center="median").pvalue < 0.05
        c += perm_p(x, y, rng) < 0.05
    print(f"{lab:>12s} {a / B:8.4f} {b / B:14.4f} {c / B:16.4f}")

두 집단의 분산이 실제로 같음, n=25 씩, 명목 0.05
          분포     F 검정       브라운-포사이드       순열(로그 분산비)
          정규   0.0520         0.0387           0.0462
        t(5)   0.1658         0.0398           0.0488
          지수   0.2732         0.0460           0.0565
        로그정규   0.5020         0.0430           0.0580

z = [x - median(x)] ++ [y - median(y)]

from scipy import stats

def n_var(ratio, power=0.80, alpha=0.05):
    for n in range(4, 100_000):
        d = n - 1
        lo = stats.f.ppf(alpha / 2, d, d)
        hi = stats.f.ppf(1 - alpha / 2, d, d)
        pw = stats.f.cdf(lo / ratio, d, d) + stats.f.sf(hi / ratio, d, d)
        if pw >= power:
            return n

print(f"{'σ1²/σ2²':>8s} {'집단당 n':>9s}")
for r in [1.5, 2.0, 2.5, 3.0, 4.0]:
    print(f"{r:8.2f} {n_var(r):9d}")

 σ1²/σ2²     집단당 n
    1.50       193
    2.00        68
    2.50        40
    3.00        28
    4.00        19

두 집단의 분산을 비교하려 한다
    │
    ├─ 목적이 "분산분석을 써도 되나" 확인이라면
    │     └─→ 하지 말 것. 웰치를 쓴다
    │
    └─ 분산 자체가 연구 질문이라면
          │
          ├─ 정규성이 보장됨 ──→ F 검정 + 분산비 신뢰구간
          │
          └─ 그 외 ──→ 브라운-포사이드  또는  순열검정
                          (중앙값으로 중심을 맞춘 뒤 섞는다)

공장 A: n = 30, S² = 2.10
공장 B: n = 30, S² = 3.80
분산비 = 0.553

자료가 오른쪽으로 치우쳐 F 검정은 부적절.
브라운-포사이드 검정: W = ..., p = ...
분산비의 95% 부트스트랩 구간: [..., ...]