대수정규에서의 비교: 치우침이 극단일 때¶

import numpy as np
from scipy import stats

# 1부: 제1종 오류율. 세 집단을 모두 같은 로그정규에서 뽑으므로 분산은
# 참으로 같다. 기각 비율이 0.05 를 크게 넘는 검정은 정규성에 매여 있다는 뜻이다.
rng = np.random.default_rng(42)
n_sims, n, alpha = 2000, 30, 0.05
results = {"Bartlett": 0, "Levene (mean)": 0,
           "Brown-Forsythe": 0, "Fligner-Killeen": 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 = stats.bartlett(g1, g2, g3)
    if p < alpha:
        results["Bartlett"] += 1

    _, p = stats.levene(g1, g2, g3, center='mean')
    if p < alpha:
        results["Levene (mean)"] += 1

    _, p = stats.levene(g1, g2, g3, center='median')
    if p < alpha:
        results["Brown-Forsythe"] += 1

    _, p = stats.fligner(g1, g2, g3)
    if p < alpha:
        results["Fligner-Killeen"] += 1

for name, count in results.items():
    print(f"{name:20s}: false-positive rate = {count/n_sims:.4f}")

# ---- (1) 몬테카를로 오차로 재고, 바틀렛의 극한값을 계산한다 ----
# 반복 2000 회이므로 0.05 근처의 표준오차는 sqrt(.05*.95/2000) = 0.0049 다.
print(f"\n명목 0.05, 반복 {n_sims} 회의 몬테카를로 표준오차 = "
      f"{np.sqrt(0.05 * 0.95 / n_sims):.5f}")
print(f"{'검정':>18}{'크기':>9}{'SE':>9}{'0.05 로부터':>14}{'명목의 배수':>12}")
for name, count in results.items():
    r = count / n_sims
    se = np.sqrt(r * (1 - r) / n_sims)
    print(f"{name:>18}{r:>9.4f}{se:>9.5f}{(r - 0.05) / se:>11.1f} SE{r / 0.05:>12.2f}")

# 로그정규(0,1) 의 첨도와, 그것이 정하는 바틀렛의 극한 크기
b2 = np.exp(4) + 2 * np.exp(3) + 3 * np.exp(2) - 3
crit = stats.chi2(2).ppf(1 - alpha)
print(f"\n로그정규(0,1) 의 첨도 beta_2 = {b2:.4f}")
print(f"예측: T / ((beta_2-1)/2) -> chi2_2 이므로 E[T] -> beta_2 - 1 = {b2 - 1:.2f}")
print(f"극한 크기 = P(chi2_2 > 2 * {crit:.4f} / (beta_2 - 1)) = "
      f"{stats.chi2(2).sf(2 * crit / (b2 - 1)):.4f}")
print(f"검산 (정규, beta_2 = 3): P(chi2_2 > 2*{crit:.4f}/2) = "
      f"{stats.chi2(2).sf(2 * crit / 2):.4f}")

# ---- (2) 표본을 키우면 나아지는가. 세 검정을 n 으로 훑는다 ----
# rng 와 따로 쓰는 생성기여야 위쪽 출력이 바뀌지 않는다.
rng_b = np.random.default_rng(2024)
B = 2000
print(f"\n로그정규 세 집단, 등분산, 반복 {B} 회")
print(f"{'n':>6}{'Bartlett':>20}{'Brown-Forsythe':>22}{'Fligner-Killeen':>22}")
for nn in (30, 120, 480):
    rej = {"ba": 0, "bf": 0, "fk": 0}
    for _ in range(B):
        g = [rng_b.lognormal(0, 1, nn) for _ in range(3)]
        rej["ba"] += stats.bartlett(*g)[1] < alpha
        rej["bf"] += stats.levene(*g, center='median')[1] < alpha
        rej["fk"] += stats.fligner(*g)[1] < alpha
    cells = ""
    for k in ("ba", "bf", "fk"):
        r = rej[k] / B
        cells += f"{r:>14.4f} (SE {np.sqrt(r * (1 - r) / B):.4f})"
    print(f"{nn:>6}{cells}")

Bartlett            : false-positive rate = 0.7320
Levene (mean)       : false-positive rate = 0.2605
Brown-Forsythe      : false-positive rate = 0.0295
Fligner-Killeen     : false-positive rate = 0.1140

명목 0.05, 반복 2000 회의 몬테카를로 표준오차 = 0.00487
                검정       크기       SE      0.05 로부터      명목의 배수
          Bartlett   0.7320  0.00990       68.9 SE       14.64
     Levene (mean)   0.2605  0.00981       21.4 SE        5.21
    Brown-Forsythe   0.0295  0.00378       -5.4 SE        0.59
   Fligner-Killeen   0.1140  0.00711        9.0 SE        2.28

로그정규(0,1) 의 첨도 beta_2 = 113.9364
예측: T / ((beta_2-1)/2) -> chi2_2 이므로 E[T] -> beta_2 - 1 = 112.94
극한 크기 = P(chi2_2 > 2 * 5.9915 / (beta_2 - 1)) = 0.9483
검산 (정규, beta_2 = 3): P(chi2_2 > 2*5.9915/2) = 0.0500

로그정규 세 집단, 등분산, 반복 2000 회
     n            Bartlett        Brown-Forsythe       Fligner-Killeen
    30        0.7145 (SE 0.0101)        0.0400 (SE 0.0044)        0.1150 (SE 0.0071)
   120        0.8025 (SE 0.0089)        0.0430 (SE 0.0045)        0.1495 (SE 0.0080)
   480        0.8715 (SE 0.0075)        0.0450 (SE 0.0046)        0.1700 (SE 0.0084)

# 2부: 검정력. 이번에는 정규모집단이고 표준편차를 1, 1.5, 2 로 실제로
# 다르게 두었다. 오류율만 보고 검정을 고를 수 없는 까닭이 여기 있다 —
# 보수적인 검정은 오류율은 낮지만 차이도 잘 못 잡는다.
results_power = {"Bartlett": 0, "Levene (mean)": 0,
                 "Brown-Forsythe": 0, "Fligner-Killeen": 0}

for _ in range(n_sims):
    g1 = rng.normal(0, 1.0, n)
    g2 = rng.normal(0, 1.5, n)
    g3 = rng.normal(0, 2.0, n)

    _, p = stats.bartlett(g1, g2, g3)
    if p < alpha:
        results_power["Bartlett"] += 1

    _, p = stats.levene(g1, g2, g3, center='mean')
    if p < alpha:
        results_power["Levene (mean)"] += 1

    _, p = stats.levene(g1, g2, g3, center='median')
    if p < alpha:
        results_power["Brown-Forsythe"] += 1

    _, p = stats.fligner(g1, g2, g3)
    if p < alpha:
        results_power["Fligner-Killeen"] += 1

for name, count in results_power.items():
    print(f"{name:20s}: power = {count/n_sims:.4f}")

# ---- 크기를 맞춘 뒤의 검정력. 따로 쓰는 생성기여야 위쪽이 바뀌지 않는다 ----
# 정규 H0 에서 각 검정의 p-값 5% 분위수를 문턱으로 다시 잡으면, 네 검정의
# 실제 크기가 정확히 0.05 가 된다. 그 문턱으로 검정력을 다시 센다.
rng_c = np.random.default_rng(777)
names = ["Bartlett", "Levene (mean)", "Brown-Forsythe", "Fligner-Killeen"]


def four_pvalues(g):
    return [stats.bartlett(*g)[1],
            stats.levene(*g, center='mean')[1],
            stats.levene(*g, center='median')[1],
            stats.fligner(*g)[1]]


P0 = np.array([four_pvalues([rng_c.normal(0, 1, n) for _ in range(3)])
               for _ in range(n_sims)])                      # 정규 H0
P1 = np.array([four_pvalues([rng_c.normal(0, s, n) for s in (1.0, 1.5, 2.0)])
               for _ in range(n_sims)])                      # 정규 H1
thresh = np.quantile(P0, alpha, axis=0)

print(f"\n{'검정':>18}{'정규 크기':>11}{'보정 문턱':>11}{'명목 검정력':>13}{'보정 검정력':>13}")
for i, nm in enumerate(names):
    size = (P0[:, i] < alpha).mean()
    pw = (P1[:, i] < alpha).mean()
    pw_cal = (P1[:, i] < thresh[i]).mean()
    print(f"{nm:>18}{size:>11.4f}{thresh[i]:>11.4f}{pw:>13.4f}{pw_cal:>13.4f}")
print(f"\n검정력 0.85 쯤의 몬테카를로 표준오차 = "
      f"{np.sqrt(0.85 * 0.15 / n_sims):.5f}")

Bartlett            : power = 0.9215
Levene (mean)       : power = 0.8490
Brown-Forsythe      : power = 0.8155
Fligner-Killeen     : power = 0.7810

                검정      정규 크기      보정 문턱       명목 검정력       보정 검정력
          Bartlett     0.0485     0.0511       0.9220       0.9240
     Levene (mean)     0.0515     0.0471       0.8505       0.8450
    Brown-Forsythe     0.0420     0.0597       0.8185       0.8480
   Fligner-Killeen     0.0420     0.0613       0.7845       0.8160

검정력 0.85 쯤의 몬테카를로 표준오차 = 0.00798

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n_sims, n, alpha = 5000, 100, 0.05
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 = stats.bartlett(g1, g2, g3)
    if p < alpha:
        rej += 1

print(f"Bartlett false-positive rate (n=100): {rej/n_sims:.4f}")

Bartlett false-positive rate (n=100): 0.8134

import numpy as np
from scipy import stats

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

for _ in range(n_sims):
    g1 = rng.lognormal(0, 1, n)
    g2 = rng.lognormal(0, 1, n)
    F = np.var(g1, ddof=1) / np.var(g2, ddof=1)
    p_f = 2 * min(stats.f(n-1, n-1).cdf(F), stats.f(n-1, n-1).sf(F))
    _, p_b = stats.bartlett(g1, g2)
    if p_f < alpha:
        rej_f += 1
    if p_b < alpha:
        rej_b += 1

print(f"F-test:   {rej_f/n_sims:.4f}")
print(f"Bartlett: {rej_b/n_sims:.4f}")

F-test:   0.5108
Bartlett: 0.5108

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n_sims, n, alpha = 5000, 30, 0.05
power = {"Bartlett": 0, "Levene": 0, "BF": 0, "FK": 0}

for _ in range(n_sims):
    g1 = stats.t(df=5).rvs(n, random_state=rng) * 1.0
    g2 = stats.t(df=5).rvs(n, random_state=rng) * 1.5
    g3 = stats.t(df=5).rvs(n, random_state=rng) * 2.0
    _, p = stats.bartlett(g1, g2, g3)
    if p < alpha:
        power["Bartlett"] += 1
    _, p = stats.levene(g1, g2, g3, center='mean')
    if p < alpha:
        power["Levene"] += 1
    _, p = stats.levene(g1, g2, g3, center='median')
    if p < alpha:
        power["BF"] += 1
    _, p = stats.fligner(g1, g2, g3)
    if p < alpha:
        power["FK"] += 1

for name, count in power.items():
    print(f"{name:10s}: {count/n_sims:.4f}")

Bartlett  : 0.8578
Levene    : 0.7096
BF        : 0.6584
FK        : 0.6616