두 분산에 대한 F-검정¶

from scipy.stats import f

def test_ratio_two_variances(n1, s1, n2, s2, theta0=1.0,
                              alt="two-sided", alpha=0.05):
    """H0: sigma1^2 / sigma2^2 = theta0. **정규모집단**을 가정한다.

    F-검정은 정규성에서 벗어나는 데 특히 약하다. 꼬리가 조금만 두꺼워도
    제1종 오류율이 크게 부풀기 때문에, 실무에서 등분산 사전검정으로 쓰는 것은
    권장되지 않는다(그냥 Welch를 쓰는 편이 낫다).
    s1, s2에는 ddof=1로 계산한 표본표준편차를 넣는다.
    """
    # 자유도의 **순서**가 중요하다. 분자 쪽이 df1이다.
    # s1과 s2를 바꿔 넣으면 자유도도 함께 바꿔야 하며, 그러지 않으면
    # 오류 없이 조용히 틀린 p-값이 나온다.
    df1, df2 = n1 - 1, n2 - 1
    F_stat = (s1 ** 2 / s2 ** 2) / theta0
    if alt == "two-sided":
        p = 2 * min(f.cdf(F_stat, df1, df2),
                     1 - f.cdf(F_stat, df1, df2))
    elif alt == "less":
        p = f.cdf(F_stat, df1, df2)
    else:
        p = 1 - f.cdf(F_stat, df1, df2)
    return F_stat, p, (p < alpha)

import numpy as np

# (1) 1/F_{d1,d2} ~ F_{d2,d1} 을 꼬리확률로 확인한다.
d1, d2 = 14, 11
print("     x    P(F_{14,11} >= x)   P(F_{11,14} <= 1/x)")
for x in (0.5, 1.0, 2.0864197530864197, 5.0):
    print(f"{x:7.4f}        {f.sf(x, d1, d2):.12f}      {f.cdf(1 / x, d2, d1):.12f}")

# 집단을 바꿀 때 세 대립가설이 모두 짝을 맞추는가
print("\n집단을 바꾸고 자유도와 방향까지 함께 바꾼 경우")
for alt_a, alt_b in (("greater", "less"), ("less", "greater"),
                     ("two-sided", "two-sided")):
    _, pa, _ = test_ratio_two_variances(15, 1.3, 12, 0.9, alt=alt_a)
    _, pb, _ = test_ratio_two_variances(12, 0.9, 15, 1.3, alt=alt_b)
    print(f"  {alt_a:>9} / {alt_b:<9}  p = {pa:.15f}  {pb:.15f}"
          f"   차 = {abs(pa - pb):.2e}")

# 자유도를 함께 바꾸는 것을 잊으면
print("\n자유도를 바꾸지 않고 s1, s2 만 맞바꾼 경우 (틀린 계산)")
F_wrong = 0.9**2 / 1.3**2
for alt in ("less", "two-sided"):
    _, p_ok, _ = test_ratio_two_variances(12, 0.9, 15, 1.3, alt=alt)
    p_bad = (f.cdf(F_wrong, d1, d2) if alt == "less"
             else 2 * min(f.cdf(F_wrong, d1, d2), f.sf(F_wrong, d1, d2)))
    print(f"  {alt:>9}:  올바른 p = {p_ok:.6f},  틀린 p = {p_bad:.6f}"
          f"   (상대오차 {abs(p_bad - p_ok) / p_ok:.1%})")

     x    P(F_{14,11} >= x)   P(F_{11,14} <= 1/x)
 0.5000        0.888783358782      0.888783358782
 1.0000        0.509190497951      0.509190497951
 2.0864        0.112914221518      0.112914221518
 5.0000        0.005436708133      0.005436708133

집단을 바꾸고 자유도와 방향까지 함께 바꾼 경우
    greater / less       p = 0.112914221518176  0.112914221518176   차 = 4.16e-17
       less / greater    p = 0.887085778481824  0.887085778481824   차 = 0.00e+00
  two-sided / two-sided  p = 0.225828443036351  0.225828443036351   차 = 8.33e-17

자유도를 바꾸지 않고 s1, s2 만 맞바꾼 경우 (틀린 계산)
       less:  올바른 p = 0.112914,  틀린 p = 0.098065   (상대오차 13.2%)
  two-sided:  올바른 p = 0.225828,  틀린 p = 0.196131   (상대오차 13.2%)

F_stat, p, reject = test_ratio_two_variances(
    n1=15, s1=1.3, n2=12, s2=0.9, theta0=1.0, alt="greater"
)
print("F:", F_stat, "p:", p, "reject:", reject)

# 두 집단의 역할을 바꾸면 F는 역수가 되고 대립가설의 방향도 뒤집힌다.
# 제대로 바꾸면 p-값은 같아야 한다.
F2, p2, reject2 = test_ratio_two_variances(
    n1=12, s1=0.9, n2=15, s2=1.3, theta0=1.0, alt="less"
)
print("F:", F2, "p:", p2, "reject:", reject2)

F: 2.0864197530864197 p: 0.11291422151817565 reject: False
F: 0.47928994082840237 p: 0.11291422151817561 reject: False

import numpy as np
from scipy import stats

n1, n2 = 15, 12
d1, d2 = n1 - 1, n2 - 1

# 기각에 필요한 관측 분산비
up_one = stats.f.ppf(0.95, d1, d2)
hi = stats.f.ppf(0.975, d1, d2)
lo = stats.f.ppf(0.025, d1, d2)
print(f"관측된 s1/s2 = {1.3 / 0.9:.6f},  s1^2/s2^2 = {1.3**2 / 0.9**2:.6f}")
print(f"단측 5% 임계값  F_0.95(14,11) = {up_one:.6f}"
      f"  →  s1/s2 가 {np.sqrt(up_one):.6f} 이상이어야 기각")
print(f"양측 5% 채택역  ({lo:.6f}, {hi:.6f})"
      f"  →  s1/s2 로는 ({np.sqrt(lo):.6f}, {np.sqrt(hi):.6f})")

# 정규성이 깨지면 이 임계값이 무슨 뜻을 잃는가
rng = np.random.default_rng(11)
B = 200_000
gens = (
    ("정규", 3.0, lambda size: rng.standard_normal(size)),
    ("균등", 1.8, lambda size: rng.uniform(0, 1, size)),
    ("t(5)", 9.0, lambda size: rng.standard_t(5, size)),
    ("지수", 9.0, lambda size: rng.exponential(1.0, size)),
    ("로그정규", 113.9, lambda size: rng.lognormal(0.0, 1.0, size)),
)
print(f"\nH0 가 참인 자료로 재는 실제 오류율 (반복 {B:,} 회, 오차 약 ±0.001)")
print(f"{'모집단':>9} {'beta2':>7} {'양측 실제':>10} {'양측 극한':>10}"
      f" {'단측 실제':>10} {'단측 극한':>10}")
for name, b2, gen in gens:
    x1, x2 = gen((B, n1)), gen((B, n2))
    F = x1.var(1, ddof=1) / x2.var(1, ddof=1)
    r = np.sqrt(2 / (b2 - 1))
    print(f"{name:>9} {b2:7.1f} {np.mean((F < lo) | (F > hi)):10.4f} "
          f"{2 * stats.norm.sf(1.959964 * r):10.4f} "
          f"{np.mean(F > up_one):10.4f} {stats.norm.sf(1.644854 * r):10.4f}")

관측된 s1/s2 = 1.444444,  s1^2/s2^2 = 2.086420
단측 5% 임계값  F_0.95(14,11) = 2.738648  →  s1/s2 가 1.654886 이상이어야 기각
양측 5% 채택역  (0.323145, 3.358810)  →  s1/s2 로는 (0.568458, 1.832706)

H0 가 참인 자료로 재는 실제 오류율 (반복 200,000 회, 오차 약 ±0.001)
      모집단   beta2      양측 실제      양측 극한      단측 실제      단측 극한
       정규     3.0     0.0496     0.0500     0.0504     0.0500
       균등     1.8     0.0076     0.0019     0.0119     0.0047
     t(5)     9.0     0.1403     0.3271     0.1088     0.2054
       지수     9.0     0.2465     0.3271     0.1723     0.2054
     로그정규   113.9     0.4405     0.7942     0.2743     0.4134

import numpy as np
from scipy import stats

print(f"{'ν1':>4s} {'ν2':>4s} {'F(0.025)':>10s} {'F(0.975)':>10s} "
      f"{'1/F(0.975; ν2,ν1)':>19s}")
for v1, v2 in [(9, 7), (19, 24), (15, 20), (5, 5), (30, 10)]:
    lo = stats.f.ppf(0.025, v1, v2)
    hi = stats.f.ppf(0.975, v1, v2)
    check = 1 / stats.f.ppf(0.975, v2, v1)
    print(f"{v1:4d} {v2:4d} {lo:10.5f} {hi:10.5f} {check:19.5f}")

def f_test_var_ratio(s1, n1, s2, n2, alpha=0.05):
    """양측 F 검정. (F, p, 아래 임계값, 위 임계값) 을 돌려준다."""
    d1, d2 = n1 - 1, n2 - 1
    F = s1**2 / s2**2
    p = 2 * min(stats.f.cdf(F, d1, d2), stats.f.sf(F, d1, d2))
    return (F, p,
            stats.f.ppf(alpha / 2, d1, d2),
            stats.f.ppf(1 - alpha / 2, d1, d2))

F, p, lo, hi = f_test_var_ratio(4.2, 20, 3.1, 25)
print(f"\n연습문제 1: F={F:.4f}  p={p:.4f}  기각역 밖: ({lo:.4f}, {hi:.4f})")
print(f"  → {'기각' if not (lo < F < hi) else '기각하지 못함'}")

_, _, lo, hi = f_test_var_ratio(1, 10, 1, 8, alpha=0.10)
print(f"연습문제 5: α=0.10, ν=(9,7)  F_L={lo:.4f}  F_U={hi:.4f}")

  ν1   ν2   F(0.025)   F(0.975)   1/F(0.975; ν2,ν1)
   9    7    0.23826    4.82322             0.23826
  19   24    0.40778    2.34515             0.40778
  15   20    0.36286    2.57310             0.36286
   5    5    0.13993    7.14638             0.13993
  30   10    0.39822    3.31102             0.39822

연습문제 1: F=1.8356  p=0.1600  기각역 밖: (0.4078, 2.3452)
  → 기각하지 못함
연습문제 5: α=0.10, ν=(9,7)  F_L=0.3037  F_U=3.6767

import numpy as np
from scipy import stats

def ci_F(x, y, alpha=0.05):
    """고전적 F 구간."""
    n1, n2 = len(x), len(y)
    f = x.var(ddof=1) / y.var(ddof=1)
    return (f / stats.f.ppf(1 - alpha / 2, n1 - 1, n2 - 1),
            f / stats.f.ppf(alpha / 2, n1 - 1, n2 - 1))

def ci_log_normal(x, y, alpha=0.05):
    """정규성을 가정한 로그 분산비 구간 (κ=3 고정)."""
    n1, n2 = len(x), len(y)
    c = np.log(x.var(ddof=1) / y.var(ddof=1))
    se = np.sqrt(2 / (n1 - 1) + 2 / (n2 - 1))
    z = stats.norm.ppf(1 - alpha / 2)
    return np.exp(c - z * se), np.exp(c + z * se)

def ci_kurtosis(x, y, alpha=0.05):
    """첨도를 자료에서 추정해 보정한 구간."""
    def kur(a):
        m = a.mean()
        return ((a - m)**4).mean() / (((a - m)**2).mean())**2

    n1, n2 = len(x), len(y)
    c = np.log(x.var(ddof=1) / y.var(ddof=1))
    se = np.sqrt((kur(x) - (n1 - 3) / (n1 - 1)) / n1
                 + (kur(y) - (n2 - 3) / (n2 - 1)) / n2)
    z = stats.norm.ppf(1 - alpha / 2)
    return np.exp(c - z * se), np.exp(c + z * se)

rng = np.random.default_rng(99)
M, n = 10_000, 30
gens = {"정규":     lambda k: rng.standard_normal(k),
        "균등":     lambda k: rng.uniform(-1, 1, k),
        "t(5)":     lambda k: rng.standard_t(5, k),
        "지수":     lambda k: rng.exponential(1, k),
        "로그정규": lambda k: rng.lognormal(0, 1, k)}

print(f"{'분포':>8s} {'F 구간':>9s} {'로그(κ=3)':>11s} {'첨도 보정':>11s}")
for name, gen in gens.items():
    a = b = c = 0
    for _ in range(M):
        x, y = gen(n), gen(n)            # 참 분산비 = 1
        lo, hi = ci_F(x, y);          a += lo <= 1 <= hi
        lo, hi = ci_log_normal(x, y); b += lo <= 1 <= hi
        lo, hi = ci_kurtosis(x, y);   c += lo <= 1 <= hi
    print(f"{name:>8s} {a / M:9.4f} {b / M:11.4f} {c / M:11.4f}")

      분포      F 구간     로그(κ=3)       첨도 보정
      정규    0.9469      0.9419      0.9205
      균등    0.9973      0.9963      0.9512
    t(5)    0.8319      0.8232      0.8912
      지수    0.7260      0.7163      0.8508
    로그정규    0.4742      0.4665      0.7415

rng = np.random.default_rng(1234)

def ci_boot(x, y, B, alpha, rng):
    """로그 분산비의 백분위 붓스트랩 구간."""
    n1, n2 = len(x), len(y)
    bs = np.empty(B)
    for b in range(B):
        xb = x[rng.integers(0, n1, n1)]
        yb = y[rng.integers(0, n2, n2)]
        bs[b] = np.log(xb.var(ddof=1) / yb.var(ddof=1))
    lo, hi = np.percentile(bs, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return np.exp(lo), np.exp(hi)

M, n = 2_000, 30
for name, gen in [("정규",     lambda k: rng.standard_normal(k)),
                  ("t(5)",     lambda k: rng.standard_t(5, k)),
                  ("로그정규", lambda k: rng.lognormal(0, 1, k))]:
    c = 0
    for _ in range(M):
        x, y = gen(n), gen(n)
        lo, hi = ci_boot(x, y, 399, 0.05, rng)
        c += lo <= 1 <= hi
    print(f"{name:>8s} 붓스트랩 포함률 {c / M:.4f}")

      정규 붓스트랩 포함률 0.9250
    t(5) 붓스트랩 포함률 0.9175
    로그정규 붓스트랩 포함률 0.8860

import numpy as np
from scipy import stats
from itertools import combinations

rng = np.random.default_rng(2718)
M, n = 10_000, 20

print(f"{'집단 수':>7s} {'쌍':>4s} {'쌍별 F(무보정)':>15s} {'본페로니':>10s} "
      f"{'바틀렛':>8s} {'B-F':>8s}")
for k in [2, 3, 5, 8]:
    n_pair = k * (k - 1) // 2
    a = b = c = d = 0
    for _ in range(M):
        g = [rng.standard_normal(n) for _ in range(k)]   # 모든 분산이 같다
        ps = []
        for i, j in combinations(range(k), 2):
            F = g[i].var(ddof=1) / g[j].var(ddof=1)
            ps.append(2 * min(stats.f.cdf(F, n - 1, n - 1),
                              stats.f.sf(F, n - 1, n - 1)))
        ps = np.array(ps)
        a += (ps < 0.05).any()
        b += (ps < 0.05 / n_pair).any()
        c += stats.bartlett(*g).pvalue < 0.05
        d += stats.levene(*g, center='median').pvalue < 0.05
    print(f"{k:7d} {n_pair:4d} {a / M:15.4f} {b / M:10.4f} "
          f"{c / M:8.4f} {d / M:8.4f}")

   집단 수    쌍       쌍별 F(무보정)       본페로니      바틀렛      B-F
      2    1          0.0489     0.0489   0.0489   0.0390
      3    3          0.1237     0.0449   0.0538   0.0392
      5   10          0.2795     0.0390   0.0498   0.0325
      8   28          0.4994     0.0310   0.0498   0.0284

import numpy as np
from scipy import stats

LSL, USL = 9.5, 10.5
processes = [("기존 공정", 10.00, 0.12, 30),
             ("신 공정", 10.02, 0.09, 30)]

for name, mu, s, n in processes:
    cp = (USL - LSL) / (6 * s)
    cpk = min(USL - mu, mu - LSL) / (3 * s)
    ppm = (stats.norm.cdf(LSL, mu, s) + stats.norm.sf(USL, mu, s)) * 1e6
    print(f"{name}: 평균 {mu:.3f}  SD {s:.3f}  "
          f"Cp={cp:.4f}  Cpk={cpk:.4f}  불량 {ppm:.1f} ppm")

s1, n1 = 0.12, 30
s2, n2 = 0.09, 30
d1, d2 = n1 - 1, n2 - 1
F = s1**2 / s2**2
p = 2 * min(stats.f.cdf(F, d1, d2), stats.f.sf(F, d1, d2))
lo = F / stats.f.ppf(0.975, d1, d2)
hi = F / stats.f.ppf(0.025, d1, d2)
print(f"\n분산비 F = {F:.4f},  p = {p:.4f}")
print(f"  σ1²/σ2² 의 95% CI ({lo:.4f}, {hi:.4f})")
print(f"  σ1/σ2   의 95% CI ({np.sqrt(lo):.4f}, {np.sqrt(hi):.4f})")
print(f"  Cp 비(기존/신) = {1 / np.sqrt(F):.4f}, "
      f"95% CI ({1 / np.sqrt(hi):.4f}, {1 / np.sqrt(lo):.4f})")

기존 공정: 평균 10.000  SD 0.120  Cp=1.3889  Cpk=1.3889  불량 30.9 ppm
신 공정: 평균 10.020  SD 0.090  Cp=1.8519  Cpk=1.7778  불량 0.1 ppm

분산비 F = 1.7778,  p = 0.1271
  σ1²/σ2² 의 95% CI (0.8462, 3.7351)
  σ1/σ2   의 95% CI (0.9199, 1.9326)
  Cp 비(기존/신) = 0.7500, 95% CI (0.5174, 1.0871)