S₁²/S₂²의 표본분포 (비정규 모집단)¶

(1) 퍼짐의 배율은 첨도 하나로 정해진다. 5.2절의 델타법 결과

에서 척도 \(\sigma^2\)이 사라지고 첨도만 남는다. 정규모집단의 \(\beta_2 = 3\)을 기준으로 삼으면 표준편차의 배율은

이다. 네 모집단에 넣으면 이렇다.

그림에서 볼 것이 이것으로 정해졌다. 균등은 \(F(19,19)\)보다 좁고 높아야 하고, 지수는 두 배 넓어야 하며, 로그정규는 일곱 배 넘게 넓어야 한다. 정규만 검은 곡선에 얹힌다.

로그정규가 일곱 배 넓다는 것은 \(0\)에서 \(6\)까지 그린 가로 범위를 크게 넘어선다는 뜻이기도 하다. 그림에서는 봉우리가 왼쪽 끝으로 눌리고 꼬리가 화면 밖으로 빠져나가는 모습으로 나타날 것이다.

(2) 오른쪽 패널의 방향도 미리 정해진다. 임계값은 \(F(d,d)\)가 정하고 실제 산포는 모집단이 정하는데, 둘 다 \(1/\sqrt n\)로 줄어 비가 바뀌지 않는다(보기 2에서 식으로 쓴다). 그러므로 오류율은 \(0.05\)로 수렴하지 않고 각자의 상수로 간다. 작은 \(n\)에서는 \(F\) 분포 자체가 넓어 왜곡이 가려지므로, 꼬리가 무거운 쪽은 아래에서 위로 올라가고 가벼운 쪽은 위에서 아래로 내려가야 한다.

그려서 확인한다.

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

rng = np.random.default_rng(11)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3.5))

specs = [("Normal", lambda s: rng.normal(size=s), 3.0),
         ("Uniform", lambda s: rng.uniform(size=s), 1.8),
         ("Exp(1)", lambda s: rng.exponential(size=s), 9.0),
         ("LogN(0,1)", lambda s: rng.lognormal(0, 1, size=s), 113.94)]

# 왼쪽: n=20에서 네 모집단의 F 통계량 분포를 겹쳐 그린다.
# 네 경우 모두 귀무가설이 참(같은 모집단에서 두 표본)이므로
# F(19,19)를 따라야 하는데, 정규만 곡선에 얹힌다.
n = 20
for name, rvs, kap in specs:
    F = (rvs((60_000, n)).var(axis=1, ddof=1)
         / rvs((60_000, n)).var(axis=1, ddof=1))
    ax1.hist(F, bins=np.linspace(0, 6, 80), density=True,
             histtype='step', lw=2, label=name)
g = np.linspace(0.01, 6, 400)
ax1.plot(g, stats.f(n - 1, n - 1).pdf(g), "--k", lw=2, label="F(19,19) PDF")
ax1.set_xlim(0, 6)
ax1.set_xlabel("$S_1^2/S_2^2$")
ax1.set_ylabel("Density")
ax1.set_title("Same F curve, four populations ($n_1=n_2=20$)")
ax1.legend(fontsize=8)

# 오른쪽: n을 키우면 오류율이 어떻게 되는가.
# 정규는 0.05에 머물고, 비정규는 각자의 극한값으로 **올라간다**.
ns = [10, 20, 50, 100, 200]
for name, rvs, kap in specs:
    rate = []
    for nn in ns:
        dd = nn - 1
        lo, hi = stats.f(dd, dd).ppf([0.025, 0.975])
        F = (rvs((20_000, nn)).var(axis=1, ddof=1)
             / rvs((20_000, nn)).var(axis=1, ddof=1))
        rate.append(np.mean((F < lo) | (F > hi)))
    ax2.plot(ns, rate, 'o-', lw=2, label=f"{name} (κ={kap:g})")
ax2.axhline(0.05, color='k', ls=':', lw=1)
ax2.set_xscale('log')
ax2.set_xticks(ns)
ax2.set_xticklabels(ns)
ax2.set_xlabel("n (= $n_1$ = $n_2$)")
ax2.set_ylabel("actual type I error")
ax2.set_title("Nominal 5% test: error rate vs sample size")
ax2.legend(fontsize=8)

plt.tight_layout()
plt.show()

왼쪽에서 정규(파란 계단)만 검은 곡선에 얹힌다. 균등은 더 좁고 높으며, 지수와 로그정규는 봉우리가 왼쪽으로 밀리고 오른쪽 꼬리가 두껍다. 오른쪽에서는 정규만 0.05 선에 붙어 있고 나머지는 각자의 극한값으로 향한다.

(1)의 배율이 왼쪽 패널에 그대로 보인다. 정규(계단 하나)만 검은 \(F(19,19)\) 곡선에 얹힌다. 균등은 봉우리가 더 높고 좁은데 표준편차 배율 \(0.632\)가 그 말이다. 지수와 로그정규는 봉우리가 왼쪽으로 밀리고 오른쪽 꼬리가 두꺼워진다. 배율 \(7.5\)인 로그정규는 그릴 수 있는 범위를 벗어나 화면 안에서는 왼쪽 벽에 눌린 모양으로만 보인다.

봉우리가 왼쪽으로 밀리는 까닭도 적어 둘 만하다. \(\log F\)의 퍼짐이 커지면 원래 눈금에서는 위쪽으로 길게 늘어나는데, 전체 넓이가 1로 묶여 있으므로 밀도의 봉우리는 반대로 0 쪽으로 눌린다. 중심이 옮겨 간 것이 아니라 꼬리가 길어진 결과다. 두 모집단이 같은 모양이면 \(E[\log F] = 0\)이 그대로 유지되므로 로그 척도에서는 중심이 움직이지 않는다.

오른쪽 패널의 방향도 예측대로다. 정규만 \(0.05\) 점선에 붙어 있고, 지수와 로그정규는 아래에서 위로 올라가며, 균등은 내려간다. 네 곡선 가운데 수평인 것은 정규 하나뿐이고 나머지는 모두 \(n\)과 함께 명목에서 멀어진다.

이 그림의 한 줄 요약. 네 모집단에 같은 \(F(19,19)\) 자를 들이댔는데 맞는 것은 하나뿐이고, 자를 더 길게(표본을 더 크게) 만들어도 나머지 셋은 맞지 않는다. 오히려 더 어긋난다.

이 그림이 가리는 것도 있다. 왼쪽 패널의 가로 범위가 \(0\)에서 \(6\)으로 잘려 있어 로그정규의 꼬리가 얼마나 긴지 보이지 않는다. \(F(19,19)\)의 \(97.5\%\) 분위수가 \(2.53\)인데 로그정규 자료에서 같은 꼬리확률을 주는 자리는 그보다 한 자릿수 멀다. 그 거리가 바로 오류율 \(0.478\)의 정체이며, 그림이 아니라 보기 2의 표로 읽어야 한다.

(1) 로그 척도로 옮기면 \(\sqrt n\)이 약분된다. \(n_1 = n_2 = n\)이면 기각역이 \(F < F_{0.025}(d,d)\) 또는 \(F > F_{0.975}(d,d)\)이고, 두 경계가 서로 역수이므로 로그 척도에서는 한 줄로 쓰인다.

임계값 \(c\)는 정규 이론이 정한다. \(n\)이 크면 \(\operatorname{Var}(\log F) \to 4/n\)(정규모집단의 델타근사)이므로

이다. 한편 실제 산포는 모집단이 정한다. 5.2절에서 본 델타법 결과

에서 실제 표준편차가 \(s_n = \sqrt{2(\beta_2-1)/n}\)이고, \(\log F\)는 \(n\)이 크면 근사적으로 정규다. 그러므로

가 된다. 분자와 분모에서 \(\sqrt n\)이 함께 약분되어 \(n\)이 사라졌다. 임계값도 \(1/\sqrt n\)로 줄고 실제 산포도 \(1/\sqrt n\)로 줄기 때문이며, 둘의 비는 처음부터 끝까지 \(\sqrt{(\beta_2-1)/2}\)로 고정되어 있다. 표본을 아무리 키워도 이 비가 바뀌지 않는다.

경계는 \(\beta_2 = 3\)이다. 그보다 꼬리가 무거우면 실제 산포가 임계값보다 커져 기각이 쏟아지고, 가벼우면 반대로 거의 기각하지 않는다. 로그정규의 \(\beta_2 = e^4+2e^3+3e^2-3 = 113.94\)에서는 산포가 정규가 믿는 것의 \(7.5\)배다.

(2) 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
B = 20_000

print(f"{'모집단':>14}{'n=20':>9}{'n=50':>9}{'n=200':>9}{'극한(이론)':>12}")
for name, rvs, kap in [("정규", lambda s: rng.normal(size=s), 3.0),
                       ("균등", lambda s: rng.uniform(size=s), 1.8),
                       ("지수", lambda s: rng.exponential(size=s), 9.0),
                       ("로그정규", lambda s: rng.lognormal(0, 1, size=s), 113.94)]:
    row = []
    for n in (20, 50, 200):
        d = n - 1
        lo, hi = stats.f(d, d).ppf([0.025, 0.975])
        F = rvs((B, n)).var(axis=1, ddof=1) / rvs((B, n)).var(axis=1, ddof=1)
        row.append(np.mean((F < lo) | (F > hi)))
    # 극한 오류율 = 2[1 - Phi(1.96 sqrt(2/(kappa-1)))]
    lim = 2 * (1 - stats.norm.cdf(1.96 * np.sqrt(2 / (kap - 1))))
    print(f"{name:>14}{row[0]:>9.3f}{row[1]:>9.3f}{row[2]:>9.3f}{lim:>12.3f}")

출력:

           모집단     n=20     n=50    n=200      극한(이론)
            정규    0.052    0.051    0.051       0.050
            균등    0.005    0.002    0.002       0.002
            지수    0.264    0.298    0.320       0.327
          로그정규    0.478    0.562    0.642       0.794

지수모집단이 0.264 → 0.298 → 0.320으로 극한값 0.327에 다가간다. 이론과 모의실험이 맞는다. 로그정규는 수렴이 느려 \(n = 200\)에서 0.642이지만 방향은 같다.

극한값은 맞고, 거기 닿는 속도가 모집단마다 다르다.

정규 줄은 \(n\)과 무관하게 \(0.05\)다. 극한이 아니라 등식이기 때문이다. 세 줄이 \(0.052,\ 0.051,\ 0.051\)인데 반복 2만 회에서 몬테카를로 오차가 \(0.0015\)이므로 모두 \(0.05\)와 구별되지 않는다.

지수 줄은 \(0.264 \to 0.298 \to 0.320\)으로 극한 \(0.327\)에 또박또박 다가간다. 아래에서 올라가는 것이 핵심이다. 작은 \(n\)에서는 로그 임계값 \(c = \log F_{0.975}(d,d)\)가 자신의 점근값 \(3.92/\sqrt n\)보다 커서(\(n = 20\)에서 \(0.9268\) 대 \(0.8766\)) 기각역이 좁고, 그 완충이 \(n\)과 함께 사라진다. \(n\)이 작을 때 덜 나쁜 것은 \(F\) 분포 자체가 넓어 왜곡을 가려 주기 때문이고, 표본을 키우면 그 가림막이 걷힌다.

로그정규 줄만 \(n = 200\)에서도 \(0.152\) 모자란다. 수렴이 느린 까닭은 델타근사가 선 자리에 있다. \(\operatorname{Var}(\log S^2) \approx (\beta_2-1)/n\)은 \(n\)이 \(\beta_2\)보다 훨씬 커야 쓸 만한데, 로그정규의 \(\beta_2\)가 \(114\)라 \(n = 200\)에서도 \(n/\beta_2 = 1.8\)에 불과하다. 점근식이 뜻을 가지려면 표본이 수천 개는 되어야 한다. 그래도 방향은 분명하다. \(0.478 \to 0.562 \to 0.642\)로 표본을 키울수록 나빠진다.

이 표가 5장에서 가장 역설적인 자리다. 다른 어디서나 표본이 크면 근사가 좋아진다. 여기서는 반대다. 자료를 열 배 모으면 명목 \(5\%\) 검정의 실제 오류율이 지수모집단에서 \(0.264\)에서 \(0.320\)으로 올라간다. 까닭은 중심극한정리가 1차 적률에만 선물을 주기 때문이다. \(F\)는 분자와 분모가 둘 다 2차 적률이라 그 선물을 받지 못하고, 대신 4차 적률에 위아래로 두 번 노출된다.

균등 줄은 반대 방향의 같은 고장이다. \(\beta_2 = 1.8 < 3\)이라 실제 산포가 정규가 믿는 것의 \(0.632\)배뿐이고, 오류율이 \(0.002\)로 내려앉는다. 안전해 보이지만 검정력을 통째로 버린 것이며, 첨도가 3에서 벗어나면 어느 쪽이든 명목이 지켜지지 않는다는 것이 요점이다.

(1) 두 사건이 같은 사건이다. 분산비의 명목 95% 구간은 정리 1을 뒤집어 만든다. \(\sigma_1^2/\sigma_2^2 = \theta\)라 하면 \((S_1^2/S_2^2)/\theta \sim F_{d_1,d_2}\)이므로

이고, 가운데를 \(\theta\)에 대해 풀면 구간

가 나온다. 참값 \(\theta = 1\)이 이 구간에 들어갈 조건은

인데, 이것은 명목 5% 검정이 기각하지 않을 조건과 글자 그대로 같다. 그러므로

가 모집단이 무엇이든 정확히 성립한다. 어림이 아니라 같은 사건을 반대로 센 것이다.

예측은 보기 2의 \(n = 20\) 열을 1에서 빼면 된다.

두 보기가 난수를 따로 뽑으므로 모의값끼리는 몬테카를로 오차만큼 어긋날 수 있다. 반복 2만 회에서 \(p \approx 0.5\)의 오차가 \(\sqrt{0.5 \times 0.5/20000} = 0.0035\)다.

(2) 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
B, n = 20_000, 20
d = n - 1
lo, hi = stats.f(d, d).ppf([0.025, 0.975])

print("명목 95% 구간, n1 = n2 = 20 (참 분산비 = 1)")
for name, rvs in [("정규", lambda s: rng.normal(size=s)),
                  ("균등", lambda s: rng.uniform(size=s)),
                  ("지수", lambda s: rng.exponential(size=s)),
                  ("로그정규", lambda s: rng.lognormal(0, 1, size=s))]:
    F = rvs((B, n)).var(axis=1, ddof=1) / rvs((B, n)).var(axis=1, ddof=1)
    print(f"  {name:>8}: 포함률 {np.mean((F / hi <= 1) & (1 <= F / lo)):.3f}")

출력:

명목 95% 구간, n1 = n2 = 20 (참 분산비 = 1)
        정규: 포함률 0.948
        균등: 포함률 0.994
        지수: 포함률 0.731
      로그정규: 포함률 0.515

로그정규 자료에서 명목 95% 구간이 참값을 절반만 담는다. 신뢰구간이라는 이름이 무의미해지는 수준이다.

네 줄이 예측과 몬테카를로 오차 안에서 맞는다.

가장 크게 어긋난 로그정규의 \(0.007\)이 몬테카를로 오차 \(0.0035\)의 두 배다. 같은 사건을 두 번 센 것이니 이보다 더 맞을 수도 더 틀릴 수도 없다.

로그정규 줄이 이 절의 결론이다. 명목 \(95\%\)를 내건 구간이 참값을 절반만 담는다. 동전을 던져 맞히는 것보다 조금 나은 수준이고, "신뢰구간"이라는 이름이 뜻을 잃는다. 지수모집단에서도 \(0.731\)로 명목에서 \(22\)퍼센트포인트 떨어져 있다.

균등모집단의 \(0.994\)도 고장이다. 방향만 반대일 뿐이다. 꼬리가 가벼워 \(S^2\)이 정규가 믿는 것보다 덜 흔들리므로 구간이 필요 이상으로 넓고, 명목 \(95\%\)를 내걸고 사실상 \(99.4\%\) 구간을 보고한다. 안전하지만 아무것도 배제하지 못한다.

그래서 어긋남의 부호를 정하는 것은 첨도 하나다. \(\beta_2 > 3\)이면 구간이 좁아 위험한 쪽으로, \(\beta_2 < 3\)이면 넓어 쓸모없는 쪽으로 틀린다. 정규분포의 \(\beta_2 = 3\)이 그 경계이고, 그 경계에서만 명목이 지켜진다. \(F\) 구간은 정규성을 가정하는 것이 아니라 첨도가 정확히 3임을 가정한다.

그리고 이 모든 것이 \(n\)을 키워도 나아지지 않는다. 보기 2가 보였듯 \(n = 200\)에서 지수가 \(0.320\), 로그정규가 \(0.642\)로 오류율이 더 올라가므로 포함률은 \(0.680\)과 \(0.358\)로 더 내려간다. 자료를 열 배 모으면 구간이 더 못 미더워진다.

지수분포의 첨도는 \(\beta_2 = 9\)이므로

이고 극한 오류율은

이다.

모의실험에서 \(n = 200\)일 때 0.320이었으므로 방향과 크기가 맞는다. 첨도가 9라는 사실 하나만으로 "명목 5% 검정이 실제로는 33%"라는 결론이 나온다는 점이 이 공식의 쓸모다.

\(g(t) = \log t\)에 델타법을 적용한다. \(g'(t) = 1/t\)이므로 \(S^2\)이 \(\sigma^2\) 근처에 모여 있을 때

이다. 5.6절에서 \(\text{Var}(S^2) = \frac1n\left(\beta_2 - \frac{n-3}{n-1}\right)\sigma^4 \approx \frac{(\beta_2-1)\sigma^4}{n}\)이므로

이다. \(\square\)

\(\sigma^4\)이 약분된다는 점이 요점이다. 로그를 씌우면 척도가 사라지고 첨도만 남는다. 그래서 \(\log F\)의 분포가 \(\sigma\)에 의존하지 않고 첨도에만 의존하며, 검정의 수준이 어긋나는 정도가 첨도만으로 결정된다.

이 계산은 왜 분산 관련 추론에서 첨도가 그토록 중요한지를 설명한다. 평균 추론에서는 2차 적률(분산)이 필요하고, 분산 추론에서는 4차 적률(첨도)이 필요하다. 한 단계 위의 적률이 필요해지는 것이 근본적인 차이다.

로그 척도에서 보면 명확하다. \(\log F = \log S_1^2 - \log S_2^2\)이고 두 항이 독립이므로

분산은 차를 잡아도 더해진다. 상쇄되는 것은 기댓값(편향)이고 분산은 누적된다.

구체적으로, 두 모집단이 같은 모양이면 \(E[\log S_1^2] = E[\log S_2^2]\)이므로 \(E[\log F] \approx 0\)이 되어 중심은 제자리에 있다. 실제로 보기 1의 왼쪽 그림에서 네 분포의 중심이 크게 다르지 않다. 어긋나는 것은 폭이다.

같은 구조가 두 평균 차에도 있다. \(\bar X_1 - \bar X_2\)에서 두 모집단의 치우침이 같으면 차의 분포는 대칭에 가까워지지만(편향 상쇄), 분산은 각자의 분산이 더해진다. 차를 잡는 것은 위치의 오류를 없애 주지만 산포의 오류는 없애 주지 못한다.

1종오류만 보면 안전하다. 실제로 등분산인 두 집단을 "다르다"고 잘못 판정할 확률이 0.5%에 그친다.

그러나 대가가 있다. 검정의 수준이 실제로 0.005라는 것은 기각역이 명목보다 훨씬 좁다는 뜻이고, 그러면 분산이 진짜 다를 때도 잡아내지 못한다. 앞 쪽에서 정규모집단·\(n=20\)에서 분산비 2.25의 검정력이 0.402였는데, 균등모집단에서는 수준이 0.005로 내려간 만큼 검정력도 크게 떨어진다.

더 근본적인 문제는 검정이 자기가 무엇을 하는지 모른다는 점이다. 우리는 5% 검정을 한다고 생각했는데 실제로는 0.5% 검정을 하고 있었다. 보고한 \(p\)-값도 참값이 아니다. 방향이 안전한 쪽이라는 것은 우연이고, 같은 절차가 지수 자료에서는 위험한 쪽으로 어긋난다.

명목 수준이 지켜지지 않는다는 것 자체가 문제이며, 어느 방향이든 절차를 신뢰할 수 없다는 뜻이다.

절차. 두 표본을 합쳐 표본첨도 \(\hat\beta_2\)를 계산하고 극한 오류율 공식에 넣어 본다.

이 값이 0.05에서 크게 벗어나면 \(F\) 검정을 쓰지 않는다. \(\hat\beta_2\)가 4를 넘으면 이미 실제 오류율이 0.1에 가까워진다.

한계가 세 가지 있다.

첫째, 표본첨도가 매우 불안정하다. 4차 적률의 추정이라 표본크기에 따라 크게 흔들린다. \(n = 20\)짜리 자료에서 \(\hat\beta_2\)를 믿기 어렵다.

둘째, 표본첨도는 위쪽으로 편향이 제한된다. 2장에서 본 대로 표본크기 \(n\)으로 계산한 표본첨도는 \(n-2+1/(n-1)\)을 넘을 수 없다. \(n = 20\)이면 약 18.05인데, 로그정규의 참 첨도 113.9를 잡아낼 방법이 아예 없다. 가장 위험한 경우를 가장 못 알아본다.

셋째, 첨도가 유한하지 않을 수도 있다. \(t_4\) 이하의 꼬리에서는 \(\beta_2\)가 존재하지 않으므로 공식 자체가 무의미하다.

그래서 실무의 결론은 "첨도를 재서 판단하라"가 아니라 "애초에 첨도에 의존하지 않는 방법을 쓰라"가 된다. 다음 쪽에서 그 방법들을 본다.

검정의 극한 오류율이 \(2[1-\Phi(z_{1-\alpha/2}\sqrt{2/(\beta_2-1)})]\)이므로 포함률은 그 여사건이다.

로그정규에서 \(\beta_2 = 113.9\)이므로

즉 극한 포함률이 0.206이다. 보기 3의 \(n = 20\)에서 얻은 0.515보다 훨씬 낮다.

왜 차이가 나는가. 로그정규의 수렴이 극도로 느리기 때문이다. 첨도가 113.9라는 것은 4차 적률이 매우 무겁다는 뜻이고, 그런 분포에서 \(\log S^2\)의 정규근사가 통하려면 \(n\)이 수천은 되어야 한다. 보기 2에서 오류율이 \(n = 200\)에서 0.642(포함률 0.358)로 아직 극한 0.794(포함률 0.206)에 못 미치는 것과 같은 이야기다.

이 사실이 상황을 더 나쁘게 만든다. 유한표본에서는 극한값만큼 나쁘지 않지만, 표본을 늘릴수록 극한값을 향해 계속 나빠진다. 분석자에게는 최악의 조합이다. 자료를 더 모으면 결론이 더 신뢰할 만해진다는 통상의 기대가 여기서는 성립하지 않는다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
print(f"{'n':>6}{'참':>7}{'표본 beta2 평균':>16}{'표준편차':>10}{'2.5~97.5%':>22}")
for n in (20, 50, 200, 1000):
    b2 = stats.kurtosis(rng.exponential(1, (20_000, n)), axis=1, fisher=False)
    q = np.percentile(b2, [2.5, 97.5])
    print(f"{n:>6}{9.0:>7.1f}{b2.mean():>16.3f}{b2.std():>10.3f}"
          f"   [{q[0]:.2f}, {q[1]:.2f}]")

출력:

     n      참     표본 beta2 평균      표준편차             2.5~97.5%
    20    9.0           4.478     2.383   [1.80, 10.94]
    50    9.0           6.136     3.392   [2.53, 15.48]
   200    9.0           7.932     3.694   [4.08, 17.36]
  1000    9.0           8.747     2.476   [5.83, 14.81]

\(n=20\)에서 표본첨도의 평균이 \(4.48\)이다. 참값의 절반이다. 게다가 \(95\%\) 범위가 \([1.80,\ 10.94]\)로, 정규분포의 \(3\)은 물론이고 균등분포의 \(1.8\)까지 포함한다.

따라서 "첨도를 재서 판단하라"가 통하지 않는다. \(n=20\)인 지수 자료에서 표본첨도 \(3.2\)를 보고 "정규와 비슷하니 \(F\) 검정을 써도 되겠다"고 판단하는 일이 충분히 일어난다. 실제 오류율은 \(26\%\)인데도 말이다.

왜 이렇게 나쁜가 — 첨도는 4차 적률이고, 그 추정은 8차 적률에 기댄다. 2장에서 본 대로 표본첨도에는 상한 \(n-2+1/(n-1)\)도 있어서, \(n=20\)이면 아무리 극단적인 표본이라도 \(18.05\)를 넘을 수 없다. 꼬리가 두꺼운 모집단일수록 추정이 아래쪽으로 치우치며, 위 표의 \(4.48\)이 정확히 그 편향이다.

\(n=1000\)에서도 \(8.75\)로 아직 못 미친다. 범위도 \([5.83,\ 14.81]\)로 넓다. 표본을 아무리 키워도 "첨도를 알고 있다"고 말하기 어렵다.

실무적 결론. 첨도 진단은 매우 큰 값이 나왔을 때만 정보가 된다. 표본첨도가 \(15\)라면 꼬리가 두껍다는 강한 증거지만, \(3\) 근처라고 해서 정규와 비슷하다는 뜻은 아니다. 한쪽으로만 읽을 수 있는 진단이며, 그래서 연습문제 5의 절차는 "첨도가 작으니 안전"이라는 결론을 내려서는 안 된다.

순열검정의 귀무가설은 "두 분산이 같다"가 아니다. 집단 표지를 섞는 것이 타당하려면 두 집단의 분포 전체가 같아야 한다. 즉 귀무가설은

이다. 이것은 \(F\) 검정이 묻는 \(H_0: \sigma_1^2 = \sigma_2^2\)보다 강한 가설이다.

import numpy as np

rng = np.random.default_rng(0)
n1 = n2 = 25
REP, B = 1000, 299

def ptest(a, b, center):
    if center:                       # 각 집단을 중앙값으로 중심화한 뒤 섞는다
        a, b = a - np.median(a), b - np.median(b)
    obs = np.log(a.var(ddof=1) / b.var(ddof=1))
    pool = np.concatenate([a, b])
    cnt = 0
    for _ in range(B):
        p = rng.permutation(pool)
        r = np.log(p[:n1].var(ddof=1) / p[n1:].var(ddof=1))
        cnt += abs(r) >= abs(obs)
    return (cnt + 1) / (B + 1)

res = {k: 0 for k in ('same_raw', 'shift_raw', 'shift_cent')}
for _ in range(REP):
    x, y = rng.exponential(1, n1), rng.exponential(1, n2)
    res['same_raw']   += ptest(x, y,       False) < 0.05
    res['shift_raw']  += ptest(x, y + 3.0, False) < 0.05
    res['shift_cent'] += ptest(x, y + 3.0, True)  < 0.05

print(f"  분포 동일, 원자료 순열       {res['same_raw']/REP:.3f}")
print(f"  평균 이동, 원자료 순열       {res['shift_raw']/REP:.3f}")
print(f"  평균 이동, 중앙값 중심화 후  {res['shift_cent']/REP:.3f}")

출력:

  분포 동일, 원자료 순열       0.052
  평균 이동, 원자료 순열       0.406
  평균 이동, 중앙값 중심화 후  0.052

분포가 같을 때는 \(0.052\)로 정확하다. \(F\) 검정이 \(\text{Exp}(1)\)에서 \(0.26\)까지 치솟았던 것과 대조적이다. 순열검정은 정규성을 전혀 쓰지 않으므로 첨도에 영향받지 않는다.

그런데 평균만 \(3\) 이동시키자 \(0.406\)으로 무너진다. 두 모분산은 여전히 완전히 같은데도 열 번 중 네 번 기각한다. \(F\) 검정보다도 나쁘다.

원인을 정확히 짚어야 한다. 관측된 통계량 \(S_1^2/S_2^2\)은 각 집단의 중심을 뺀 편차로 계산되므로 위치 이동에 영향받지 않는다. 무너지는 것은 귀무분포 쪽이다. 원자료를 한데 모아 섞으면 평균이 다른 두 덩어리가 뒤섞여, 섞인 표본의 분산이 집단 간 평균 차이만큼 부풀려진다. 그 결과 귀무분포가 실제보다 오른쪽으로 밀려나고, 관측값이 엉뚱하게 극단적으로 보인다.

고치는 방법은 간단하다. 섞기 전에 각 집단을 중심화한다. 위 표의 셋째 줄이 그것이며, \(0.052\)로 완전히 회복된다.

일반 교훈 — 순열검정의 귀무가설은 "교환가능성"이다. 표지를 바꿔도 분포가 같아야 섞는 것이 타당한데, 평균이 다르면 그 전제가 깨진다. 분산만 비교하고 싶다면 비교하지 않을 측면(위치)을 먼저 제거한 뒤 섞어야 한다.

"순열검정은 가정이 없다"는 말이 틀린 이유가 이것이다. 가정이 정규성에서 교환가능성으로 옮겨갔을 뿐이며, 위에서 보듯 그 가정이 깨지면 결과가 \(F\) 검정보다 나쁠 수도 있다. 다만 교환가능성은 정규성보다 다루기 쉽고, 중심화 같은 간단한 처리로 복구할 수 있다는 점이 실질적인 이득이다. 17장에서 이 논점을 자세히 다룬다.

원리 — 산포 문제를 위치 문제로 바꾼다. 각 관측을 자기 집단 중심으로부터의 절대편차로 바꾼다.

그러면 "집단 간 산포가 다른가"라는 질문이 "\(z\)의 평균이 집단마다 다른가"로 바뀌고, 이것은 통상의 일원배치 분산분석으로 답할 수 있다.

왜 첨도에 덜 민감한가. \(F\) 검정이 무너진 이유는 \(S^2\)의 표본분포가 4차 적률에 기대기 때문이었다(\(\operatorname{Var}(S^2) \approx \sigma^4(\beta_2-1)/n\)). 레빈 검정은 \(z\)의 평균을 비교하므로, 필요한 것이 \(z\)의 2차 적률뿐이다. 두 차수 낮은 적률만 쓰는 것이 강건성의 근원이다.

게다가 \(z\)의 평균에는 중심극한정리가 듣는다. 5.6절에서 본 대로 \(\bar X\)는 모집단이 무엇이든 정규로 수렴하지만 \(S^2\)은 훨씬 느리다. 레빈은 그 빠른 쪽에 문제를 얹은 셈이다.

브라운–포사이드는 중심을 바꾼다. 평균 대신 중앙값을 쓴다.

왜 더 나은가. 모집단이 심하게 치우쳐 있으면 평균이 꼬리에 끌려가고, 그러면 \(\lvert x - \bar x\rvert\)가 분산이 아니라 왜도를 반영하게 된다. 중앙값은 꼬리에 끌리지 않으므로 이 오염이 줄어든다. 지수분포처럼 왜도가 큰 모집단에서 차이가 뚜렷하다.

레빈이라고만 쓰면 어느 쪽인지 알 수 없다. 중심을 평균으로 잡은 원래 형태는 치우친 자료에서 함께 무너진다. \(\text{Exp}(1)\)에서 두 집단에 각 \(50\)개씩 뽑아 재면 평균 중심 레빈의 제1종 오류율이 \(0.139\)인 반면 중앙값 중심, 곧 브라운–포사이드는 \(0.049\)다. scipy.stats.levene의 기본값이 center='median'이라 그냥 쓰면 브라운–포사이드가 되며, 이 책에서 특별히 밝히지 않고 "레빈"이라 적은 곳은 그 기본값을 뜻한다.

보정의 발상. 이 쪽의 핵심 결과는

였다(연습문제 2). 정규라면 \(\beta_2 = 3\)이라 \(2/n\)인데, 일반 모집단에서는 \((\beta_2-1)/n\)이다. 따라서

이고, 이것을 표준화해

로 검정하면 된다. \(\widehat\beta_2\)는 두 표본을 합쳐 추정한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n1 = n2 = 50
REP = 20_000
rej_F = rej_adj = 0
for _ in range(REP):
    x, y = rng.exponential(1, n1), rng.exponential(1, n2)
    F = x.var(ddof=1) / y.var(ddof=1)
    # (1) 보통의 F 검정
    d1, d2 = n1 - 1, n2 - 1
    rej_F += (F > stats.f.ppf(0.975, d1, d2)) or (F < stats.f.ppf(0.025, d1, d2))
    # (2) 첨도 보정
    pooled = np.concatenate([x - x.mean(), y - y.mean()])
    b2 = stats.kurtosis(pooled, fisher=False)
    se = np.sqrt(max(b2 - 1, 0.1) * (1/n1 + 1/n2))
    rej_adj += abs(np.log(F)) / se > 1.959964

print(f"  n={n1}, Exp(1), 명목 0.05")
print(f"    보통의 F 검정   기각률 {rej_F/REP:.4f}")
print(f"    첨도 보정 검정  기각률 {rej_adj/REP:.4f}")

출력:

  n=50, Exp(1), 명목 0.05
    보통의 F 검정   기각률 0.2963
    첨도 보정 검정  기각률 0.0708

보정이 대부분을 고친다(\(0.296 \to 0.071\)). 명목값 \(0.05\)까지는 못 가지만, 오류율 초과분이 \(0.246\)에서 \(0.021\)로 \(92\%\) 줄었다. 예상보다 잘 듣는다.

그래도 남는 \(0.071\)의 원인은 연습문제 7이 설명한다. \(n=50\)에서 표본첨도의 평균이 \(6.14\)로 참값 \(9\)에 못 미친다. \(\widehat\beta_2\)를 과소추정하면 표준오차가 작아지고, 그러면 기각을 조금 더 많이 하게 된다. 보정에 쓰는 재료 자체가 아래로 치우쳐 있는 것이 남은 오차의 원인이다.

실용성 평가. 보정은 쓸 만하다. \(F\) 검정을 그대로 쓰는 것보다 훨씬 낫고, 계산도 간단하다. 큰 표본에서는 \(\widehat\beta_2\)가 개선되므로 더 나아진다.

그럼에도 레빈이나 순열이 낫다고 보는 이유가 있다.

• 보정은 추정한 모수에 기댄다. \(\widehat\beta_2\)의 오차가 그대로 새 문제로 들어오며, 표본이 작을수록(연습문제 7의 \(n=20\)에서 평균 \(4.48\)) 보정이 오히려 어긋날 수 있다.
• 점근 근사를 한 겹 더 쌓는다. \(\log F\)의 정규근사는 \(n\)이 커야 타당한데, 그 \(n\)이 얼마여야 하는지는 또 첨도에 달려 있다.
• 대안은 애초에 첨도를 알 필요가 없는 통계량으로 바꾼다. 모르는 것을 추정해 메우는 대신, 모르는 것에 의존하지 않도록 문제를 옮기는 편이 안전하다.

모르는 모수를 추정해 끼워 넣는 것과, 그 모수가 등장하지 않는 방법으로 바꾸는 것 — 이 둘의 대비는 강건 통계 전반에서 반복되는 선택이다. 여기서는 후자가 조금 더 낫지만, 전자도 충분히 실용적이라는 것이 이 모의실험의 결론이다. \(\square\)