Levene 검정¶

import inspect

import numpy as np
import scipy.stats as stats
from scipy.stats import levene

# 기본값이 무엇인지 함수 서명에서 직접 확인한다.
print("levene 의 center 기본값 =", inspect.signature(levene).parameters["center"].default)

# 표준편차가 1.0과 1.5로 다른 두 집단
group1 = stats.norm(0, 1.0).rvs(50, random_state=0)
group2 = stats.norm(0, 1.5).rvs(50, random_state=1)

stat, pval = levene(group1, group2, center='median')
print(f"F = {stat:.4f}, p = {pval:.4f}")

# 세 중심을 나란히 둔다. 'mean' 이 원래의 Levene(1960), 'median' 이 Brown-Forsythe.
print(f"\n{'center':>10}{'W':>10}{'p':>9}{'Zbar_1':>9}{'Zbar_2':>9}")
for center in ("mean", "median", "trimmed"):
    kw = {"proportiontocut": 0.1} if center == "trimmed" else {}
    W, p = levene(group1, group2, center=center, **kw)
    if center == "mean":
        locs = [g.mean() for g in (group1, group2)]
    elif center == "median":
        locs = [np.median(g) for g in (group1, group2)]
    else:
        locs = [stats.trim_mean(g, 0.1) for g in (group1, group2)]
    zbar = [np.abs(g - c).mean() for g, c in zip((group1, group2), locs)]
    print(f"{center:>10}{W:>10.4f}{p:>9.4f}{zbar[0]:>9.4f}{zbar[1]:>9.4f}")

# (1) 의 부등식: 중앙값이 sum |y - c| 를 최소화하므로 Zbar 는 중앙값에서 가장 작다.
print("\n중심을 바꿔 가며 group2 의 Zbar 를 재 본다 (중앙값에서 최소여야 한다)")
med2 = np.median(group2)
for label, c in [("mean", group2.mean()), ("median", med2),
                 ("median-0.3", med2 - 0.3), ("median+0.3", med2 + 0.3)]:
    print(f"  {label:>12}: c = {c:+.4f}   Zbar = {np.abs(group2 - c).mean():.6f}")

# 치우친 자료에서는 두 중심이 크게 갈린다. 로그정규 세 집단, 모분산은 모두 같다.
rng = np.random.default_rng(3)
g = [np.exp(rng.normal(size=40)) for _ in range(3)]
print("\n로그정규 세 집단 (모분산 동일)")
for center in ("mean", "median"):
    W, p = levene(*g, center=center)
    print(f"  center={center:>6}:  W = {W:.4f}   p = {p:.4f}")

levene 의 center 기본값 = median
F = 2.8007, p = 0.0974

    center         W        p   Zbar_1   Zbar_2
      mean    3.3525   0.0701   0.9017   1.1842
    median    2.8007   0.0974   0.9015   1.1695
   trimmed    3.7073   0.0578   0.9016   1.1800

중심을 바꿔 가며 group2 의 Zbar 를 재 본다 (중앙값에서 최소여야 한다)
          mean: c = -0.0383   Zbar = 1.184206
        median: c = -0.2732   Zbar = 1.169549
    median-0.3: c = -0.5732   Zbar = 1.199370
    median+0.3: c = +0.0268   Zbar = 1.191244

로그정규 세 집단 (모분산 동일)
  center=  mean:  W = 2.0636   p = 0.1316
  center=median:  W = 0.6888   p = 0.5042

import numpy as np
import scipy.stats as stats

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

sigma_y=1.00: F=0.00, p=1.000
sigma_y=1.05: F=0.20, p=0.652
sigma_y=1.10: F=0.78, p=0.379
sigma_y=1.15: F=1.67, p=0.198
sigma_y=1.20: F=2.82, p=0.095

import scipy.optimize as opt

a = np.abs(x - np.median(x))
C = (size - 1) * size * a.mean() ** 2 / ((a - a.mean()) ** 2).sum()
print(f"C = {C:.5f}")

print(f"\n{'sigma_y':>8}{'levene':>12}{'공식':>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)
    W, p = stats.levene(x, y)
    W_formula = C * (scale - 1) ** 2 / (1 + scale ** 2)
    print(f"{scale:>8.2f}{W:>12.6f}{W_formula:>12.6f}{p:>9.4f}")

# p = 0.05 가 되는 sigma 를 공식을 거꾸로 풀어 찾는다. 격자 탐색이 아니다.
crit = stats.f(1, 2 * size - 2).ppf(0.95)
root = opt.brentq(lambda s: C * (s - 1) ** 2 / (1 + s ** 2) - crit, 1.0001, 5.0)
print(f"\n임계값 F_0.95(1,198) = {crit:.4f}")
print(f"W = 임계값이 되는 sigma_y = {root:.4f}")
print(f"sigma_y -> inf 에서 W 의 상한 = C = {C:.4f}  (해가 존재할 조건)")

C = 171.94090

 sigma_y      levene          공식        p
    1.00    0.000000    0.000000   1.0000
    1.05    0.204448    0.204448   0.6516
    1.10    0.778013    0.778013   0.3788
    1.15    1.665735    1.665735   0.1983
    1.20    2.818703    2.818703   0.0947

임계값 F_0.95(1,198) = 3.8889
W = 임계값이 되는 sigma_y = 1.2395
sigma_y -> inf 에서 W 의 상한 = C = 171.9409  (해가 존재할 조건)

import numpy as np
from scipy import stats

def levene_manual(groups, center="median"):
    """레빈 검정 = |관측 - 중심| 에 대한 일원배치 분산분석."""
    c = np.median if center == "median" else np.mean
    Z = [np.abs(g - c(g)) for g in groups]
    return stats.f_oneway(*Z)

rng = np.random.default_rng(1212)
gs = [rng.normal(0, 1, 12), rng.normal(0, 1.8, 15), rng.normal(0, 2.5, 10)]
for c in ["mean", "median"]:
    m = levene_manual(gs, c)
    s = stats.levene(*gs, center=c)
    print(f"  center={c:7s}  자작 F = {m.statistic:.8f}, p = {m.pvalue:.8f}")
    print(f"  {'':15s}  scipy F = {s.statistic:.8f}, p = {s.pvalue:.8f}")

  center=mean     자작 F = 3.13477119, p = 0.05630381
                   scipy F = 3.13477119, p = 0.05630381
  center=median   자작 F = 3.05861559, p = 0.06005033
                   scipy F = 3.05861559, p = 0.06005033

import numpy as np
from scipy import stats

rng = np.random.default_rng(1212)
B = 8_000
print("분산이 실제로 같음, k=3, n=20, 명목 0.05")
print(f"{'분포':>20s} {'mean':>8s} {'median':>8s} {'trimmed(0.1)':>13s}")
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)),
               "균등": lambda n: rng.uniform(-1, 1, n)}.items():
    a = b = c = 0
    for _ in range(B):
        gs = [f(20) for _ in range(3)]
        a += stats.levene(*gs, center="mean").pvalue < 0.05
        b += stats.levene(*gs, center="median").pvalue < 0.05
        c += stats.levene(*gs, center="trimmed",
                          proportiontocut=0.1).pvalue < 0.05
    print(f"{lab:>20s} {a / B:8.4f} {b / B:8.4f} {c / B:13.4f}")

a = b = c = 0
for _ in range(B):
    gs = [rng.normal(0, 1, 20), rng.normal(0, 1, 20), rng.normal(0, 2, 20)]
    a += stats.levene(*gs, center="mean").pvalue < 0.05
    b += stats.levene(*gs, center="median").pvalue < 0.05
    c += stats.levene(*gs, center="trimmed", proportiontocut=0.1).pvalue < 0.05
print(f"{'검정력 정규 σ=(1,1,2)':>20s} {a / B:8.4f} {b / B:8.4f} {c / B:13.4f}")

분산이 실제로 같음, k=3, n=20, 명목 0.05
                  분포     mean   median  trimmed(0.1)
                  정규   0.0517   0.0319        0.1424
                t(5)   0.0566   0.0362        0.1646
                  지수   0.1850   0.0466        0.3031
                로그정규   0.2409   0.0341        0.3713
                  균등   0.0553   0.0267        0.1356
    검정력 정규 σ=(1,1,2)   0.8125   0.7649        0.7814

center='trimmed' 일 때 scipy 는
  1. 표본 자체를 양끝에서 proportiontocut 만큼 잘라 내고
  2. 남은 자료의 평균을 중심으로 쓴다

import numpy as np
from scipy import stats

rng = np.random.default_rng(3434)
B = 6_000
print("분산이 실제로 같음, 명목 0.05")
for lab, ns in [("(15,15,15)", [15, 15, 15]),
                ("(15,15,200)", [15, 15, 200])]:
    for dl, f in [("정규", lambda n: rng.normal(0, 1, n)),
                  ("t(5)", lambda n: rng.standard_t(5, n))]:
        a = b = 0
        for _ in range(B):
            gs = [f(n) for n in ns]
            a += stats.bartlett(*gs).pvalue < 0.05
            b += stats.levene(*gs, center="median").pvalue < 0.05
        print(f"  n={lab:12s} {dl:5s}  바틀렛 {a / B:.4f}   "
              f"브라운-포사이드 {b / B:.4f}")

분산이 실제로 같음, 명목 0.05
  n=(15,15,15)   정규     바틀렛 0.0465   브라운-포사이드 0.0288
  n=(15,15,15)   t(5)   바틀렛 0.2018   브라운-포사이드 0.0333
  n=(15,15,200)  정규     바틀렛 0.0555   브라운-포사이드 0.0487
  n=(15,15,200)  t(5)   바틀렛 0.2075   브라운-포사이드 0.0430

import numpy as np
from scipy import stats

rng = np.random.default_rng(5656)
B = 6_000
print("정규, k=3, σ=(1, 1, r), 브라운-포사이드의 검정력")
print(f"{'r':>5s} " + " ".join(f"{'n=' + str(n):>8s}" for n in [10, 20, 50, 100]))
for r in [1.5, 2.0, 3.0]:
    row = []
    for n in [10, 20, 50, 100]:
        a = sum(stats.levene(rng.normal(0, 1, n), rng.normal(0, 1, n),
                             rng.normal(0, r, n)).pvalue < 0.05
                for _ in range(B))
        row.append(a / B)
    print(f"{r:5.1f} " + " ".join(f"{x:8.4f}" for x in row))

정규, k=3, σ=(1, 1, r), 브라운-포사이드의 검정력
    r     n=10     n=20     n=50    n=100
  1.5   0.1290   0.3255   0.7697   0.9763
  2.0   0.3700   0.7703   0.9965   1.0000
  3.0   0.7462   0.9897   1.0000   1.0000

집단별 표준편차: 1.17, 1.09, 2.31  (비 2.1배)
브라운-포사이드 검정: W = 5.42, p = 0.007
→ 등분산을 가정하지 않고 웰치 분산분석 + 게임스-하웰 사후검정