Levene 검정 정규 대 치우친 분포 모의실험¶

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

rng = np.random.default_rng(0)


def simulate_once(n=20, dist="normal"):
    """주어진 분포에서 세 집단을 만들어 검정하고 p-값을 돌려준다.

    셋 다 같은 모수를 쓰므로 분산은 참으로 같다. 그러니 기각 비율이
    0.05 근처로 나와야 옳다.
    """
    if dist == "normal":
        g1 = rng.normal(0, 1.0, size=n)
        g2 = rng.normal(0, 1.0, size=n)
        g3 = rng.normal(0, 1.0, size=n)
    else:
        g1 = rng.lognormal(0, 1.0, size=n)
        g2 = rng.lognormal(0, 1.0, size=n)
        g3 = rng.lognormal(0, 1.0, size=n)
    _, p = levene(g1, g2, g3, center='median')
    return p


# 정규와 로그정규에서 각각 돌려 본다. Bartlett 이라면 로그정규에서
# 오류율이 크게 부풀지만, 중앙값 중심 방법은 0.05 근처를 지킨다.
alpha = 0.05
n_sims = 5000

rate = {}
for dist in ["normal", "lognormal"]:
    pvals = [simulate_once(20, dist) for _ in range(n_sims)]
    type1 = np.mean(np.array(pvals) < alpha)
    rate[dist] = type1
    print(f"Type I error (median-centered) under {dist}: {type1:.4f}")

# (1) 의 두 자로 잰다. 명목값과의 거리, 그리고 두 수치 사이의 거리.
print(f"\n몬테카를로 표준오차 (p = 0.05, B = {n_sims}): "
      f"{np.sqrt(0.05 * 0.95 / n_sims):.5f}")
for dist in ["normal", "lognormal"]:
    p = rate[dist]
    se = np.sqrt(p * (1 - p) / n_sims)
    print(f"  {dist:>10}: {p:.4f}  SE {se:.5f}  "
          f"->  0.05 로부터 {(p - 0.05) / se:+.2f} SE")
d = rate["lognormal"] - rate["normal"]
se_d = np.sqrt(sum(rate[k] * (1 - rate[k]) / n_sims for k in rate))
print(f"  두 수치의 차 {d:+.4f},  그 SE {se_d:.5f}  ->  {d / se_d:+.2f} SE")

# 모집단을 다섯으로 늘려 Bartlett 과 나란히 둔다. 같은 자료, 같은 반복수.
rng2 = np.random.default_rng(7)
n = 20
pops = [
    ("정규 N(0,1)", lambda m: rng2.normal(0, 1, (m, n)), 3.0),
    ("균등 U(0,1)", lambda m: rng2.uniform(0, 1, (m, n)), 1.8),
    ("지수 Exp(1)", lambda m: rng2.exponential(1, (m, n)), 9.0),
    ("t(5)", lambda m: rng2.standard_t(5, (m, n)), 9.0),
    ("로그정규 LN(0,1)", lambda m: rng2.lognormal(0, 1, (m, n)), 113.94),
]

print(f"\n{'모집단':>18}{'beta2':>9}{'BF 크기':>10}{'Bartlett':>10}")
bf, ba = [], []
for name, draw, b2 in pops:
    A, B, C = draw(n_sims), draw(n_sims), draw(n_sims)
    r_bf = r_ba = 0
    for i in range(n_sims):
        g = (A[i], B[i], C[i])
        if stats.levene(*g, center="median")[1] < alpha:
            r_bf += 1
        if stats.bartlett(*g)[1] < alpha:
            r_ba += 1
    bf.append(r_bf / n_sims)
    ba.append(r_ba / n_sims)
    print(f"{name:>18}{b2:>9.1f}{r_bf / n_sims:>10.4f}{r_ba / n_sims:>10.4f}")
print(f"{'폭 (최대 - 최소)':>27}{max(bf) - min(bf):>10.4f}{max(ba) - min(ba):>10.4f}")

Type I error (median-centered) under normal: 0.0370
Type I error (median-centered) under lognormal: 0.0374

몬테카를로 표준오차 (p = 0.05, B = 5000): 0.00308
      normal: 0.0370  SE 0.00267  ->  0.05 로부터 -4.87 SE
   lognormal: 0.0374  SE 0.00268  ->  0.05 로부터 -4.70 SE
  두 수치의 차 +0.0004,  그 SE 0.00379  ->  +0.11 SE

               모집단    beta2     BF 크기  Bartlett
         정규 N(0,1)      3.0    0.0370    0.0518
         균등 U(0,1)      1.8    0.0286    0.0032
         지수 Exp(1)      9.0    0.0458    0.3922
              t(5)      9.0    0.0422    0.2280
      로그정규 LN(0,1)    113.9    0.0382    0.6866
                폭 (최대 - 최소)    0.0172    0.6834

import numpy as np
from scipy.stats import levene

rng = np.random.default_rng(0)
n_sims, n, alpha = 5000, 20, 0.05

for center in ['mean', 'median']:
    rej = 0
    for _ in range(n_sims):
        g1 = rng.lognormal(0, 1, n)
        g2 = rng.lognormal(0, 1, n)
        g3 = rng.lognormal(0, 1, n)
        _, p = levene(g1, g2, g3, center=center)
        if p < alpha:
            rej += 1
    print(f"Levene ({center:6s}): Type I error = {rej/n_sims:.4f}")

Levene (mean  ): Type I error = 0.2470
Levene (median): Type I error = 0.0374

import numpy as np
from scipy.stats import levene
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n_sims, alpha = 5000, 0.05
sizes = [10, 20, 50, 100]
rates = []

for n in sizes:
    rej = 0
    for _ in range(n_sims):
        g1 = rng.lognormal(0, 1, n)
        g2 = rng.lognormal(0, 1, n)
        g3 = rng.lognormal(0, 1, n)
        _, p = levene(g1, g2, g3, center='median')
        if p < alpha:
            rej += 1
    rates.append(rej / n_sims)
    print(f"n={n:4d}: Type I error = {rej/n_sims:.4f}")

plt.figure(figsize=(6, 3))
plt.plot(sizes, rates, "o-")
plt.axhline(0.05, ls="--", color="red")
plt.xlabel("Group size n")
plt.ylabel("Type I error rate")
plt.title("Brown-Forsythe under lognormal data")
plt.tight_layout()
plt.show()

n=  10: Type I error = 0.0378
n=  20: Type I error = 0.0382
n=  50: Type I error = 0.0420
n= 100: Type I error = 0.0462

import numpy as np
from scipy.stats import levene, t as tdist

rng = np.random.default_rng(0)
n_sims, n, alpha = 5000, 20, 0.05

for center in ['mean', 'median']:
    rej = 0
    for _ in range(n_sims):
        g1 = tdist(df=3).rvs(n, random_state=rng)
        g2 = tdist(df=3).rvs(n, random_state=rng)
        g3 = tdist(df=3).rvs(n, random_state=rng)
        _, p = levene(g1, g2, g3, center=center)
        if p < alpha:
            rej += 1
    print(f"t(3), Levene ({center:6s}): {rej/n_sims:.4f}")

t(3), Levene (mean  ): 0.0640
t(3), Levene (median): 0.0340

import numpy as np
from scipy.stats import levene

rng = np.random.default_rng(0)
n_sims, n, alpha = 5000, 30, 0.05

for dist in ["normal", "lognormal"]:
    rej = 0
    for _ in range(n_sims):
        if dist == "normal":
            g1 = rng.normal(0, 1.0, n)
            g2 = rng.normal(0, 1.0, n)
            g3 = rng.normal(0, 1.5, n)
        else:
            g1 = rng.lognormal(0, 1.0, n)
            g2 = rng.lognormal(0, 1.0, n)
            g3 = rng.lognormal(0, 1.5, n)
        _, p = levene(g1, g2, g3, center='median')
        if p < alpha:
            rej += 1
    print(f"{dist:10s}: power = {rej/n_sims:.4f}")

normal    : power = 0.5122
lognormal : power = 0.2086