Brown-Forsythe 검정¶

import numpy as np
from scipy import stats

# 2번 집단만 유난히 퍼져 있다. 5 와 30 처럼 멀리 떨어진 값이 섞여 있어
# 평균 중심 방법이 흔들리기 좋은 자료다.
group1 = [8, 10, 12, 9, 11]
group2 = [5, 30, 18, 22, 15]
group3 = [14, 16, 15, 17, 14]

print("variances:", [round(np.var(g, ddof=1), 2)
                     for g in (group1, group2, group3)])

# (1) 어느 집단에서 평균과 중앙값이 갈리는가.
print(f"\n{'집단':>5}{'평균':>9}{'중앙값':>9}{'차':>8}")
for i, g in enumerate((group1, group2, group3), start=1):
    a = np.array(g, dtype=float)
    print(f"{i:>5}{a.mean():>9.2f}{np.median(a):>9.2f}{a.mean() - np.median(a):>8.2f}")

# Brown-Forsythe 는 각 값에서 제 집단의 중앙값을 뺀 절대편차에 분산분석을
# 돌린다. 중심을 평균이 아니라 중앙값으로 잡는 것이 전부인데, 그 한 가지
# 차이로 치우친 자료와 이상치에 훨씬 잘 버틴다.
stat, p_value = stats.levene(group1, group2, group3, center='median')
print(f"\nBrown-Forsythe W*:      {stat:.4f}, p = {p_value:.4f}")

# 원래 Levene 은 평균을 중심으로 쓴다. 이상치에 끌려간다.
s_m, p_m = stats.levene(group1, group2, group3, center='mean')
print(f"Levene (mean-centered): {s_m:.4f}, p = {p_m:.4f}")

# Bartlett 은 정규성을 전제하므로 이런 자료에서 가장 많이 어긋난다.
s_b, p_b = stats.bartlett(group1, group2, group3)
print(f"Bartlett:               {s_b:.4f}, p = {p_b:.6f}")

# 집단 3 의 다섯째 값만 14 -> 13 으로 바꾸면 세 집단 모두 평균 = 중앙값이 된다.
group3b = [14, 16, 15, 17, 13]
print(f"\n집단 3 을 {group3b} 로 바꾸면 "
      f"평균 {np.mean(group3b):.1f} = 중앙값 {np.median(group3b):.1f}")
a = stats.levene(group1, group2, group3b, center='median')[0]
b = stats.levene(group1, group2, group3b, center='mean')[0]
print(f"  W*(중앙값) = {a:.12f}")
print(f"  W (평균)   = {b:.12f}   차 = {abs(a - b):.2e}")

# (2) Bartlett 통계량을 정의대로 조립한다.
groups = [np.array(g, dtype=float) for g in (group1, group2, group3)]
nu = np.array([len(g) - 1 for g in groups])
s2 = np.array([g.var(ddof=1) for g in groups])
N, k = sum(len(g) for g in groups), len(groups)
sp2 = (nu * s2).sum() / nu.sum()
num = nu.sum() * np.log(sp2) - (nu * np.log(s2)).sum()
C = 1 + ((1 / nu).sum() - 1 / nu.sum()) / (3 * (k - 1))
print(f"\n합동분산 s_p^2 = {sp2:.6f}")
print(f"분자 = {nu.sum()} * ln(s_p^2) - sum nu_i ln(s_i^2) = {num:.6f}")
print(f"보정계수 C = {C:.6f}")
print(f"T = {num:.6f} / {C:.6f} = {num / C:.6f}   (scipy {s_b:.6f}, 차 {abs(num / C - s_b):.1e})")
print(f"p = chi2(2).sf(T) = {stats.chi2(k - 1).sf(num / C):.6e}")

# 두 검정이 보는 '신호'의 척도가 다르다.
Z = [np.abs(g - np.median(g)) for g in groups]
zbar = np.array([z.mean() for z in Z])
print(f"\n표본분산        {np.round(s2, 2)}   최대/최소 = {s2.max() / s2.min():.2f}")
print(f"Zbar (중앙값)   {np.round(zbar, 2)}   최대/최소 = {zbar.max() / zbar.min():.2f}")
print(f"  로그 폭:  ln(분산비) = {np.log(s2.max() / s2.min()):.4f},  "
      f"ln(Zbar 비) = {np.log(zbar.max() / zbar.min()):.4f}")
print(f"  sqrt(분산비) = {np.sqrt(s2.max() / s2.min()):.4f}  (Zbar 비 {zbar.max() / zbar.min():.4f} 과 견주라)")
ss_w = sum(((z - z.mean()) ** 2).sum() for z in Z)
print(f"\nBF 의 분모 SS_within(Z) = {ss_w:.3f}  "
      f"(집단별 {[round(((z - z.mean()) ** 2).sum(), 2) for z in Z]})")
n_i = np.array([len(z) for z in Z])
grand = np.concatenate(Z).mean()
ss_b = (n_i * (zbar - grand) ** 2).sum()
print(f"BF 의 분자 SS_between(Z) = {ss_b:.4f}  ->  W* = 12*{ss_b:.4f}/(2*{ss_w:.1f}) = {12 * ss_b / (2 * ss_w):.4f}")
print(f"집단 2 가 분모에서 차지하는 몫 = {((Z[1] - Z[1].mean()) ** 2).sum() / ss_w:.4f}")
print(f"F(0.95; 2, 12) = {stats.f(2, 12).ppf(0.95):.4f}")

variances: [2.5, 84.5, 1.7]

   집단       평균      중앙값       차
    1    10.00    10.00    0.00
    2    18.00    18.00    0.00
    3    15.20    15.00    0.20

Brown-Forsythe W*:      4.0754, p = 0.0446
Levene (mean-centered): 4.0610, p = 0.0450
Bartlett:               15.3946, p = 0.000454

집단 3 을 [14, 16, 15, 17, 13] 로 바꾸면 평균 15.0 = 중앙값 15.0
  W*(중앙값) = 3.896253602305
  W (평균)   = 3.896253602305   차 = 0.00e+00

합동분산 s_p^2 = 29.566667
분자 = 12 * ln(s_p^2) - sum nu_i ln(s_i^2) = 17.105089
보정계수 C = 1.111111
T = 17.105089 / 1.111111 = 15.394580   (scipy 15.394580, 차 0.0e+00)
p = chi2(2).sf(T) = 4.540560e-04

표본분산        [ 2.5 84.5  1.7]   최대/최소 = 49.71
Zbar (중앙값)   [1.2 6.4 1. ]   최대/최소 = 6.40
  로그 폭:  ln(분산비) = 3.9061,  ln(Zbar 비) = 1.8563
  sqrt(분산비) = 7.0502  (Zbar 비 6.4000 과 견주라)

BF 의 분모 SS_within(Z) = 138.000  (집단별 [2.8, 133.2, 2.0])
BF 의 분자 SS_between(Z) = 93.7333  ->  W* = 12*93.7333/(2*138.0) = 4.0754
집단 2 가 분모에서 차지하는 몫 = 0.9652
F(0.95; 2, 12) = 3.8853

import numpy as np
from scipy import stats

rng = np.random.default_rng(4)
R, n = 5000, 25

cases = [
    ("Normal",         lambda: rng.normal(0, 1, n)),
    ("Lognormal(0,1)", lambda: rng.lognormal(0, 1, n)),
    ("t(3)",           lambda: rng.standard_t(3, n)),
]

for name, gen in cases:
    c_mean = c_med = 0
    for _ in range(R):
        g = [gen() for _ in range(3)]
        c_mean += stats.levene(*g, center='mean')[1] < 0.05
        c_med += stats.levene(*g, center='median')[1] < 0.05
    print(f"{name:>16}: mean-centered {c_mean/R:.4f}, "
          f"median-centered {c_med/R:.4f}")

          Normal: mean-centered 0.0592, median-centered 0.0418
  Lognormal(0,1): mean-centered 0.2658, median-centered 0.0376
            t(3): mean-centered 0.0622, median-centered 0.0402

import numpy as np
from scipy import stats

# 2집단에는 극단값이 하나 있어 평균이 나머지 무리에서 멀리 끌려간다
g1 = [10, 11, 12, 13, 14, 11, 12]
g2 = [10, 11, 12, 13, 14, 11, 60]   # 60 is the outlier
g3 = [10, 11, 12, 13, 14, 11, 12]

for g in (g1, g2, g3):
    print(f"mean = {np.mean(g):7.3f}, median = {np.median(g):5.1f}, "
          f"var = {np.var(g, ddof=1):8.2f}")

print(f"\nLevene (mean):   p = {stats.levene(g1, g2, g3, center='mean')[1]:.4f}")
print(f"Brown-Forsythe:  p = {stats.levene(g1, g2, g3, center='median')[1]:.4f}")
print(f"Bartlett:        p = {stats.bartlett(g1, g2, g3)[1]:.4g}")

mean =  11.857, median =  12.0, var =     1.81
mean =  18.714, median =  12.0, var =   333.24
mean =  11.857, median =  12.0, var =     1.81

Levene (mean):   p = 0.0224
Brown-Forsythe:  p = 0.3723
Bartlett:        p = 2.018e-09