분산에 대한 가능도비 검정¶

1단계. 합동분산:

2단계. LRT 통계량:

3단계. 보정인자:

4단계. Bartlett 통계량: \(T = 0.757 / 1.0327 = 0.733\).

5단계. \(\chi^2_{0.95, 1} = 3.841\)과 비교한다. \(0.733 < 3.841\)이므로 \(H_0\)을 기각하지 못한다(\(p = 0.392\)).

(1) \(C\) 는 자료를 보지 않는다.

에 들어 있는 것은 \(k\) 와 \(\nu_1, \ldots, \nu_k\) 뿐이고 \(S_i^2\) 은 어디에도 없다. 그러므로 설계가 정해지면 \(C\) 도 정해지고, 자료를 보기 전에 계산할 수 있다. 보정은 통계량을 자료와 무관한 상수로 나누는 일, 곧 눈금을 한 번 다시 매기는 일이다. 따라서

가 모든 자료에서 같은 값이고, 두 검정의 \(p\) 값 순서도 결코 뒤바뀌지 않는다.

늘 \(C > 1\) 이다. \(\nu_i \ge 1\) 이고 \(\sum_i \nu_i = N - k\) 이므로, \(k \ge 2\) 일 때

이다(가운데 부등호는 \(k \ge 2\) 라 \(\max_i \nu_i < \sum_i \nu_i\) 이기 때문이다). 괄호 안이 양수이므로 \(C > 1\) 이고, 따라서

가 언제나 성립한다. 보정은 늘 통계량을 줄이고, 그러므로 늘 기각을 덜 하게 만든다. 보정하지 않은 LRT 가 작은 표본에서 지나치게 자주 기각한다는 말과 방향이 맞는다.

등자유도일 때. \(\nu_i = \nu\) 이면 \(\sum 1/\nu_i = k/\nu\) 이고 \(N - k = k\nu\) 이므로

이다. \(k = 2\) 면 \(C = 1 + \dfrac{1}{2\nu}\) 로 아주 간단하다.

이것이 최솟값이다. 산술–조화평균 부등식(또는 코시–슈바르츠)에서

이고 등호는 \(\nu_i\) 가 모두 같을 때만 성립한다. 이를 넣으면

인데, \(N - k = k\nu\) 를 쓰면 이 하한이 바로 위의 등자유도 공식이다. 총자유도가 같다면 집단 크기가 고를 때 보정이 가장 작다. 거꾸로 말하면 집단 크기가 불균형할수록 카이제곱 근사가 더 나쁘고 보정이 더 커야 한다.

(2) 수로 확인한다. \(H_0\) 이 참이고 자료가 정규라면 \(\nu_i S_i^2/\sigma^2 \sim \chi^2_{\nu_i}\) 이 서로 독립이다. 통계량이 \(S_i^2\) 의 비에만 의존하므로 \(\sigma^2 = 1\) 로 두어도 일반성을 잃지 않고, 원자료를 만들 필요 없이 카이제곱에서 바로 뽑으면 된다.

import numpy as np
from scipy import stats

# 원자료 없이 집단별 크기와 표본분산만으로 계산할 수 있다.
n = np.array([15, 18])
s2 = np.array([22.4, 35.1])
k = len(n)
nu = n - 1
N = n.sum()

# 귀무가설(분산이 모두 같다) 아래에서의 공통 분산 추정값
s2_pooled = np.sum(nu * s2) / np.sum(nu)

# 가능도비 검정통계량. 제약 있는 최대가능도와 없는 것의 비에 로그를 씌운
# 값이며, 각 집단의 분산이 공통값에서 얼마나 벗어났는지를 재는 셈이다.
lrt = np.sum(nu * np.log(s2_pooled / s2))

# 그대로 두면 카이제곱 근사가 작은 표본에서 어긋난다. Bartlett 보정은
# 통계량을 이 인자로 나눠 근사를 개선한다.
C = 1 + (1 / (3 * (k - 1))) * (np.sum(1 / nu) - 1 / np.sum(nu))

# 보정한 통계량
T = lrt / C

# p-값
p_value = stats.chi2.sf(T, k - 1)

print(f"Pooled variance: {s2_pooled:.4f}")
print(f"LRT statistic (uncorrected): {lrt:.4f}")
print(f"Correction factor: {C:.4f}")
print(f"Bartlett statistic (corrected): {T:.4f}")
print(f"P-value: {p_value:.4f}")

# === (1) C 가 자료를 보지 않는다. 불균형이 C 를 키우는가 ===
def corr_factor(nu):
    nu = np.asarray(nu, dtype=float)
    kk = len(nu)
    return 1 + (1 / (3 * (kk - 1))) * (np.sum(1 / nu) - 1 / np.sum(nu))

print(f"\n총자유도 N-k = {nu.sum()} 를 2 집단에 어떻게 쪼개는가")
print(f"  균형 하한 1+(k+1)/(3(N-k)) = {1 + (k + 1) / (3 * nu.sum()):.6f}")
for split in ([15, 16], [14, 17], [10, 21], [5, 26], [1, 30]):
    print(f"  nu = {str(split):10s} -> C = {corr_factor(split):.6f}")

# === (2) 보정이 1종 오류율을 고치는가 ===
def sizes(nu, M=200_000, alpha=0.05, rng=None):
    """H0 아래에서 nu_i S_i^2 ~ chi2(nu_i) 이므로 바로 뽑는다."""
    nu = np.asarray(nu, dtype=float)
    kk = len(nu)
    W = rng.chisquare(nu, size=(M, kk))       # nu_i * S_i^2  (sigma^2 = 1)
    s2b = W / nu                              # S_i^2
    sp = W.sum(axis=1) / nu.sum()             # 합동분산
    stat = (nu * np.log(sp[:, None] / s2b)).sum(axis=1)
    crit = stats.chi2(kk - 1).ppf(1 - alpha)
    Cc = corr_factor(nu)
    return np.mean(stat > crit), np.mean(stat / Cc > crit), Cc

rng = np.random.default_rng(0)
M = 200_000
print(f"\n정규 자료, H0 참, alpha = 0.05, 반복 {M} 회 "
      f"(0.05 근처 표준오차 +-{np.sqrt(0.05 * 0.95 / M):.5f})")
print("  nu                    C      보정 전    보정 후   등자유도 공식")
for nus in ([4, 4], [3, 3, 3, 3], [4, 4, 4, 4, 4], [9, 9], [9, 9, 9],
            [14, 17], [49, 49]):
    a, b, Cc = sizes(nus, M=M, rng=rng)
    kk = len(nus)
    eq = (f"{1 + (kk + 1) / (3 * kk * nus[0]):.4f}"
          if len(set(nus)) == 1 else "-")
    print(f"  {str(nus):18s}{Cc:7.4f}   {a:.4f}     {b:.4f}   {eq:>8s}")

출력:

Pooled variance: 29.3645
LRT statistic (uncorrected): 0.7571
Correction factor: 1.0327
Bartlett statistic (corrected): 0.7332
P-value: 0.3919

총자유도 N-k = 31 를 2 집단에 어떻게 쪼개는가
  균형 하한 1+(k+1)/(3(N-k)) = 1.032258
  nu = [15, 16]   -> C = 1.032303
  nu = [14, 17]   -> C = 1.032665
  nu = [10, 21]   -> C = 1.038454
  nu = [5, 26]    -> C = 1.068734
  nu = [1, 30]    -> C = 1.333692

정규 자료, H0 참, alpha = 0.05, 반복 200000 회 (0.05 근처 표준오차 +-0.00049)
  nu                    C      보정 전    보정 후   등자유도 공식
  [4, 4]             1.1250   0.0650     0.0500     1.1250
  [3, 3, 3, 3]       1.1389   0.0736     0.0475     1.1389
  [4, 4, 4, 4, 4]    1.1000   0.0697     0.0493     1.1000
  [9, 9]             1.0556   0.0561     0.0496     1.0556
  [9, 9, 9]          1.0494   0.0572     0.0494     1.0494
  [14, 17]           1.0327   0.0531     0.0494          -
  [49, 49]           1.0102   0.0507     0.0494     1.0102

보기 1 의 손계산과 코드가 맞는다. 합동분산 \(29.3645\), 보정 전 \(0.7571\), \(C = 1.0327\), 보정 후 \(0.7332\), \(p = 0.3919\) 로 다섯 값 모두 보기 1 과 같다.

(1)의 등자유도 공식이 정확히 맞는다. 마지막 열의 \(1 + (k+1)/(3k\nu)\) 가 실제 \(C\) 와 네 자리까지 모두 같다. \(\nu = 4\), \(k = 2\) 에서 \(1 + 1/(2\cdot4) = 1.1250\), \(\nu = 3\), \(k = 4\) 에서 \(1 + 5/36 = 1.1389\) 이다.

균형이 \(C\) 를 최소로 한다는 것도 확인된다. 총자유도 \(31\) 을 두 집단에 쪼갤 때 하한이 \(1.032258\) 인데, 가장 고른 \([15,16]\) 이 \(1.032303\) 으로 거의 거기에 붙고, 불균형을 키우면 \([10,21]\) \(1.038454\), \([5,26]\) \(1.068734\), \([1,30]\) \(1.333692\) 로 올라간다. 한 집단이 자유도 1 밖에 없으면 보정이 \(33\%\) 나 된다. 보기 1 의 \([14,17]\) 은 거의 균형이라 보정이 \(3.3\%\) 로 작다.

보정은 정말로 1종 오류율을 고친다. 보정하지 않은 \(-2\ln\Lambda\) 의 크기가 \(\nu = 4\) 두 집단에서 \(0.0650\), 네 집단 \(\nu = 3\) 에서 \(0.0736\), 다섯 집단 \(\nu = 4\) 에서 \(0.0697\) 으로 모두 명목값을 넘는다. 같은 자료에 \(C\) 로 나누면 \(0.0500\), \(0.0475\), \(0.0493\) 이 된다. 20 만 번 반복에서 표준오차가 \(\pm0.00049\) 이므로 보정 전의 초과(\(0.0650\) 은 명목값보다 \(31\) 표준오차 위다)는 의심의 여지 없이 실재하고, 보정 후에는 거의 다 명목값으로 돌아온다.

자유도가 커지면 보정할 것도 줄어든다. \(\nu = 49\) 두 집단에서는 보정 전이 이미 \(0.0507\) 이고 보정 후 \(0.0494\) 이며 \(C\) 도 \(1.0102\) 로 1 에 가깝다. \(C\) 는 \(\nu \to \infty\) 에서 1 로 가고, Wilks 정리의 \(\chi^2_{k-1}\) 근사가 그 극한에서 정확해진다. 보정은 그 수렴을 앞당기는 유한표본 손질이다.

다만 완벽하지는 않다. 가장 작은 자유도인 \([3,3,3,3]\) 에서 보정 후가 \(0.0475\) 로 명목값보다 \(5\) 표준오차 낮다. 보정이 살짝 지나친 것이다. \(\nu = 3\) 은 집단당 관측값 4 개이니 그 정도는 당연하고, 실무에서 문제가 될 크기도 아니다(보정 전 \(0.0736\) 과 견주라). 그러나 "보정했으니 정확하다"가 아니라 "보정했으니 훨씬 낫다"가 맞는 말이다. 보정 후 값들이 \(0.0475\)–\(0.0500\) 으로 모두 명목값 쪽에서 조금 아래에 모이는 것을 보면, 이 보정은 평균적으로 아주 조금 보수적인 편이다.

그리고 이 모의실험은 정규성을 가정한 것이다. \(\nu_i S_i^2 \sim \chi^2_{\nu_i}\) 자체가 정규모집단에서만 참이기 때문이다. 정규성이 깨지면 \(C\) 로 나누는 일은 아무 도움이 되지 않는다. Bartlett 검정이 꼬리가 무거운 자료에서 무너지는 것은 보정의 실패가 아니라 가능도의 실패다. 15.4절 Bartlett 검정의 한계가 그 이야기다.

가능도비 통계량은

여기서 \(\hat{\sigma}^2 = \frac{1}{n}\sum(x_i - \bar{x})^2\)은 최대가능도추정량이다.

양변에 로그를 취하고 \(-2\)를 곱하면 \(r = \hat{\sigma}^2/\sigma_0^2\)에 대해

이며 \(H_0\) 아래에서 점근적으로 \(-2\ln\Lambda \sim \chi^2_1\)이다. 본문의 유도와 일치한다.

검산. \(\Lambda \leq 1\)이어야 하므로 \(\ln \Lambda \leq 0\), 곧 \(-2\ln\Lambda \geq 0\)이다. \(g(r) = r - 1 - \ln r\)가 \(r>0\)에서 항상 음이 아니므로 조건이 충족된다. \(\square\)

고전적 카이제곱 검정은 \(\chi^2 = (n-1)s^2/\sigma_0^2 \sim \chi^2_{n-1}\)을 쓰며 정규성 아래에서 정확하다. 가능도비 검정은 \(-2\ln\Lambda \sim \chi^2_1\)을 점근적으로 쓴다.

동등하지 않다. 유한표본에서 카이제곱 검정은 정확하고 LRT는 근사에 기댄다. 큰 \(n\)에서는 거의 같은 결과를 준다.

모의실험으로 확인해 보자(\(H_0\) 참, \(\alpha = 0.05\), 반복 20,000회).

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
R = 20000

for n in [10, 25, 50, 100, 500]:
    a = b = 0
    for _ in range(R):
        x = rng.normal(0, 1, n)
        r = x.var(ddof=0)                      # MLE / sigma0^2 with sigma0=1
        lrt = n * (r - 1 - np.log(r))
        a += lrt > stats.chi2.ppf(0.95, 1)

        c = (n - 1) * x.var(ddof=1)
        p = 2 * min(stats.chi2.cdf(c, n - 1), stats.chi2.sf(c, n - 1))
        b += p < 0.05
    print(f"n = {n:>4}: LRT size {a/R:.4f}, exact chi2 size {b/R:.4f}")

출력:

n =   10: LRT size 0.0731, exact chi2 size 0.0500
n =   25: LRT size 0.0586, exact chi2 size 0.0500
n =   50: LRT size 0.0558, exact chi2 size 0.0519
n =  100: LRT size 0.0515, exact chi2 size 0.0490
n =  500: LRT size 0.0537, exact chi2 size 0.0521

정확 카이제곱 검정은 모든 \(n\)에서 정확히 0.05이다(당연하다, 정확한 검정이므로). LRT는 \(n = 10\)에서 0.073으로 46% 부풀려져 있고 \(n\)이 커지면서 0.05로 수렴한다.

그렇다면 왜 LRT를 쓰는가. 일표본 분산 문제에서는 쓸 이유가 없다. 정확한 검정이 존재하기 때문이다. LRT의 가치는 정확한 검정이 존재하지 않는 복잡한 가설(\(k\)개 분산의 동시 검정 등)로 쉽게 일반화된다는 점에 있다. \(\square\)

Bartlett 검정통계량은 정규 모집단에서 \(k\)개 분산의 동일성에 대한 가능도비 검정이다. LRT는 \(H_0\) 아래의 합동분산 추정값과 개별 분산 추정값들을 비교한다.

여기서 \(s_p^2 = \sum(n_j-1)s_j^2/\sum(n_j-1)\)이다. Bartlett 검정은 카이제곱 근사를 개선하기 위해 보정인자 \(C\)를 적용한다. 보정된 통계량은 \(-2\ln\Lambda / C \sim \chi^2_{k-1}\)이다.

일표본의 경우와의 차이. 연습문제 1에서 일표본 LRT는 \(n(r - 1 - \ln r)\)로 \(r-1\) 항을 포함했다. 다표본에서는 그 항이 사라진다. 이유는

이기 때문이다. 곧 합동분산이 개별 분산의 가중평균이라는 사실이 선형항을 정확히 소거한다.

이는 일표본에서 \(\sigma_0^2\)이 외부에서 주어진 고정값인 반면 다표본에서는 \(S_p^2\)이 자료에서 추정된 값이라는 차이에서 온다. 추정하면서 이미 최적화 조건을 만족시켰으므로 1차 항이 0이 된다. \(\square\)

LRT는 자료가 정규분포를 따른다고 가정한다. 가능도함수 \(L(\sigma^2)\)이 정규 가능도이기 때문이다. 비정규성 아래에서는 가능도가 잘못 설정되고 점근 \(\chi^2\) 분포가 성립하지 않는다. 두꺼운 꼬리가 분산 추정값을 부풀려 검정통계량을 왜곡한다.

더 정확히는 정보행렬 등식이 깨진다. 올바르게 설정된 모형에서는 Fisher 정보의 두 표현(스코어의 분산과 헤시안의 음수 기댓값)이 일치하고, 그 덕분에 \(-2\ln\Lambda\)가 \(\chi^2\)을 따른다. 모형이 잘못 설정되면 두 값이 달라지고 \(-2\ln\Lambda\)는 \(\chi^2\)이 아니라 가중 카이제곱의 혼합을 따른다. 정규 모형에서 그 가중치는 \((\gamma_2+2)/2\)에 비례하며, 이것이 15.1절부터 반복해 온 팽창 인자와 정확히 같은 양이다.

대안.

1. Levene 검정(절대편차 기반, 비정규성에 로버스트)
2. Brown-Forsythe 검정(평균 대신 중앙값 사용)
3. 붓스트랩 LRT(\(\chi^2\) 기준분포를 붓스트랩 기준분포로 대체). LRT 통계량은 그대로 쓰되 기준분포만 자료에서 얻는다. 모형 오설정의 영향을 대부분 제거한다.
4. Fligner-Killeen 검정(순위 기반, 분포무관)
5. 샌드위치(로버스트) 보정. 정보행렬 등식이 깨진 것을 명시적으로 보정하는 방법으로, \(-2\ln\Lambda\)를 추정된 가중치로 나눈다. 준가능도(quasi-likelihood) 이론의 표준 도구이다.

실무적으로는 3번(붓스트랩)이 가장 간단하면서도 효과적이다. 기존 코드에서 \(p\)값 계산 부분만 바꾸면 되기 때문이다. \(\square\)