Brown-Forsythe 검정 (scipy)¶

import numpy as np
from scipy.stats import levene

# 세 집단의 평균은 크게 다르지만 퍼짐은 거의 같다.
g1 = np.array([12, 15, 14, 10, 13, 14, 12, 11], dtype=float)
g2 = np.array([22, 25, 20, 18, 24, 23, 19, 21], dtype=float)
g3 = np.array([32, 35, 34, 30, 33, 34, 32, 31], dtype=float)

print("variances:", [round(np.var(g, ddof=1), 4) for g in (g1, g2, g3)])

# scipy 에는 brown_forsythe 라는 함수가 따로 없다. levene 에
# center='median' 을 주는 것이 곧 Brown-Forsythe 이고, 이것이 기본값이다.
W, p = levene(g1, g2, g3, center='median')
print(f"Brown-Forsythe W = {W:.6f}, p-value = {p:.6f}")

W_mean, p_mean = levene(g1, g2, g3, center='mean')
print(f"Levene (mean)  W = {W_mean:.6f}, p-value = {p_mean:.6f}")

# (1) g3 가 g1 의 평행이동인가.
print(f"\ng3 - g1 = {g3 - g1}   ->  g3 = g1 + {int((g3 - g1)[0])} 인가: "
      f"{np.all(g3 - g1 == 20)}")

# (1)(2) 닫힌 꼴 W = 56 (zbar1 - zbar2)^2 / (2 S1 + S2) 를 두 중심에서 조립한다.
print(f"\n{'중심':>7}{'zbar1':>9}{'zbar2':>9}{'S1':>9}{'S2':>9}"
      f"{'분자 56(dz)^2':>15}{'분모 2S1+S2':>13}{'W':>11}")
for name, c in [("median", np.median), ("mean", np.mean)]:
    z1, z2, z3 = (np.abs(g - c(g)) for g in (g1, g2, g3))
    assert np.array_equal(z1, z3)          # g3 = g1 + 20 이므로 편차가 같다
    a, b = z1.mean(), z2.mean()
    S1 = ((z1 - a) ** 2).sum()
    S2 = ((z2 - b) ** 2).sum()
    num, den = 56 * (a - b) ** 2, 2 * S1 + S2
    print(f"{name:>7}{a:>9.4f}{b:>9.4f}{S1:>9.4f}{S2:>9.4f}"
          f"{num:>15.5f}{den:>13.4f}{num / den:>11.6f}")

# 왜 zbar1 이 두 중심에서 같은가. 중앙값이 sum|y-c| 를 최소화하는 구간을 본다.
s = np.sort(g1)
print(f"\ng1 정렬 = {s}")
print(f"  중앙값 = {np.median(g1):.3f},  평균 = {g1.mean():.3f}")
print(f"  최소화 구간 [x_(4), x_(5)] = [{s[3]:.0f}, {s[4]:.0f}] 에 평균이 드는가: "
      f"{s[3] <= g1.mean() <= s[4]}")
for c in (s[3], np.median(g1), g1.mean(), s[4]):
    print(f"    c = {c:7.3f}  ->  sum|g1 - c| = {np.abs(g1 - c).sum():.6f}")
print(f"\ng2 의 평균 {g2.mean():.3f} = 중앙값 {np.median(g2):.3f} 인가: "
      f"{g2.mean() == np.median(g2)}")

# (2) 분자가 같으니 두 W 의 비는 분모의 비와 같아야 한다.
print(f"\nW(mean)/W(median) = {W_mean / W:.6f}")
print(f"분모의 비 19.75/19.50 = {19.75 / 19.50:.6f}")

variances: [2.8393, 6.0, 2.8393]
Brown-Forsythe W = 1.107595, p-value = 0.348894
Levene (mean)  W = 1.121795, p-value = 0.344444

g3 - g1 = [20. 20. 20. 20. 20. 20. 20. 20.]   ->  g3 = g1 + 20 인가: True

     중심    zbar1    zbar2       S1       S2    분자 56(dz)^2    분모 2S1+S2          W
 median   1.3750   2.0000   4.8750  10.0000       21.87500      19.7500   1.107595
   mean   1.3750   2.0000   4.7500  10.0000       21.87500      19.5000   1.121795

g1 정렬 = [10. 11. 12. 12. 13. 14. 14. 15.]
  중앙값 = 12.500,  평균 = 12.625
  최소화 구간 [x_(4), x_(5)] = [12, 13] 에 평균이 드는가: True
    c =  12.000  ->  sum|g1 - c| = 11.000000
    c =  12.500  ->  sum|g1 - c| = 11.000000
    c =  12.625  ->  sum|g1 - c| = 11.000000
    c =  13.000  ->  sum|g1 - c| = 11.000000

g2 의 평균 21.500 = 중앙값 21.500 인가: True

W(mean)/W(median) = 1.012821
분모의 비 19.75/19.50 = 1.012821

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

x1 = np.array([2, 4, 6.])
x2 = np.array([1, 5, 9, 13.])
print(levene(x1, x2, center='median'))
print(f"F critical (1, 5) at alpha=0.05: {f.ppf(0.95, 1, 5):.4f}")

LeveneResult(statistic=3.2653061224489797, pvalue=0.13057297707695592)
F critical (1, 5) at alpha=0.05: 6.6079

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

rng = np.random.default_rng(42)
rej_bart, rej_bf = 0, 0
n_sims = 3000

for _ in range(n_sims):
    g1 = rng.lognormal(0, 1, 30)
    g2 = rng.lognormal(0, 1, 30)
    g3 = rng.lognormal(0, 1, 30)
    _, p_b = bartlett(g1, g2, g3)
    _, p_bf = levene(g1, g2, g3, center='median')
    if p_b < 0.05:
        rej_bart += 1
    if p_bf < 0.05:
        rej_bf += 1

print(f"Bartlett false-positive rate:       {rej_bart/n_sims:.4f}")
print(f"Brown-Forsythe false-positive rate: {rej_bf/n_sims:.4f}")

Bartlett false-positive rate:       0.7303
Brown-Forsythe false-positive rate: 0.0307

import numpy as np
from scipy.stats import levene

rng = np.random.default_rng(42)

for n in [25, 100]:
    g1 = rng.normal(0, 1.0, n)
    g2 = rng.normal(0, 1.5, n)
    g3 = rng.normal(0, 2.0, n)
    g4 = rng.normal(0, 3.0, n)
    W, p = levene(g1, g2, g3, g4, center='median')
    print(f"n={n:3d}: W={W:.4f}, p={p:.6g}")

n= 25: W=8.3766, p=5.29567e-05
n=100: W=29.2350, p=4.28901e-17

import numpy as np
from scipy.stats import levene, trim_mean, f_oneway

rng = np.random.default_rng(0)
rej_med, rej_trim = 0, 0
n_sims = 3000

for _ in range(n_sims):
    g1 = rng.standard_t(3, 30)
    g2 = rng.standard_t(3, 30)
    g3 = rng.standard_t(3, 30)

    # 중앙값 중심 — Brown-Forsythe
    _, p_med = levene(g1, g2, g3, center='median')
    if p_med < 0.05:
        rej_med += 1

    # 절사평균 중심. 편차는 잘라 내지 않은 전체 표본에서 잰다
    # from the trimmed mean (data are NOT discarded)
    z_groups = [np.abs(g - trim_mean(g, 0.1)) for g in (g1, g2, g3)]
    _, p_trim = f_oneway(*z_groups)
    if p_trim < 0.05:
        rej_trim += 1

print(f"Median centering FPR:       {rej_med/n_sims:.4f}")
print(f"Trimmed mean centering FPR: {rej_trim/n_sims:.4f}")

Median centering FPR:       0.0383
Trimmed mean centering FPR: 0.0447