X̄₁ - X̄₂의 표본분포 (비정규 모집단)¶

(1) 대칭은 점근적 성질이 아니다. \(X_1,\ldots,X_n\)과 \(Y_1,\ldots,Y_n\)을 같은 분포에서 같은 개수로 뽑았으므로 두 표본의 이름을 맞바꾸어도 자료 전체의 결합분포가 변하지 않는다. 이름을 바꾸면

이다. 분모와 자유도는 두 집단을 대칭으로 쓰고 분자만 부호가 뒤집히므로 \(T \longmapsto -T\)다. 따라서 \(T\)와 \(-T\)가 같은 분포를 따르고, 임계값도 \(\nu\)만으로 정해지므로

가 모든 \(n\)과 모든 모집단에서 정확히 성립한다. 왜도가 \(6.18\)이든 \(0\)이든 상관없다.

이것을 누적률로 다시 보면 본문의 식이 된다. \(\sigma_1 = \sigma_2\), \(n_1 = n_2\)에서

이고, 홀수 차수 누적률이 모두 같은 이유로 지워진다. 에지워스 전개의 맨 앞 \(O(n^{-1/2})\) 항이 왜도에 비례하므로 그 항이 통째로 사라진다. 남는 것은 첨도가 만드는 \(O(n^{-1})\) 항이다.

(2) 합계는 \(0.05\)보다 작아야 한다. 꼬리가 무거운 모집단에서는 어쩌다 한 집단에 큰 값이 들어와 \(S_i^2\)이 몇 배로 뛴다. 그러면 Welch 통계량의 분모가 커져 \(|T|\)가 \(0\) 쪽으로 눌리고, 동시에 \(\nu\)가 작아져 임계값이 커진다. 두 효과가 모두 기각을 줄이는 쪽이다. 따라서

를 예측한다. 첨도가 클수록 어긋남이 크므로 정규 \(0\), 지수 \(6\), 로그정규 \(110.9\) 순으로 벌어져야 한다.

이제 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000
alpha = 0.05


def draw(pop, size):
    if pop == "Exp":
        return rng.exponential(1.0, size=size)          # 왜도 2
    if pop == "LogN":
        return rng.lognormal(0.0, 1.0, size=size)       # 왜도 6.18
    return rng.normal(0.0, 1.0, size=size)              # 왜도 0


def welch_reject(pop, n):
    """두 집단을 같은 모집단에서 뽑아 H0가 참인 상황을 만든다."""
    x1, x2 = draw(pop, (B, n)), draw(pop, (B, n))
    a = x1.var(axis=1, ddof=1) / n
    b = x2.var(axis=1, ddof=1) / n
    t = (x1.mean(axis=1) - x2.mean(axis=1)) / np.sqrt(a + b)
    nu = (a + b) ** 2 / (a ** 2 / (n - 1) + b ** 2 / (n - 1))
    crit = stats.t(nu).ppf(1 - alpha / 2)
    return np.mean(t < -crit), np.mean(t > crit)


print("모집단     n    왼쪽꼬리  오른쪽꼬리   합계")
for pop in ("Normal", "Exp", "LogN"):
    for n in (10, 30, 100):
        left, right = welch_reject(pop, n)
        print(f"{pop:<8} {n:>4}     {left:.3f}     {right:.3f}    {left + right:.3f}")

출력:

모집단     n    왼쪽꼬리  오른쪽꼬리   합계
Normal     10     0.025     0.024    0.049
Normal     30     0.025     0.025    0.050
Normal    100     0.024     0.025    0.049
Exp        10     0.019     0.019    0.037
Exp        30     0.024     0.022    0.046
Exp       100     0.025     0.023    0.048
LogN       10     0.013     0.014    0.027
LogN       30     0.020     0.020    0.040
LogN      100     0.022     0.023    0.045

세 가지를 읽는다.

첫째, 두 꼬리가 같다. 어느 줄에서도 왼쪽과 오른쪽이 0.001 안쪽으로 붙어 있다. 이론에서 본 왜도의 상쇄가 그대로 나타난 것이다. 왜도 6.18짜리 모집단에서도 차의 분포는 대칭이다.

둘째, 어긋남이 보수적인 방향이다. 로그정규 \(n = 10\)에서 0.027로 명목값의 절반 남짓이다. 위험한 쪽이 아니라 검정력을 잃는 쪽이다.

셋째, \(n\)을 키우면 사라진다. 0.05와의 차이가 로그정규에서 0.023 → 0.010 → 0.005로 줄어든다. \(n\)이 3배가 될 때마다 절반 남짓으로 줄어드는 셈이어서, \(O(n^{-1})\)이 예측하는 3분의 1보다는 조금 느리다(연습문제 4).

그러나 \(n = 100\)에서도 0.045다. 완전히 회복되지는 않았다. 첨도 111짜리 모집단에서는 수백 개가 필요하다.

(1)이 그대로 확인된다. 아홉 줄 어디에서도 두 꼬리가 \(0.002\) 넘게 벌어지지 않는다. 반복 5만 회에서 두 꼬리 차의 몬테카를로 오차는

이므로 모든 줄이 \(2\) 오차 안이다. 왜도 \(6.18\)짜리 모집단에서 \(n = 10\)으로도 \(0.013\)과 \(0.014\)다. 대칭은 표본크기가 벌어 주는 것이 아니라 설계가 거저 주는 것이다.

(2)의 방향도 맞는다. 합계가 정규 \(0.049\), 지수 \(0.037\), 로그정규 \(0.027\)로 첨도 순서대로 아래로 내려간다. 모두 \(0.05\)보다 작다.

속도는 따로 재야 한다. 표의 합계는 몬테카를로 오차가 \(0.001\)인데 어긋남 자체가 \(n = 100\)에서 \(0.002 \sim 0.005\)라 자릿수가 모자란다. 반복을 50만 번으로 올리고 \(n = 300\)을 더해 다시 잰다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
alpha, B = 0.05, 500_000

def gap(pop, n):
    """0.05 와의 차이를 몬테카를로 오차 아래로 끌어내리려 반복을 50만 번으로 올린다."""
    f = (rng.exponential if pop == "Exp" else rng.lognormal)
    arg = (1.0,) if pop == "Exp" else (0.0, 1.0)
    x1, x2 = f(*arg, size=(B, n)), f(*arg, size=(B, n))
    a, b = x1.var(axis=1, ddof=1) / n, x2.var(axis=1, ddof=1) / n
    t = (x1.mean(axis=1) - x2.mean(axis=1)) / np.sqrt(a + b)
    nu = (a + b) ** 2 / (a ** 2 / (n - 1) + b ** 2 / (n - 1))
    p = np.mean(np.abs(t) > stats.t(nu).ppf(1 - alpha / 2))
    return p, np.sqrt(p * (1 - p) / B)

print(f"{'모집단':>6}{'n':>6}{'실제':>9}{'MC오차':>9}{'0.05-실제':>11}{'국소 기울기':>12}")
for pop in ("Exp", "LogN"):
    prev = None
    for n in (10, 30, 100, 300):
        p, se = gap(pop, n)
        g = 0.05 - p
        slope = "" if prev is None else f"{np.log(prev[1] / g) / np.log(n / prev[0]):.2f}"
        print(f"{pop:>6}{n:>6}{p:>9.4f}{se:>9.4f}{g:>11.4f}{slope:>12}".rstrip())
        prev = (n, g)

출력:

   모집단     n       실제     MC오차    0.05-실제      국소 기울기
   Exp    10   0.0376   0.0003     0.0124
   Exp    30   0.0467   0.0003     0.0033        1.21
   Exp   100   0.0493   0.0003     0.0007        1.27
   Exp   300   0.0499   0.0003     0.0001        2.28
  LogN    10   0.0275   0.0002     0.0225
  LogN    30   0.0394   0.0003     0.0106        0.68
  LogN   100   0.0451   0.0003     0.0049        0.64
  LogN   300   0.0477   0.0003     0.0023        0.68

"국소 기울기"는 \(\log(\Delta_{\text{앞}}/\Delta_{\text{뒤}})\big/\log(n_{\text{뒤}}/n_{\text{앞}})\)이므로 \(O(n^{-1})\)이면 \(1\), \(O(n^{-1/2})\)이면 \(0.5\)가 나와야 한다.

지수는 예측대로 \(O(n^{-1})\)이다. 기울기가 \(1.21\)과 \(1.27\)로 \(1\) 둘레에 있고, \(n = 300\)의 \(2.28\)은 어긋남 \(0.0001\)이 몬테카를로 오차 \(0.0003\)보다 작아 읽을 수 없는 값이다. 어긋남이 \(0.0124 \to 0.0001\)로 두 자릿수 줄었다.

로그정규는 기울기가 \(0.64 \sim 0.68\)에 머문다. \(O(n^{-1})\)보다 뚜렷이 느리다. 왜도 항이 지워진 것은 맞지만 첨도가 만드는 항이 아직 점근식에 들어가지 못했기 때문이다. 전개의 크기를 재는 양이 초과첨도를 \(n\)으로 나눈 \(\gamma_2/n = 110.9/n\)인데, \(n = 100\)에서도 \(1.1\)로 아직 \(1\)을 넘는다. \(O(n^{-1})\)이라는 말이 뜻을 가지려면 이 양이 \(1\)보다 한참 작아야 한다.

그래서 두 결론을 함께 적어야 정직하다. 비정규성의 어긋남은 분명히 사라진다 — 지수에서 \(n = 300\)이면 \(0.0499\)로 명목과 구별되지 않는다. 그러나 빨리 사라진다는 보장은 없다 — 로그정규는 \(n = 300\)에서도 \(0.0477\)로 \(0.0023\) 모자란다. "\(n \ge 30\)이면 중심극한정리"라는 어림은 왜도 \(1\) 안팎의 온건한 모집단을 두고 한 말이며, 왜도 \(6.18\)·초과첨도 \(110.9\)짜리 모집단에는 통하지 않는다(연습문제 4).

(1) 상쇄 조건이 깨졌다. 보기 1에서 두 항이 지워진 것은 \(\sigma_1^3/n_1^2 = \sigma_2^3/n_2^2\)였기 때문이다. 여기서는 \(n_1 = n_2 = n\)이지만 \(\sigma_2 = 2\sigma_1\)이므로

로 지워지기는커녕 일곱 배가 되어 남는다. 두 집단이 같은 모양이라 왜도의 부호도 같은데, 큰 쪽이 세제곱으로 들어가니 뺄셈이 오히려 키운 셈이다. 분산은

이므로 표준화 왜도는

이다. \(\sigma_1\)이 약분되어 사라지고 모집단의 왜도와 \(n\)만 남는다. \(\text{Exp}(1)\)은 \(\gamma_1 = 2\), \(\text{LogNormal}(0,1)\)은 \(\gamma_1 = (e+2)\sqrt{e-1} = 6.1849\)이므로

이고 \(1/\sqrt n\) 속도로만 준다. 모의로 재어 닫힌 꼴이 맞는지 확인해 둔다.

import numpy as np
from scipy import stats

def closed(g1, s1, s2, n1, n2):
    """차의 표준화 왜도. 누적률은 더해지고 상수배에서 세제곱으로 늘어난다."""
    k3 = g1 * (s1 ** 3 / n1 ** 2 - s2 ** 3 / n2 ** 2)
    v = s1 ** 2 / n1 + s2 ** 2 / n2
    return k3 / v ** 1.5

rng = np.random.default_rng(11)
s_exp, g_exp = 1.0, 2.0
s_log, g_log = np.sqrt(np.e ** 2 - np.e), (np.e + 2) * np.sqrt(np.e - 1)

print(f"{'모집단':>6}{'n':>5}{'닫힌 꼴':>10}{'모의':>10}")
for name, s, g, f in (("Exp", s_exp, g_exp, lambda sz: rng.exponential(1.0, sz)),
                      ("LogN", s_log, g_log, lambda sz: rng.lognormal(0.0, 1.0, sz))):
    for n in (10, 30, 100):
        d = f((400_000, n)).mean(1) - 2.0 * f((400_000, n)).mean(1)
        print(f"{name:>6}{n:>5}{closed(g, s, 2 * s, n, n):>10.4f}{stats.skew(d):>10.4f}")

출력:

   모집단    n      닫힌 꼴        모의
   Exp   10   -0.3960   -0.3991
   Exp   30   -0.2286   -0.2239
   Exp  100   -0.1252   -0.1263
  LogN   10   -1.2245   -1.2428
  LogN   30   -0.7070   -0.6955
  LogN  100   -0.3872   -0.3976

여섯 줄이 모두 셋째 자리까지 맞는다. 40만 번에서 왜도 추정의 몬테카를로 오차가 대략 \(\sqrt{6/400000} = 0.0039\)이므로 그 서너 배 안이며, 꼬리가 무거운 분포에서 표본왜도 자체가 느리게 수렴한다는 것을 감안하면 닫힌 꼴이 맞다.

(2) 부호가 방향을 정한다. 차의 왜도가 음수이므로 \(\bar X_1 - \bar X_2\)의 분포는 왼쪽 꼬리가 길고 봉우리는 평균보다 오른쪽에 있다. 그러면 대부분의 표본에서 \(\bar X_1-\bar X_2\)가 참값 \(\delta\)보다 크게 나오고, 그 둘레에 세운 구간도 \(\delta\) 위쪽에 놓이기 쉽다.

왼쪽 꼬리는 길지만 거기에 해당하는 표본에서는 큰 값 하나가 \(S_2^2\)도 함께 키워 구간이 넓어지므로 빗나가기보다 덮어 버린다. 그래서 모자라는 포함률의 거의 전부가 "위쪽으로 빗나감"으로 간다.

이제 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000
alpha = 0.05

# 모집단 2는 모집단 1을 2배로 늘린 것이다. 따라서 mu2 = 2*mu1, sigma2 = 2*sigma1 이고
# 참 차이는 delta = mu1 - mu2 = -mu1 이다. 비정규성과 이분산이 함께 들어 있다.
pops = {"Exp": 1.0, "LogN": np.exp(0.5)}     # 각 모집단의 평균 mu1

print("모집단     n    포함률   왼쪽으로 벗어남  오른쪽으로 벗어남   mean nu")
for pop, mu1 in pops.items():
    for n in (10, 30, 100):
        if pop == "Exp":
            x1 = rng.exponential(1.0, size=(B, n))
            x2 = 2.0 * rng.exponential(1.0, size=(B, n))
        else:
            x1 = rng.lognormal(0.0, 1.0, size=(B, n))
            x2 = 2.0 * rng.lognormal(0.0, 1.0, size=(B, n))
        delta = -mu1
        a = x1.var(axis=1, ddof=1) / n
        b = x2.var(axis=1, ddof=1) / n
        se = np.sqrt(a + b)
        nu = (a + b) ** 2 / (a ** 2 / (n - 1) + b ** 2 / (n - 1))
        crit = stats.t(nu).ppf(1 - alpha / 2)
        d = x1.mean(axis=1) - x2.mean(axis=1)
        lo, hi = d - crit * se, d + crit * se
        cover = np.mean((lo <= delta) & (delta <= hi))
        below = np.mean(hi < delta)      # 구간이 참값보다 아래에 놓인 경우
        above = np.mean(lo > delta)
        print(f"{pop:<8} {n:>4}   {cover:.3f}        {below:.3f}            {above:.3f}       {nu.mean():6.1f}")

출력:

모집단     n    포함률   왼쪽으로 벗어남  오른쪽으로 벗어남   mean nu
Exp        10   0.932        0.007            0.062         13.5
Exp        30   0.942        0.011            0.048         43.7
Exp       100   0.948        0.016            0.036        147.2
LogN       10   0.920        0.003            0.077         13.1
LogN       30   0.931        0.005            0.065         43.4
LogN      100   0.939        0.010            0.051        148.6

여기서는 두 척도가 달라 왜도의 상쇄가 일어나지 않는다. 결과가 보기 1과 두 가지 면에서 다르다.

포함률이 명목값보다 낮다. 로그정규 \(n=10\)에서 0.920으로 3%포인트 모자란다. 보수적인 쪽이 아니라 위험한 쪽이다.

벗어나는 방향이 한쪽으로 몰린다. 두 꼬리가 각각 0.025여야 하는데 0.003과 0.077로 25배 차이가 난다. 치우친 모집단에서 \(\bar X_2\)는 참 평균보다 작게 나오는 일이 많으므로(\(\text{중앙값} < \text{평균}\)) 차 \(\bar X_1 - \bar X_2\)가 참값보다 크게 나오고, 구간이 참값 위쪽에 놓이는 일이 잦다. 명목 95% 구간이지만 한쪽으로만 틀린다.

\(n\)이 커지면 둘 다 고쳐진다. 로그정규에서 포함률이 0.920 → 0.931 → 0.939이고 두 꼬리도 0.010과 0.051로 가까워진다. 다만 그 속도가 \(O(n^{-1/2})\)이라 느리다.

예측이 맞았다. 로그정규 \(n=10\)에서 아래쪽이 \(0.003\), 위쪽이 \(0.077\)로 25배 차이 난다. 명목대로라면 둘 다 \(0.025\)여야 한다. 반복 5만 회에서 \(p \approx 0.077\)의 몬테카를로 오차가 \(0.0012\)이므로 우연이 아니다.

두 가지를 읽는다.

첫째, 보기 1과 어긋남의 방향이 반대다. 균형·동척도에서는 포함률이 명목보다 높은 쪽(검정으로 치면 보수적인 쪽)이었는데, 여기서는 \(0.920\)으로 낮다. 안전한 고장이 아니라 위험한 고장이다. 바뀐 것은 척도 하나뿐이고, 그것이 \(\kappa_3\)을 \(0\)에서 \(-7\gamma_1\sigma_1^3/n^2\)으로 옮겼다.

둘째, 비대칭의 크기가 왜도를 따라간다. 마지막 열을 왜도로 나누면 지수는 \(0.139,\ 0.162,\ 0.160\), 로그정규는 \(0.060,\ 0.085,\ 0.106\)이다. 지수 쪽은 거의 일정해 비대칭이 왜도에 비례하고 따라서 \(1/\sqrt n\)로 준다. 로그정규 쪽은 그 비가 아직 올라가는 중이어서 비대칭이 왜도보다 느리게 준다. 실제로 \(0.074 \to 0.060 \to 0.041\)은 \(n\)을 열 배 키워 절반가량 준 것이고, \(1/\sqrt n\)이 예측하는 \(1/\sqrt{10} = 0.32\)배에 한참 못 미친다.

까닭은 로그정규가 아직 점근 영역에 없기 때문이다. 초과첨도가 \(110.9\)라 전개의 크기를 재는 양 \(\gamma_2/n\)이 \(n = 100\)에서도 \(1.1\)로 \(1\)을 넘는다. 왜도 항만으로 설명되지 않는 몫이 남아 있다는 뜻이고, 그 몫은 보기 1에서 측정한 \(O(n^{-1})\) 항과 같은 자리에서 온다.

포함률은 어느 쪽이든 \(n\)과 함께 회복된다. 로그정규에서 \(0.920 \to 0.931 \to 0.939\)이고 두 꼬리도 \(0.010\)과 \(0.051\)로 가까워진다. 다만 그 속도가 \(O(n^{-1/2})\)이라 균형 설계의 \(O(n^{-1})\)보다 느리고, 척도를 맞추거나 표본크기를 \((\sigma_1/\sigma_2)^{3/2}\) 비로 배정하면(연습문제 6) 왜도 항을 되돌려 지울 수 있다.

(1) 왜도는 닫힌 꼴로 나온다. 연습문제 6에서 유도한

에 \(\sigma_1 = 4\), \(\sigma_2 = 1\)을 넣으면 표준화 왜도가 본문의 표가 된다.

\(64/100\)과 \(1/1600\)을 견주면 첫 항이 천 배 크다. 왜도의 거의 전부가 "작은 표본 + 큰 분산" 쪽 집단에서 온다. 그리고 세 줄이 \(1/\sqrt n\) 속도로 줄어든다. 균형 설계에서 통째로 사라졌던 \(O(n^{-1/2})\) 항이 되살아난 것이다.

부호를 꼬리로 옮기는 데 한 번 뒤집기가 들어간다. 차의 왜도는 양수다. 그런데 기각을 세는 것은 차가 아니라 \(t\) 통계량이고, 스튜던트화는 왜도의 부호를 뒤집는다. 큰 값 하나가 들어오면 분자 \(\bar X_1 - \bar X_2\)가 커지지만 분모 \(\sqrt{S_1^2/n_1+S_2^2/n_2}\)도 함께 커져 비가 도로 눌리고, 거꾸로 그 큰 값이 빠진 표본에서는 분자도 작고 분모는 더 작아 비가 크게 음수로 간다. 그래서

이고, 왼쪽 꼬리의 기각이 늘어야 한다. 이 쪽 연습문제 2의 일표본 표가 같은 현상을 보인다. 오른쪽으로 치우친 모집단에 일표본 \(t\)를 걸면 왼쪽 꼬리가 0.095, 오른쪽이 0.004였다.

(2) 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B = 50_000
alpha = 0.05


def rates(pop, n1, n2, scale1=4.0):
    """모집단 1은 척도가 4배(분산 16배)이고 표본은 적다. 두 평균은 0으로 맞춘다."""
    if pop == "Exp":
        m = 1.0
        x1 = scale1 * (rng.exponential(1.0, size=(B, n1)) - m)
        x2 = rng.exponential(1.0, size=(B, n2)) - m
    else:
        m = np.exp(0.5)
        x1 = scale1 * (rng.lognormal(0.0, 1.0, size=(B, n1)) - m)
        x2 = rng.lognormal(0.0, 1.0, size=(B, n2)) - m
    d = x1.mean(axis=1) - x2.mean(axis=1)
    v1, v2 = x1.var(axis=1, ddof=1), x2.var(axis=1, ddof=1)

    sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2)
    t_pool = d / np.sqrt(sp2 * (1 / n1 + 1 / n2))
    rej_pool = np.mean(np.abs(t_pool) > stats.t(n1 + n2 - 2).ppf(1 - alpha / 2))

    a, b = v1 / n1, v2 / n2
    t_w = d / np.sqrt(a + b)
    nu = (a + b) ** 2 / (a ** 2 / (n1 - 1) + b ** 2 / (n2 - 1))
    crit = stats.t(nu).ppf(1 - alpha / 2)
    return rej_pool, np.mean(np.abs(t_w) > crit), np.mean(t_w < -crit), np.mean(t_w > crit)


print("치우침 + 이분산 + 불균형: sigma1/sigma2 = 4,  n1 : n2 = 1 : 4")
print("모집단    (n1,  n2)   pooled t   Welch t   (Welch 왼쪽 / 오른쪽)")
for pop in ("Exp", "LogN"):
    for n1, n2 in ((10, 40), (30, 120), (100, 400)):
        rp, rw, left, right = rates(pop, n1, n2)
        print(f"{pop:<8} ({n1:>3}, {n2:>3})    {rp:.3f}      {rw:.3f}      {left:.3f} / {right:.3f}")

출력:

치우침 + 이분산 + 불균형: sigma1/sigma2 = 4,  n1 : n2 = 1 : 4
모집단    (n1,  n2)   pooled t   Welch t   (Welch 왼쪽 / 오른쪽)
Exp      ( 10,  40)    0.316      0.100      0.096 / 0.004
Exp      ( 30, 120)    0.293      0.072      0.064 / 0.008
Exp      (100, 400)    0.282      0.058      0.045 / 0.013
LogN     ( 10,  40)    0.339      0.154      0.153 / 0.001
LogN     ( 30, 120)    0.310      0.111      0.109 / 0.002
LogN     (100, 400)    0.294      0.081      0.076 / 0.005

이 표가 이 절의 두 교훈을 한 화면에 담고 있다.

합동 \(t\) 열을 보라. 0.316 → 0.293 → 0.282, 0.339 → 0.310 → 0.294다. 표본을 10배로 늘려도 거의 움직이지 않는다. 앞 쪽에서 본 이분산·불균형의 어긋남이며, 표본크기가 고쳐 주지 않는다.

Welch 열을 보라. 0.100 → 0.072 → 0.058, 0.154 → 0.111 → 0.081이다. 느리지만 또박또박 0.05로 내려간다. Welch는 이분산 문제를 이미 해결했으므로 남은 것은 비정규성뿐이고, 그것은 표본크기가 고쳐 준다. 감소 속도가 이론이 예측한 \(O(n^{-1/2})\)과 맞는다.

꼬리 분해를 보면 어긋남이 한쪽으로 완전히 몰려 있다. 로그정규 \((10,40)\)에서 0.153 대 0.001이다. 사실상 단측검정이 되어 버렸다. 균형 설계에서 두 꼬리가 0.013과 0.014로 맞아떨어지던 보기 1과 견주면 불균형이 무엇을 망가뜨리는지 분명하다.

세 가지가 예측대로다.

첫째, 어긋남이 왼쪽 꼬리에 몰렸다. 로그정규 \((10,40)\)에서 \(0.153\) 대 \(0.001\)이다. 명목 \(0.025\)씩이어야 할 두 꼬리가 150배 차이 난다. 양측검정이 사실상 단측검정이 되었다. 균형 설계의 보기 1에서 두 꼬리가 \(0.013\)과 \(0.014\)로 맞아떨어지던 것과 견주면, 불균형이 지운 것이 무엇인지 분명하다.

둘째, 어긋남의 크기가 왜도를 따라간다. \(0.05\)를 넘는 초과분을 (1)의 왜도로 나누면 이렇다.

로그정규 세 줄의 마지막 열이 \(0.054,\ 0.055,\ 0.051\)로 거의 일정하다. 초과분이 왜도에 비례한다는 뜻이고, 왜도가 \(1/\sqrt n\)으로 주니 초과분도 그 속도로 준다. 비례가 어디까지 가는지 한 칸 더 밀어 본다.

\((300,1200)\)의 왜도는 \(0.3485\)이므로 비례가 이어지면 오류율이 \(0.05 + 0.053 \times 0.3485 = 0.0685\)쯤이어야 한다. 재어 본다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
alpha, B, m = 0.05, 50_000, np.exp(0.5)

# (300, 1200) 은 한 번에 올리면 6천만 개라 메모리가 모자란다. 200개씩 끊어 돈다.
n1, n2 = 300, 1200
rej = left = 0
for _ in range(B // 200):
    x1 = 4.0 * (rng.lognormal(0.0, 1.0, (200, n1)) - m)
    x2 = rng.lognormal(0.0, 1.0, (200, n2)) - m
    a, b = x1.var(axis=1, ddof=1) / n1, x2.var(axis=1, ddof=1) / n2
    t = (x1.mean(axis=1) - x2.mean(axis=1)) / np.sqrt(a + b)
    nu = (a + b) ** 2 / (a ** 2 / (n1 - 1) + b ** 2 / (n2 - 1))
    c = stats.t(nu).ppf(1 - alpha / 2)
    rej += np.sum(np.abs(t) > c)
    left += np.sum(t < -c)

p = rej / B
print(f"(300, 1200)  Welch 오류율 = {p:.4f}  (왼쪽 {left / B:.4f})")
print(f"             MC 오차 = {np.sqrt(p * (1 - p) / B):.4f}")
print(f"             초과분/왜도 = {(p - 0.05) / 0.3485:.3f}   (예측은 0.053)")

출력:

(300, 1200)  Welch 오류율 = 0.0648  (왼쪽 0.0563)
             MC 오차 = 0.0011
             초과분/왜도 = 0.042   (예측은 0.053)

예측 \(0.0685\)에 모의가 \(0.0648\)로 \(0.004\) 못 미친다. 몬테카를로 오차의 세 배가 넘으므로 우연이 아니다. 초과분/왜도가 \(0.053\)에서 \(0.042\)로 내려앉은 것이고, 비례는 어림이며 초과분이 왜도보다 조금 더 빨리 준다는 뜻이다. 그래도 방향과 자릿수를 잡는 데에는 충분하며, 왼쪽 꼬리에 \(0.0563\) 대 나머지라는 쏠림은 \(n\)이 300이 되어도 그대로 남는다.

셋째, 두 열이 전혀 다른 일을 한다. 합동 \(t\) 열은 \(0.316 \to 0.282\), \(0.339 \to 0.294\)로 표본을 10배 늘려도 거의 움직이지 않는다. Welch 열은 \(0.100 \to 0.058\), \(0.154 \to 0.081\)로 또박또박 내려온다.

까닭은 두 열이 앓는 병이 다르기 때문이다. 합동 \(t\)는 이분산·불균형이라는 병과 비정규성이라는 병을 둘 다 앓는다. 앞의 병은 표본크기가 손댈 수 없으므로(앞 쪽 보기 3) 극한 \(0.277\) 둘레에 눌러앉는다. Welch는 분모를 불편추정량으로 바꿔 앞의 병을 이미 고쳤고, 남은 것은 왜도가 만드는 \(O(n^{-1/2})\) 항뿐이다. 그것은 표본크기가 고쳐 준다. 다만 느리게 고쳐 준다.

(1) 그림이 될 모양은 앞의 두 표에 이미 있다.

• 왼쪽 패널(균형, Welch). 세 색 모두 점선 아래에서 출발해 \(n\)과 함께 올라붙는다. 회색(정규)은 처음부터 점선 위에 있고, 파랑(지수)은 \(0.037 \to 0.046 \to 0.048\), 빨강(로그정규)은 \(0.027 \to 0.040 \to 0.045\)다. 어긋남이 보수적인 쪽이므로 막대가 점선을 뚫고 올라가는 일은 없다.
• 가운데·오른쪽 패널(불균형). 빨간 막대(합동 \(t\))는 \(0.316 \to 0.282\), \(0.339 \to 0.294\)로 거의 평평하고, 파란 막대(Welch)는 \(0.100 \to 0.058\), \(0.154 \to 0.081\)로 점선을 향해 또박또박 내려온다. 두 패널 모두 모든 막대가 점선 위에 있다.

한 그림 안에 세 가지가 함께 있어야 한다. 올라가는 막대(균형·비정규), 내려오는 막대(불균형이지만 Welch), 움직이지 않는 막대(불균형에 합동 \(t\))다.

(2) 그려서 확인한다.

import matplotlib.pyplot as plt
import numpy as np

fig, axes = plt.subplots(1, 3, figsize=(12, 3.8))

# 왼쪽: 균형 설계(n1 = n2 = n). 비정규성만 있는 경우 — 보기 1의 결과.
rng = np.random.default_rng(1)          # 보기 1과 같은 난수열을 쓴다
ns = [10, 30, 100]
xs = np.arange(len(ns))
colors = {"Normal": "#8c8c8c", "Exp": "#4c72b0", "LogN": "#c44e52"}
for k, pop in enumerate(("Normal", "Exp", "LogN")):
    vals = [sum(welch_reject(pop, n)) for n in ns]     # 보기 1의 함수
    off = (k - 1) * 0.27
    axes[0].bar(xs + off, vals, 0.26, color=colors[pop], edgecolor="white", label=pop)
    for x, v in zip(xs + off, vals):
        axes[0].text(x, v + 0.001, f"{v:.3f}", ha="center", fontsize=7)
axes[0].axhline(0.05, color="black", ls="--", lw=1)
axes[0].set_xticks(xs)
axes[0].set_xticklabels([f"n = {n}" for n in ns])
axes[0].set_title(r"Balanced ($n_1=n_2$): Welch $t$")
axes[0].set_ylabel("actual type I error rate")
axes[0].set_ylim(0, 0.080)
axes[0].legend(fontsize=8, loc="upper left", ncol=3)

# 가운데·오른쪽: 치우침 + 이분산 + 불균형 — 보기 3의 결과.
rng = np.random.default_rng(1)          # 보기 3과 같은 난수열을 쓴다
pairs = [(10, 40), (30, 120), (100, 400)]
xp = np.arange(len(pairs))
for ax, pop in zip(axes[1:], ("Exp", "LogN")):
    pool, welch = [], []
    for n1, n2 in pairs:
        rp, rw, *_ = rates(pop, n1, n2)                # 보기 3의 함수
        pool.append(rp)
        welch.append(rw)
    ax.bar(xp - 0.19, pool, 0.36, color="#c44e52", edgecolor="white", label="pooled $t$")
    ax.bar(xp + 0.19, welch, 0.36, color="#4c72b0", edgecolor="white", label="Welch $t$")
    for x, v in zip(xp - 0.19, pool):
        ax.text(x, v + 0.008, f"{v:.3f}", ha="center", fontsize=7)
    for x, v in zip(xp + 0.19, welch):
        ax.text(x, v + 0.008, f"{v:.3f}", ha="center", fontsize=7)
    ax.axhline(0.05, color="black", ls="--", lw=1)
    ax.set_xticks(xp)
    ax.set_xticklabels([f"({a}, {b})" for a, b in pairs], fontsize=8)
    ax.set_xlabel(r"$(n_1,\ n_2)$")
    ax.set_ylim(0, 0.42)
    ax.set_title(rf"{pop}, $\sigma_1/\sigma_2=4$, $n_1{{:}}n_2=1{{:}}4$")
axes[1].set_ylabel("actual type I error rate")
axes[2].legend(fontsize=8, loc="upper right")

plt.tight_layout()
plt.show()

왼쪽 패널에서 세 색깔의 막대가 \(n\)이 커지면서 점선으로 올라붙는다. 오른쪽 두 패널에서 파란 막대(Welch)는 점선을 향해 내려오지만 빨간 막대(합동)는 꼼짝하지 않는다.

같은 그림 안에 고쳐지는 실패와 고쳐지지 않는 실패가 함께 있다. 빨간 막대의 높이를 정하는 것은 분산비와 표본크기비이고, 그 둘은 \(n\)을 키워도 변하지 않는다.

예측한 세 가지가 그대로 나왔다. 막대 위의 수가 보기 1과 보기 3의 표와 한 자리도 다르지 않다. 같은 씨앗으로 같은 계산을 다시 돈 것이니 그래야 맞다.

가리는 것 하나 — 세로축의 눈금이 패널마다 다르다. 왼쪽은 \(0\)에서 \(0.08\)까지이고 가운데·오른쪽은 \(0\)에서 \(0.42\)까지다. 왼쪽 패널에서 점선과 \(0.027\) 막대 사이의 간격이 화면에서 커 보이지만 실제 거리는 \(0.023\)이고, 오른쪽 패널에서 비슷해 보이는 간격은 \(0.1\)이 넘는다. 같은 눈금으로 그리면 왼쪽 패널의 세 막대는 점선에 붙어 구별되지 않는다. 두 실패의 크기가 한 자릿수 다르다는 사실이 축의 재조정에 묻힌다.

가리는 것 둘 — 점선의 위와 아래가 같은 색으로 그려져 있다. 왼쪽 패널의 막대는 점선 아래에 있어 보수적이고, 오른쪽 두 패널의 막대는 점선 위에 있어 관대하다. 전자는 검정력을 잃는 고장이고 후자는 거짓 발견을 쏟는 고장이다. 그림은 "점선에서 얼마나 멀리 있는가"만 보일 뿐 어느 쪽이 위험한 어긋남인지를 말해 주지 않는다. 보기 1에서 로그정규 \(n = 10\)의 \(0.027\)과 보기 3에서 로그정규 \((10,40)\)의 \(0.154\)는 점선에서의 거리가 \(0.023\) 대 \(0.104\)로 네 배 차이지만, 실무에서 감당해야 할 위험으로 보면 종류가 아예 다르다.

그리고 평평해 보이는 빨간 막대가 정말로 평평하지는 않다. \(0.316 \to 0.293 \to 0.282\)는 분명히 내려가고 있다. 다만 \(0.05\)가 아니라 \(0.277\)로 가는 중이라 그림에서 움직임이 보이지 않을 뿐이다(앞 쪽 보기 3). "변하지 않는다"가 아니라 "엉뚱한 값으로 수렴한다"가 맞는 말이다.

적률생성함수 \(E[e^{tZ}] = e^{t^2/2}\)를 쓴다.

이므로

이다. 변동계수는 \(2.1612/1.6487 = 1.311\)이다.

지수분포는 변동계수가 정확히 1이고 왜도가 2다. 로그정규는 변동계수 1.31, 왜도 6.18로 더 치우쳐 있다. 차이는 꼬리에서 더 크다. 초과첨도가 6 대 111로 18배이며, 이것이 보기 1에서 로그정규의 회복이 더 느린 이유다.

꼬리가 두꺼우면 \(S^2\)이 불안정해진다. 5.6절에서 보았듯 \(\text{Var}(S^2)\)이 4차 적률에 비례하므로, 초과첨도 111은 표본분산이 표본마다 크게 널뛴다는 뜻이다. Welch 통계량의 분모가 그 \(S^2\)이다.

두 집단을 같은 모집단에서 같은 크기로 뽑았으므로 \(\bar X_1\)과 \(\bar X_2\)가 같은 분포를 따르고 서로 독립이다. 그러면 \(\bar X_1 - \bar X_2\)와 \(-(\bar X_1 - \bar X_2) = \bar X_2 - \bar X_1\)의 분포가 같다. 즉 차의 분포가 0을 중심으로 대칭이다. 분모도 두 집단을 대칭적으로 쓰므로 \(t\) 통계량 전체가 대칭이고, 두 꼬리의 기각률이 같아야 한다.

교환가능성으로 설명해도 된다. 두 집단의 이름을 바꿔 붙여도 자료의 분포가 변하지 않으므로 통계량의 부호만 뒤집힌다.

일표본이면 다르다. \(\bar X\)의 왜도가 \(\gamma_1/\sqrt n\)으로 남으므로 로그정규 \(n=10\)이면 \(6.18/\sqrt{10} = 1.96\)이다. 같은 모집단에 일표본 \(t\) 검정을 걸어 보면 이렇게 나온다.

명목 2.5%씩이어야 할 두 꼬리가 0.164와 0.001로 갈라지고 합계도 0.165로 부푼다. 보수적인 두 표본 결과(0.027)와 방향까지 반대다. 이것이 두 표본 검정이 일표본 검정보다 비정규성에 강건한 구조적인 이유다.

\(a = 0.9/10 = 0.09\), \(b = 8.0/10 = 0.8\)이므로 \(a + b = 0.89\)이고

이다. \(t_{0.975,\,11.0} = 2.201\)로 합동 자유도 18의 \(2.101\)보다 크다.

같은 모집단에서도 충분히 생긴다. 로그정규의 초과첨도가 111이므로 \(\text{Var}(S^2)\)이 매우 크고, \(n = 10\)에서 두 표본분산이 9배 차이 나는 일은 드물지 않다. 한쪽 표본에 큰 값 하나만 들어와도 \(S^2\)이 몇 배로 뛴다.

여기서 Welch의 자유도가 갖는 성격이 드러난다. \(\nu\)는 모분산이 아니라 표본분산으로 계산되므로 자료의 우연한 불균형에도 반응한다. 치우친 모집단에서는 \(\nu\) 자체가 표본마다 크게 흔들리고, 그 흔들림이 보기 1에서 본 어긋남의 일부를 만든다.

이 때문에 실무에서 \(\nu\)가 이상하게 작게 나오면 자료를 먼저 그려 보아야 한다. 이상점 하나가 만든 값일 수 있다.

차이는 다음과 같다.

\(n\)이 3배가 될 때 차이가 대략 3분의 1에서 2분의 1로 준다. \(O(n^{-1})\)이면 정확히 3분의 1이어야 하므로 대체로 맞고 조금 느리다. 유한표본에서 더 높은 차수의 항이 남아 있기 때문이다.

로그정규에서 차이를 0.001로 줄이려면 \(n = 100\)의 0.005에서 다섯 배를 더 줄여야 하므로 \(n \approx 500\)이 필요하다. 실제로 모의실험을 돌려 보면 \(n = 200\)에서 0.048, \(n = 1000\)에서 0.050이 나온다.

어림규칙에 대한 경고다. "\(n \ge 30\)이면 정규근사"라는 말은 왜도가 1 안팎인 온건한 모집단을 염두에 둔 것이다. 왜도 6, 초과첨도 111짜리 모집단에는 통하지 않는다. 필요한 표본크기는 모집단의 모양에 달려 있으며, 하나의 숫자로 정해지지 않는다.

보기 3이 보여 준 상황이다. 로그정규 \((10,40)\)에서 Welch도 0.154였다. Welch만으로는 부족하다.

1. 설계를 고친다(가능하면 최우선). 두 집단의 크기를 맞추면 왜도가 상쇄되어 오차가 \(O(n^{-1})\)로 떨어진다. 60명 집단에서 12명을 무작위로 뽑아 균형을 맞추는 것도 방법이다. 자료를 버리는 것이 아깝지만, 오류율 0.154를 0.03 근처로 바꾸는 값이 있다. 한계는 자료를 이미 다 모은 뒤라면 검정력을 잃는다는 점이다.

2. 변환한다. 소득처럼 로그정규에 가까운 자료는 로그를 씌우면 거의 정규가 된다. \(t\) 검정의 가정이 회복되고 검정력도 올라간다. 한계는 검정하는 대상이 바뀐다는 것이다. 로그 척도의 평균차는 원래 척도에서 기하평균의 비에 해당하므로, 평균을 비교하고 싶었다면 답이 달라진다.

3. 재표집을 쓴다. 부트스트랩 \(t\)(스튜던트화 부트스트랩)나 순열검정은 정규성을 가정하지 않는다. 특히 부트스트랩 \(t\)는 \(O(n^{-1})\)의 정확도를 갖는 것으로 알려져 있어 치우친 자료에서 \(t\) 근사보다 낫다. 한계는 계산이 필요하고, 순열검정의 귀무가설이 "두 분포가 같다"라서 분산이 다르면 평균만 비교하는 것이 아니게 된다는 점이다.

어느 쪽도 하지 말아야 할 것. 합동 \(t\)를 쓰는 것이다. 보기 3에서 \(n\)을 10배로 늘려도 0.294였다.

유도. 누적률은 독립인 확률변수의 합에서 더해지고, 상수배에서 \(\kappa_3(cX) = c^3\kappa_3(X)\)로 변한다. 한 집단에서

이다. 모집단의 왜도가 \(\gamma_1 = \kappa_3(X)/\sigma^3\)이므로 \(\kappa_3(X) = \gamma_1\sigma_1^3\)이고

이다. 둘째 집단은 부호가 바뀌어 들어간다. \(\kappa_3(-Y) = -\kappa_3(Y)\)이고 두 표본이 독립이므로

이다. \(\square\)

0이 될 조건은 \(\sigma_1^3/n_1^2 = \sigma_2^3/n_2^2\), 즉

이다. 같은 모집단에서 같은 크기로 뽑으면(\(\sigma_1=\sigma_2\), \(n_1=n_2\)) 자동으로 만족된다. 흥미롭게도 분산이 다를 때도 표본크기를 맞추면 왜도를 지울 수 있다. 다만 그 비 \((\sigma_1/\sigma_2)^{3/2}\)는 차의 분산을 최소화하는 최적 배분 \(\sigma_1/\sigma_2\)와 다르므로, 왜도를 지우는 배분과 정밀도를 높이는 배분이 일치하지는 않는다.

수치. \(\sigma_1 = 4\), \(\sigma_2 = 1\), \((n_1,n_2)=(10,40)\)이면

이므로 표준화한 왜도는

이다. 지수모집단은 \(\gamma_1 = 2\)이므로 \(0.617\), 로그정규는 \(\gamma_1 = 6.18\)이므로 \(1.909\)다.

보기 3과 맞춰 보자. 왜도가 큰 쪽(로그정규 1.909)의 오류율이 0.154, 작은 쪽(지수 0.617)이 0.100으로 왜도가 큰 쪽이 더 많이 어긋난다. \((100,400)\)에서는 왜도가 각각 0.604와 0.195로 줄고 오류율도 0.081과 0.058로 함께 내려간다. 왜도가 \(1/\sqrt n\)으로 줄고 오류율의 어긋남도 그 속도를 따른다.

\(64/100\)과 \(1/1600\)을 견주면 첫 항이 1000배 크다. 왜도의 거의 전부가 작은 표본·큰 분산 쪽 집단에서 온다. 그 집단이 차의 분산도 지배하고(앞 쪽에서 계산한 \(\nu = 9.3\)) 왜도도 지배한다.

로그정규의 두 평균. \(X = e^{\mu+\sigma Z}\)이면

이다. \(\mu\)가 같고 \(\sigma\)만 다르면 기하평균은 같고 산술평균은 다르다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
n, REP = 40, 20_000
mu1, s1 = 0.0, 0.5
mu2, s2 = 0.0, 1.2
print(f"  A: 기하평균 {np.exp(mu1):.4f}, 산술평균 {np.exp(mu1+s1**2/2):.4f}")
print(f"  B: 기하평균 {np.exp(mu2):.4f}, 산술평균 {np.exp(mu2+s2**2/2):.4f}\n")

A = rng.lognormal(mu1, s1, (REP, n))
B = rng.lognormal(mu2, s2, (REP, n))
p_raw = stats.ttest_ind(A, B, axis=1, equal_var=False).pvalue
p_log = stats.ttest_ind(np.log(A), np.log(B), axis=1, equal_var=False).pvalue
print(f"  원자료 t  기각률 {(p_raw<0.05).mean():.4f}")
print(f"  로그 후 t 기각률 {(p_log<0.05).mean():.4f}")

출력:

  A: 기하평균 1.0000, 산술평균 1.1331
  B: 기하평균 1.0000, 산술평균 2.0544

  원자료 t  기각률 0.3821
  로그 후 t 기각률 0.0515

둘 다 옳게 동작하고 있다. 산술평균은 \(1.13\) 대 \(2.05\)로 실제로 다르므로 원자료 \(t\)가 \(38\%\) 기각하는 것이 맞고, 기하평균은 \(1.00\)으로 같으므로 로그 후 \(t\)가 \(5\%\)만 기각하는 것도 맞다.

그러나 분석자가 무엇을 물었는지에 따라 하나는 틀린 답이다.

"정규성을 만족시키려고" 로그를 씌우는 것은 위험하다. 로그 변환은 분포를 예쁘게 만드는 기술적 조작이 아니라 비교 대상을 산술평균에서 기하평균으로 바꾸는 일이다. 로그정규에서 로그 후 \(t\)는 사실상 중앙값 비교이기도 하다(로그정규의 중앙값이 \(e^\mu\)).

실질적 함의가 큰 분야가 있다.

• 소득·자산. 총액이나 세수가 관심이면 산술평균이 맞다. 로그를 씌우면 상위 소득의 영향이 지워져 "전형적인 가구"의 비교가 된다.
• 약물 농도(생물학적 동등성). 관례적으로 로그 척도에서 비교하며, 그래서 결론이 "기하평균비가 \(0.8\sim1.25\) 안에 있다"로 표현된다.
• 비용 자료. 총비용을 예측해야 하면 산술평균이다. 로그 후 분석하고 지수를 취해 되돌리면 산술평균을 과소추정한다.

되돌릴 때의 함정. $\exp(\overline{\log X}) $는 기하평균이지 산술평균이 아니다. 산술평균을 원하면 \(\exp(\bar\ell + s^2/2)\) 같은 보정이 필요하며(5.2절 델타 방법), 이 보정을 빼먹는 것이 실무에서 흔한 오류다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
n1, n2, REP, B = 15, 15, 1000, 499
rej = {'welch': 0, 'perm': 0, 'boot': 0}

for _ in range(REP):
    a = rng.lognormal(0, 1, n1)
    b = rng.lognormal(0, 1, n2)          # 같은 분포 -> 귀무가설 참
    rej['welch'] += stats.ttest_ind(a, b, equal_var=False).pvalue < 0.05

    obs = a.mean() - b.mean()
    pool = np.concatenate([a, b])
    cnt = 0
    for _ in range(B):
        p = rng.permutation(pool)
        cnt += abs(p[:n1].mean() - p[n1:].mean()) >= abs(obs)
    rej['perm'] += (cnt + 1) / (B + 1) < 0.05

    # 부트스트랩 t (각 집단에서 따로 재표본)
    ta = a - a.mean(); tb = b - b.mean()
    ts = []
    for _ in range(B):
        aa = ta[rng.integers(0, n1, n1)]
        bb = tb[rng.integers(0, n2, n2)]
        se = np.sqrt(aa.var(ddof=1)/n1 + bb.var(ddof=1)/n2)
        ts.append(abs(aa.mean() - bb.mean()) / se)
    se0 = np.sqrt(a.var(ddof=1)/n1 + b.var(ddof=1)/n2)
    rej['boot'] += abs(obs)/se0 > np.percentile(ts, 95)

for k, v in rej.items():
    print(f"  {k:>6}: {v/REP:.4f}")

출력:

   welch: 0.0350
    perm: 0.0520
    boot: 0.0220

예상과 다른 결과가 나왔다. 셋 다 명목을 넘을 것 같았는데, 오히려 웰치와 부트스트랩이 보수적이다. 순열만 \(0.052\)로 명목에 가깝다. (모의 표준오차가 \(\sqrt{0.05\times0.95/1000} \approx 0.007\)이므로 \(0.035\)와 \(0.022\)는 우연으로 보기 어렵다.)

웰치가 보수적인 것은 이 쪽 연습문제 2의 결과 그대로다. 두 집단이 같은 분포이고 표본크기도 같으면 양쪽 꼬리의 왜곡이 서로 상쇄된다. 한쪽 꼬리는 부풀고 다른 쪽은 줄어드는데 양측검정에서 둘을 더하므로 순효과가 작고, 여기서는 오히려 아래로 간다. 불균형하거나 두 모집단이 다르면 상쇄가 깨지고 오류율이 올라간다(연습문제 10).

부트스트랩 \(t\)가 가장 보수적이다(\(0.022\)). \(n=15\)에서 재표본이 원표본의 꼬리를 넘어서지 못하므로, 귀무분포의 꼬리가 실제보다 두껍게 추정된다. 임계값이 커지고 기각이 줄어든다. 5.6절 지수모집단 연습문제 10에서 본 "부트스트랩이 꼬리 의존 문제에서 덜 듣는다"는 현상의 다른 얼굴이다.

순열이 가장 정확한 이유는 가정이 실제로 성립하기 때문이다. 두 집단이 같은 분포에서 나왔으므로 표지를 섞는 것이 정확히 타당하고, 그러면 순열검정은 유한표본에서 정확한 검정이 된다. 분포 모양과 무관하게 성립하는 성질이다.

다만 이 이점은 조건부다. 5.9절에서 본 대로 두 분포의 모양이 다르면 순열의 전제가 깨진다. 여기서는 "분산도 같고 모양도 같고 평균만 비교"하는 이상적인 상황이라 순열이 빛난 것이다.

보수적인 것이 안전하다는 뜻은 아니다. 오류율이 \(0.022\)라는 것은 검정력을 그만큼 잃었다는 뜻이기도 하다. 명목 수준을 지키는 것과 그 아래로 크게 내려가는 것은 다른 문제이며, 후자는 실제 차이를 놓치게 만든다.

만–휘트니가 검정하는 것은 평균이 아니다. 통계량은 "\(A\)에서 뽑은 값이 \(B\)에서 뽑은 값보다 클" 횟수를 세므로, 귀무가설은

이다. 이것을 확률적 우위가 없다는 진술로 읽는다(2장 ECDF 연습문제 10).

"중앙값 검정"이라는 흔한 설명은 부정확하다. 두 분포의 모양이 같고 위치만 다르다는 추가 가정(위치이동 모형)이 있을 때만 중앙값 비교가 된다. 모양이 다르면 중앙값이 같아도 \(P(X_A > X_B) \ne 1/2\)일 수 있다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(4)
# 중앙값은 같고 산포만 다른 두 분포
a = rng.lognormal(0, 0.3, 200_000)
b = rng.lognormal(0, 1.5, 200_000)
print(f"  중앙값 A={np.median(a):.4f}  B={np.median(b):.4f}")
print(f"  평균   A={a.mean():.4f}  B={b.mean():.4f}")
print(f"  P(A > B) = {(a > b).mean():.4f}")

출력:

  중앙값 A=1.0016  B=0.9930
  평균   A=1.0473  B=3.0668
  P(A > B) = 0.5028

중앙값이 같으니 \(P(A>B) \approx 0.50\)이고 만–휘트니는 차이를 찾지 못한다. 그런데 산술평균은 \(1.05\) 대 \(3.07\)로 세 배 차이 난다. 만–휘트니로 "차이 없음"을 얻고 "평균이 같다"고 읽으면 완전히 틀린다.

정리하면.

"비모수가 안전하다"는 말이 오해를 부른다. 만–휘트니는 정규성 가정을 없애는 대신 검정 대상을 바꾼다. 총액·합계가 중요한 문제(소득세, 총비용, 총수요)에서는 산술평균이 관심사이므로 만–휘트니가 답이 되지 않는다.

선택 기준은 통계적 성질이 아니라 질문이다. 무엇을 비교하고 싶은지 먼저 정하고, 그다음에 그 양에 맞는 방법 중 가정이 덜 까다로운 것을 고른다. 16장에서 비모수 방법을 다룰 때 이 원칙을 반복한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
REP = 40_000
print(f"{'n1,n2':>10}{'치우친 쪽':>12}{'실제 오류율':>14}")
for n1, n2 in [(15,15), (10,40), (40,10), (5,45), (45,5)]:
    a = rng.lognormal(0, 1, (REP, n1))       # 치우친 쪽은 늘 1번
    b = rng.normal(np.exp(0.5), np.sqrt((np.e-1)*np.e), (REP, n2))
    p = stats.ttest_ind(a, b, axis=1, equal_var=False).pvalue
    who = "작은 표본" if n1 < n2 else ("큰 표본" if n1 > n2 else "동일")
    print(f"{f'{n1},{n2}':>10}{who:>12}{(p<0.05).mean():>14.4f}")

출력:

     n1,n2       치우친 쪽        실제 오류율
     15,15          동일        0.0530
     10,40       작은 표본        0.0933
     40,10        큰 표본        0.0507
      5,45       작은 표본        0.1409
      45,5        큰 표본        0.0564

치우친 모집단이 작은 표본 쪽에 있으면 훨씬 나쁘다. \((5,45)\)에서 오류율이 \(0.141\)로 명목의 거의 세 배인데, 표본크기만 뒤집은 \((45,5)\)에서는 \(0.056\)으로 거의 정상이다.

이유는 중심극한정리가 어디에 적용되는가다. 차 \(\bar X_1 - \bar X_2\)의 왜도는 연습문제 6의

에서 오는데, 치우친 집단의 \(n\)이 분모에 제곱으로 들어간다. 그 집단의 표본이 작으면 왜도가 거의 그대로 남는다. 반대로 치우친 집단에 표본을 많이 주면 \(\bar X_1\)이 충분히 정규화되어 문제가 사라진다.

실무 지침 셋.

• 치우친 집단에 표본을 더 배정하라. 5.7절 연습문제 10의 비용·분산 기반 배분에 "왜도"라는 기준이 하나 더 붙는 셈이다.
• 왜도를 먼저 진단하라. 두 집단의 왜도와 표본크기를 함께 보고, 위험한 조합인지 판단한다.
• 위험하면 연습문제 5의 대처로 옮겨라. 균형을 맞추거나 척도를 바꾸는 쪽이 먼저다. 재표집은 작은 쪽 표본이 다섯 개밖에 안 되면 크게 도움이 되지 않는다. 연습문제 8에서 본 대로 부트스트랩 \(t\)는 원표본의 꼬리를 넘어서지 못한다.

한 줄로. 비정규성과 불균형은 따로 있을 때보다 함께 있을 때 훨씬 위험하며, 특히 치우친 집단의 표본이 작을 때 그렇다. \(\square\)