분산분석 F-통계량 모의실험¶

import numpy as np
from scipy import stats

np.random.seed(42)      # 아래 출력을 재현하려면 고정한다

def simulate_f(mu, sigma, sizes, n_sim=1000):
    """지정한 평균/표준편차/크기로 집단을 만들고 F-검정을 n_sim번 반복한다."""
    F_vals, p_vals = [], []
    for _ in range(n_sim):
        # 집단마다 평균도 표준편차도 크기도 다를 수 있게 zip으로 묶는다.
        groups = [stats.norm.rvs(m, s, n)
                  for m, s, n in zip(mu, sigma, sizes)]
        F, p = stats.f_oneway(*groups)
        F_vals.append(F)
        p_vals.append(p)
    return np.array(F_vals), np.array(p_vals)


# (1) 의 불변성: 자료를 a + c*y 로 바꾸어도 F 가 그대로인가.
rng = np.random.default_rng(7)
y = [rng.normal(3, 6, m) for m in (10, 20, 30)]
F0, _ = stats.f_oneway(*y)
F1, _ = stats.f_oneway(*[100 + 7 * g for g in y])
print(f"원자료      F = {F0!r}")
print(f"100 + 7y 로 F = {F1!r}")

# sigma 를 100배 키워도 귀무분포가 같은가.
for sigma in (6, 600):
    F_vals, p_vals = simulate_f((3, 3, 3), (sigma,) * 3, (10, 20, 30), n_sim=10000)
    ks_F = stats.kstest(F_vals, stats.f(2, 57).cdf)
    ks_p = stats.kstest(p_vals, 'uniform')
    print(f"\nsigma = {sigma}")
    print(f"  F 의 평균 {F_vals.mean():.4f} (이론 {57 / 55:.4f}),  "
          f"표준편차 {F_vals.std(ddof=1):.4f} (이론 {np.sqrt(stats.f(2, 57).var()):.4f})")
    print(f"  F ~ F(2,57) 인가:  KS = {ks_F.statistic:.4f}, p = {ks_F.pvalue:.4f}")
    print(f"  p-값이 균등인가:   KS = {ks_p.statistic:.4f}, p = {ks_p.pvalue:.4f}")
    print(f"  기각률 = {np.mean(p_vals < 0.05):.4f}")

np.random.seed(42)   # 보기 2 의 출력을 그대로 재현하려고 씨앗을 되돌린다

원자료      F = 2.0897728316088156
100 + 7y 로 F = 2.0897728316088156

sigma = 6
  F 의 평균 1.0331 (이론 1.0364),  표준편차 1.0693 (이론 1.0748)
  F ~ F(2,57) 인가:  KS = 0.0055, p = 0.9212
  p-값이 균등인가:   KS = 0.0055, p = 0.9212
  기각률 = 0.0496

sigma = 600
  F 의 평균 1.0389 (이론 1.0364),  표준편차 1.0832 (이론 1.0748)
  F ~ F(2,57) 인가:  KS = 0.0083, p = 0.4886
  p-값이 균등인가:   KS = 0.0083, p = 0.4886
  기각률 = 0.0501

scenarios = [
    ("큰 분리, 등분산",        (3, 6, 9),       (6, 6, 6),      (10, 20, 30)),
    ("아주 작은 분리, 등분산", (3, 3.1, 2.9),   (6, 6, 6),      (10, 20, 30)),
    ("중간 분리, 큰 분산, n=10",  (3, 5, 6),    (10, 10, 10),   (10, 10, 10)),
    ("중간 분리, 큰 분산, n=5000",(3, 5, 6),    (10, 10, 10),   (5000, 5000, 5000)),
    ("분리 없음, 등분산",      (3, 3, 3),       (1, 1, 1),      (10, 20, 30)),
]

print(f"{'시나리오':<26} {'F_crit':>7} {'기각률':>8}")
for name, mu, sigma, sizes in scenarios:
    F_vals, p_vals = simulate_f(mu, sigma, sizes, n_sim=1000)
    df1 = len(mu) - 1                # k - 1
    df2 = sum(sizes) - len(mu)       # N - k
    # 임계값은 자유도만으로 정해진다. 자료와 무관하다.
    F_crit = stats.f.ppf(0.95, df1, df2)
    reject_pct = np.mean(F_vals > F_crit) * 100
    print(f"{name:<26} {F_crit:>7.3f} {reject_pct:>7.1f}%")

시나리오                        F_crit      기각률
큰 분리, 등분산                    3.159    68.1%
아주 작은 분리, 등분산                3.159     4.9%
중간 분리, 큰 분산, n=10            3.354     7.1%
중간 분리, 큰 분산, n=5000          2.996   100.0%
분리 없음, 등분산                   3.159     4.8%

import warnings

np.random.seed(2024)
print(f"{'scenario':<28}{'lambda':>10}{'theory':>9}{'sim':>8}{'MCse':>8}{'z':>7}")
with warnings.catch_warnings():            # 거대한 lambda 에서 나는 overflow 경고를 끈다
    warnings.simplefilter("ignore")
    for name, mu, sigma, sizes in scenarios:
        mu_a, n_a = np.array(mu, float), np.array(sizes, float)
        N, k = n_a.sum(), len(mu)
        mu_bar = (n_a * mu_a).sum() / N       # 표본크기로 가중한 평균
        lam = (n_a * (mu_a - mu_bar) ** 2).sum() / sigma[0] ** 2
        d1, d2 = k - 1, N - k
        crit = stats.f(d1, d2).ppf(0.95)
        power = 0.05 if lam == 0 else float(stats.ncf(d1, d2, lam).sf(crit))
        B = 300 if max(sizes) > 1000 else 10000
        F_vals, _ = simulate_f(mu, sigma, sizes, n_sim=B)
        sim = np.mean(F_vals > crit)
        se = np.sqrt(power * (1 - power) / B)
        z = (sim - power) / se if se > 0 else 0.0
        print(f"{name:<28}{lam:>10.4f}{power:>9.4f}{sim:>8.4f}{se:>8.4f}{z:>7.2f}")

scenario                        lambda   theory     sim    MCse      z
큰 분리, 등분산                       8.3333   0.7120  0.7059  0.0045  -1.34
아주 작은 분리, 등분산                   0.0134   0.0510  0.0475  0.0022  -1.57
중간 분리, 큰 분산, n=10               0.4667   0.0829  0.0839  0.0028   0.38
중간 분리, 큰 분산, n=5000           233.3333   1.0000  1.0000  0.0000   0.00
분리 없음, 등분산                      0.0000   0.0500  0.0539  0.0022   1.79

for n in [10, 20, 30, 50, 100, 200, ...]:
    set mu = [0, delta, 0] = [0, 2, 0]
    set sigma = [5, 5, 5]
    set sizes = [n, n, n]
    run simulate_f(mu, sigma, sizes, n_sim=5000)
    compute rejection_rate = mean(p_vals < 0.05)
    if rejection_rate >= 0.80:
        return n

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)  # ncf 의 내부 경고

import numpy as np
from scipy import stats

rng = np.random.default_rng(2025)

def power_theory(mu, sigma, n, alpha=0.05):
    """집단당 n 명일 때 비중심 F 로 구한 검정력."""
    mu = np.asarray(mu, float)
    k = len(mu)
    N = k * n
    lam = n * ((mu - mu.mean())**2).sum() / sigma**2
    crit = stats.f.ppf(1 - alpha, k - 1, N - k)
    return stats.ncf.sf(crit, k - 1, N - k, lam), lam

def power_sim(mu, sigma, n, M=3_000, alpha=0.05):
    rej = 0
    for _ in range(M):
        groups = [rng.normal(m, sigma, n) for m in mu]
        rej += stats.f_oneway(*groups).pvalue < alpha
    return rej / M

mu, sigma = [0, 2, 0], 5.0
print(f"{'n':>5s} {'λ':>9s} {'이론 검정력':>12s} {'모의 검정력':>12s}")
for n in [10, 20, 30, 50, 75, 100, 150]:
    th, lam = power_theory(mu, sigma, n)
    print(f"{n:5d} {lam:9.4f} {th:12.4f} {power_sim(mu, sigma, n):12.4f}")

for n in range(5, 500):
    if power_theory(mu, sigma, n)[0] >= 0.80:
        print(f"\n검정력 80% 에 필요한 군당 n = {n}")
        break

    n         λ       이론 검정력       모의 검정력
   10    1.0667       0.1288       0.1237
   20    2.1333       0.2279       0.2247
   30    3.2000       0.3303       0.3320
   50    5.3333       0.5225       0.5260
   75    8.0000       0.7113       0.7130
  100   10.6667       0.8371       0.8367
  150   16.0000       0.9555       0.9543

검정력 80% 에 필요한 군당 n = 92

import numpy as np
from scipy import stats

rng = np.random.default_rng(707)

def welch_anova(*groups):
    """웰치의 일원배치 분산분석. (F, p, df2) 를 돌려준다."""
    k = len(groups)
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    w = n / v
    W = w.sum()
    m_tilde = (w * m).sum() / W
    A = (w * (m - m_tilde)**2).sum() / (k - 1)
    tmp = np.sum((1 - w / W)**2 / (n - 1))
    B = 1 + 2 * (k - 2) / (k * k - 1) * tmp
    F = A / B
    df2 = (k * k - 1) / (3 * tmp)
    return F, stats.f.sf(F, k - 1, df2), df2

M = 10_000
print(f"{'σ':>12s} {'n':>14s} {'고전 F':>9s} {'Welch F':>9s} "
      f"{'알렉산더·고번':>14s}")
for sig, ns in [((1, 1, 1), (50, 50, 5)),
                ((1, 1, 10), (50, 50, 5)),
                ((1, 1, 10), (5, 50, 50)),
                ((3, 2, 1), (10, 20, 30)),
                ((1, 1, 4), (20, 20, 20))]:
    a = b = c = 0
    for _ in range(M):
        g = [rng.normal(0, s, n) for s, n in zip(sig, ns)]   # 평균은 모두 0
        a += stats.f_oneway(*g).pvalue < 0.05
        b += welch_anova(*g)[1] < 0.05
        c += stats.alexandergovern(*g).pvalue < 0.05
    print(f"{str(sig):>12s} {str(ns):>14s} {a / M:9.4f} {b / M:9.4f} "
          f"{c / M:14.4f}")

           σ              n      고전 F   Welch F        알렉산더·고번
   (1, 1, 1)    (50, 50, 5)    0.0508    0.0585         0.0598
  (1, 1, 10)    (50, 50, 5)    0.5971    0.0568         0.0469
  (1, 1, 10)    (5, 50, 50)    0.0215    0.0558         0.0584
   (3, 2, 1)   (10, 20, 30)    0.1641    0.0559         0.0553
   (1, 1, 4)   (20, 20, 20)    0.0783    0.0500         0.0485

import numpy as np
from scipy import stats

rng = np.random.default_rng(303)
M = 5_000

print(f"{'상황':>22s} {'KS D':>8s} {'KS p':>10s} {'기각률':>8s} "
      f"{'p<0.2 비율':>11s}")
setups = [("등분산 (10,20,30)", (6, 6, 6), (10, 20, 30)),
          ("이분산 큰σ↔작은n", (1, 1, 10), (50, 50, 5)),
          ("이분산 큰σ↔큰n", (1, 1, 10), (5, 50, 50)),
          ("등분산 (5,5,5)", (1, 1, 1), (5, 5, 5))]
for label, sig, ns in setups:
    pvals = np.empty(M)
    for i in range(M):
        g = [rng.normal(0, s, n) for s, n in zip(sig, ns)]
        pvals[i] = stats.f_oneway(*g).pvalue
    ks = stats.kstest(pvals, 'uniform')
    print(f"{label:>22s} {ks.statistic:8.4f} {ks.pvalue:10.4g} "
          f"{np.mean(pvals < 0.05):8.4f} {np.mean(pvals < 0.2):11.4f}")

                    상황     KS D       KS p      기각률    p<0.2 비율
        등분산 (10,20,30)   0.0083     0.8768   0.0484      0.2000
            이분산 큰σ↔작은n   0.5584          0   0.6024      0.7174
             이분산 큰σ↔큰n   0.2790          0   0.0240      0.0976
           등분산 (5,5,5)   0.0184    0.06713   0.0590      0.2162

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

import numpy as np
from scipy import stats

sigma, n, k = 5.0, 30, 4
crit = stats.f.ppf(0.95, k - 1, k * n - k)

print("① Σ(μ-μ̄)² 를 4 로 고정")
print(f"{'배치':>18s} {'중심화 평균':>30s} {'λ':>8s} {'검정력':>9s}")
for label, raw in [("하나만 높음", [3., 0, 0, 0]),
                   ("둘씩 갈림", [1., 1, -1, -1]),
                   ("고르게 퍼짐", [-1.5, -0.5, 0.5, 1.5]),
                   ("하나만 낮음", [0., 0, 0, -3.])]:
    mu = np.array(raw) - np.mean(raw)
    mu = mu / np.sqrt((mu**2).sum()) * 2.0        # Σ(μ-μ̄)² = 4
    lam = n * (mu**2).sum() / sigma**2
    print(f"{label:>18s} {str(np.round(mu, 3).tolist()):>30s} "
          f"{lam:8.4f} {stats.ncf.sf(crit, k - 1, k * n - k, lam):9.4f}")

① Σ(μ-μ̄)² 를 4 로 고정
                배치                         중심화 평균        λ       검정력
            하나만 높음 [1.732, -0.577, -0.577, -0.577]   4.8000    0.4115
             둘씩 갈림         [1.0, 1.0, -1.0, -1.0]   4.8000    0.4115
            고르게 퍼짐 [-1.342, -0.447, 0.447, 1.342]   4.8000    0.4115
            하나만 낮음  [0.577, 0.577, 0.577, -1.732]   4.8000    0.4115

rng = np.random.default_rng(808)
sigma, n, k, R = 5.0, 10, 4, 5.0        # 범위(최대-최소) = 5 로 고정
crit = stats.f.ppf(0.95, k - 1, k * n - k)

print("② 범위(최대-최소)를 5 로 고정")
print(f"{'배치':>18s} {'중심화 평균':>28s} {'Σ(μ-μ̄)²':>10s} {'λ':>7s} "
      f"{'이론':>8s} {'모의':>8s}")
configs = {"하나만 높음": [R, 0., 0, 0],
           "둘씩 갈림": [R, R, 0., 0],
           "고르게 퍼짐": list(np.linspace(0, R, 4)),
           "가운데 둘이 같음": [0., R / 2, R / 2, R]}
for label, raw in configs.items():
    mu = np.array(raw, float)
    mu = mu - mu.mean()
    ss = (mu**2).sum()
    lam = n * ss / sigma**2
    th = stats.ncf.sf(crit, k - 1, k * n - k, lam)
    rej, M = 0, 5_000
    for _ in range(M):
        g = [rng.normal(m, sigma, n) for m in mu]
        rej += stats.f_oneway(*g).pvalue < 0.05
    print(f"{label:>18s} {str(np.round(mu, 2).tolist()):>28s} {ss:10.2f} "
          f"{lam:7.3f} {th:8.4f} {rej / M:8.4f}")

② 범위(최대-최소)를 5 로 고정
                배치                       중심화 평균   Σ(μ-μ̄)²       λ       이론       모의
            하나만 높음  [3.75, -1.25, -1.25, -1.25]      18.75   7.500   0.5711   0.5584
             둘씩 갈림       [2.5, 2.5, -2.5, -2.5]      25.00  10.000   0.7097   0.7056
            고르게 퍼짐     [-2.5, -0.83, 0.83, 2.5]      13.89   5.556   0.4396   0.4492
         가운데 둘이 같음        [-2.5, 0.0, 0.0, 2.5]      12.50   5.000   0.3992   0.4004