Bartlett 검정 비정규성 민감도¶

import numpy as np
from scipy.stats import bartlett, levene, fligner

rng = np.random.default_rng(0)


def simulate_once(n=20, sigmas=(1.0, 1.0, 1.0), skew=True):
    """치우친 로그정규 또는 정규에서 k 개 집단을 만든다.

    셋 다 같은 sigma 를 주므로 분산은 참으로 같다. 그런데도 자료가
    치우쳐 있으면 검정이 기각해 버리는지를 보려는 것이다.
    """
    if skew:
        groups = [rng.lognormal(mean=0.0, sigma=s, size=n) for s in sigmas]
    else:
        groups = [rng.normal(loc=0.0, scale=s, size=n) for s in sigmas]
    return groups


def trial(n=20, sigmas=(1.0, 1.0, 1.0), skew=True):
    """한 벌의 모의자료에 세 검정을 모두 돌린다."""
    g1, g2, g3 = simulate_once(n=n, sigmas=sigmas, skew=skew)
    _, p_bartlett = bartlett(g1, g2, g3)
    _, p_levene = levene(g1, g2, g3, center='mean')
    _, p_fligner = fligner(g1, g2, g3)
    return p_bartlett, p_levene, p_fligner


# 분산이 참으로 같으므로, 아래 비율은 모두 0.05 근처여야 옳다.
# Bartlett 만 크게 벗어난다면 그것이 이 검정의 약점이다.
n_sims, alpha = 5000, 0.05
ps_b, ps_lv, ps_fl = [], [], []

for _ in range(n_sims):
    p_b, p_l, p_f = trial(n=20, sigmas=(1.0, 1.0, 1.0), skew=True)
    ps_b.append(p_b)
    ps_lv.append(p_l)
    ps_fl.append(p_f)

fp_b = np.mean(np.array(ps_b) < alpha)
fp_lv = np.mean(np.array(ps_lv) < alpha)
fp_fl = np.mean(np.array(ps_fl) < alpha)

print("False-positive rates under skewed (lognormal) data:")
print(f"  Bartlett       : {fp_b:.4f}")
print(f"  Levene (mean)  : {fp_lv:.4f}")
print(f"  Fligner-Killeen: {fp_fl:.4f}")

False-positive rates under skewed (lognormal) data:
  Bartlett       : 0.6748
  Levene (mean)  : 0.2470
  Fligner-Killeen: 0.1028

import numpy as np
from scipy import stats

k, alpha = 3, 0.05
crit = stats.chi2.ppf(1 - alpha, k - 1)

# 로그정규(0,1) 의 초과첨도와 그에서 나오는 팽창 배율
g2 = np.exp(4) + 2 * np.exp(3) + 3 * np.exp(2) - 6
beta2_m1 = g2 + 2                      # = gamma2 + 2 = beta2 - 1
print(f"Lognormal(0,1):  gamma2 = {g2:.3f},  beta2 - 1 = {beta2_m1:.3f}")
print(f"  Var(ln S^2) 의 팽창 배율 = (beta2-1)/2 = {beta2_m1 / 2:.2f}"
      f"   (표준편차로는 {np.sqrt(beta2_m1 / 2):.2f} 배)")
print(f"  극한 거짓 양성률 = P(chi2_2 > 2*{crit:.4f}/{beta2_m1:.3f})"
      f" = {stats.chi2.sf(2 * crit / beta2_m1, k - 1):.4f}")

# 표본을 키우면 그 극한으로 가는가. 반복 2000 회.
rng = np.random.default_rng(4)
M = 2000
print(f"\n{'n':>8}{'Bartlett 거짓 양성률':>24}")
for n, ch in [(20, 250), (100, 250), (1000, 100), (10000, 25)]:
    C = 1 + (k + 1) / (3 * k * (n - 1))
    hits = 0
    for _ in range(M // ch):
        S = rng.lognormal(0, 1, (ch, k, n)).var(axis=2, ddof=1)
        T = k * (n - 1) * np.log(S.mean(axis=1)
                                 / np.exp(np.log(S).mean(axis=1))) / C
        hits += int(np.sum(T > crit))
    print(f"{n:>8}{hits / M:>24.4f}")
print(f"\n반복 {M} 회이므로 몬테카를로 표준오차는 0.67 근처 "
      f"{np.sqrt(0.67 * 0.33 / M):.4f}, 0.93 근처 {np.sqrt(0.93 * 0.07 / M):.4f}")

Lognormal(0,1):  gamma2 = 110.936,  beta2 - 1 = 112.936
  Var(ln S^2) 의 팽창 배율 = (beta2-1)/2 = 56.47   (표준편차로는 7.51 배)
  극한 거짓 양성률 = P(chi2_2 > 2*5.9915/112.936) = 0.9483

       n         Bartlett 거짓 양성률
      20                  0.6655
     100                  0.8025
    1000                  0.8890
   10000                  0.9270

반복 2000 회이므로 몬테카를로 표준오차는 0.67 근처 0.0105, 0.93 근처 0.0057

import numpy as np
from scipy.stats import bartlett, levene, fligner

rng = np.random.default_rng(0)
n_sims, alpha = 5000, 0.05
rej = {"Bartlett": 0, "Levene": 0, "Fligner": 0}

for _ in range(n_sims):
    g1 = rng.normal(0, 1, 20)
    g2 = rng.normal(0, 1, 20)
    g3 = rng.normal(0, 1, 20)
    _, p = bartlett(g1, g2, g3)
    if p < alpha: rej["Bartlett"] += 1
    _, p = levene(g1, g2, g3, center='mean')
    if p < alpha: rej["Levene"] += 1
    _, p = fligner(g1, g2, g3)
    if p < alpha: rej["Fligner"] += 1

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

Bartlett: 0.0508
Levene: 0.0556
Fligner: 0.0364

import numpy as np
from scipy.stats import bartlett

rng = np.random.default_rng(0)
n_sims = 5000
rej = 0
for _ in range(n_sims):
    g1 = rng.lognormal(0, 1, 100)
    g2 = rng.lognormal(0, 1, 100)
    g3 = rng.lognormal(0, 1, 100)
    _, p = bartlett(g1, g2, g3)
    if p < 0.05:
        rej += 1
print(f"Bartlett FP rate (n=100): {rej/n_sims:.4f}")

Bartlett FP rate (n=100): 0.8134

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

rng = np.random.default_rng(0)
n_sims = 5000
rej = 0
for _ in range(n_sims):
    g1 = tdist(df=3).rvs(20, random_state=rng)
    g2 = tdist(df=3).rvs(20, random_state=rng)
    g3 = tdist(df=3).rvs(20, random_state=rng)
    _, p = bartlett(g1, g2, g3)
    if p < 0.05:
        rej += 1
print(f"Bartlett FP rate (t(3)): {rej/n_sims:.4f}")

Bartlett FP rate (t(3)): 0.4076

import numpy as np
from scipy.stats import bartlett, kurtosis, chi2

def bartlett_corrected(*groups):
    stat, _ = bartlett(*groups)
    k = len(groups)
    all_data = np.concatenate(groups)
    kappa_hat = kurtosis(all_data, fisher=True)
    ns = [len(g) for g in groups]
    N = sum(ns)
    correction = max(1 + kappa_hat / (2 * (k + 1)) *
                     (sum(1 / (n - 1) for n in ns) - 1 / (N - k)), 0.5)
    adjusted_stat = stat / correction
    return adjusted_stat, chi2(df=k - 1).sf(adjusted_stat)

Corrected Bartlett size, Normal data:      0.0510
Corrected Bartlett size, Lognormal data:   0.6488