카이제곱 분포 (chi-squared)¶

(1) 이론값. 모집단이 정규이므로 \(V = (n-1)S^2/\sigma^2\)은 근사가 아니라 정확히 \(\chi^2_{n-1} = \chi^2_9\)을 따른다. 그러면 더 할 일이 없다. 4.2절이 구해 둔 \(\chi^2_d\)의 적률에 \(d = 9\)를 넣으면 된다.

\(S^2\) 쪽으로 옮겨 적으면 \(E[S^2] = \sigma^2 = 4\)이고

이다. 왜도가 \(0.94\)나 된다는 것도 새겨 둘 만하다. \(n = 10\)에서 \(S^2\)의 표본분포는 대칭이 아니며, 그것이 \(\sigma^2\)의 신뢰구간이 \(s^2\)을 중심으로 대칭이 아닌 까닭이다(연습문제 3).

두 모수가 왜 사라지는가. \(\mu\)는 \(S^2\)이 편차 \(X_i - \bar X\)만 쓰기 때문에 지워진다. 모든 관측값에 같은 수를 더하면 \(X_i\)와 \(\bar X\)가 나란히 밀려 편차가 그대로 남는다. \(\sigma\)는 식이 이미 \(\sigma^2\)으로 나누어져 있어 약분된다. 그러므로 \(V\)의 분포를 정하는 것은 \(n\) 하나뿐이다. 코드가 \(\mu = 1\), \(\sigma = 2\)라는 특정한 값을 쓴 것은 결과가 그 값에 매이지 않음을 드러내려는 선택이고, 아무 값을 넣어도 같은 그림이 나온다.

(2) 모의실험.

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

n, n_sim, mu, sigma = 10, 10_000, 1, 2
# 크기 10짜리 정규 표본을 1만 개 만든다. 행 하나가 표본 하나다.
samples = stats.norm(loc=mu, scale=sigma).rvs(size=(n_sim, n))
s = samples.std(axis=1, ddof=1)      # 행마다 표본표준편차(n-1로 나눔)

# 이 통계량이 정확히 chi^2(n-1) 을 따른다는 것이 정리의 내용이다.
# mu = 1 을 썼지만 결과는 mu에 의존하지 않는다.
# S^2 이 편차만 쓰므로 위치가 상쇄되기 때문이다.
data = (n - 1) * s**2 / sigma**2

fig, ax = plt.subplots(figsize=(12, 3))
_, bins, _ = ax.hist(data, bins=100, density=True, alpha=0.7)
ax.plot(bins, stats.chi2(n-1).pdf(bins), '--r', lw=3, label='χ²(n-1) PDF')
ax.set_title('(n-1)S²/σ² from Normal Population → χ² Exact')
ax.legend()
ax.spines[['top', 'right']].set_visible(False)
plt.show()

히스토그램과 \(\chi^2_9\) 곡선이 빈틈없이 겹친다. 근사가 아니라 정확한 결과이니 당연한 일이다.

그러나 그림만으로는 "겹친다"가 어느 정도인지 말할 수 없으니 적률을 직접 재어 본다. 위 블록은 난수 씨앗을 고정하지 않아 돌릴 때마다 히스토그램이 조금씩 달라지므로, 수치를 맞추려면 씨앗을 박아야 한다. 적률의 몬테카를로 오차를 줄이려 되풀이도 20만 번으로 늘린다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, B, mu, sigma = 10, 200_000, 1, 2
d = n - 1

# 위 그림과 같은 설정이지만 씨앗을 고정하고 되풀이를 20만 번으로 늘려 수치를 잰다.
s = rng.normal(mu, sigma, size=(B, n)).std(axis=1, ddof=1)
V = (n - 1) * s**2 / sigma**2

print(f"E[V]    이론 {d}        모의 {V.mean():.4f}   (MC 오차 {V.std(ddof=1)/np.sqrt(B):.4f})")
print(f"sd(V)   이론 {np.sqrt(2*d):.4f}   모의 {V.std(ddof=1):.4f}")
print(f"왜도    이론 {np.sqrt(8/d):.4f}   모의 {stats.skew(V):.4f}")
print(f"P(V > {stats.chi2(d).ppf(0.95):.3f})  이론 0.0500   모의 {np.mean(V > stats.chi2(d).ppf(0.95)):.4f}")
print(f"KS 검정 p = {stats.kstest(V, 'chi2', args=(d,)).pvalue:.4f}")

# mu 를 1 에서 1000 으로 바꿔도 분포가 그대로인지 확인한다.
V2 = (n - 1) * rng.normal(1000, sigma, size=(B, n)).std(axis=1, ddof=1) ** 2 / sigma**2
print(f"mu = 1000 으로 바꾸면   E[V] {V2.mean():.4f}   sd(V) {V2.std(ddof=1):.4f}")

출력:

E[V]    이론 9        모의 8.9982   (MC 오차 0.0095)
sd(V)   이론 4.2426   모의 4.2511
왜도    이론 0.9428   모의 0.9670
P(V > 16.919)  이론 0.0500   모의 0.0502
KS 검정 p = 0.7181
mu = 1000 으로 바꾸면   E[V] 8.9949   sd(V) 4.2462

세 적률이 모두 맞는다. 평균 \(8.9982\)는 이론값 \(9\)에서 몬테카를로 오차 \(0.0095\)의 \(0.2\)배만큼 떨어져 있고, 표준편차 \(4.2511\)도 이론값 \(4.2426\)에 붙는다. 꼬리확률 \(0.0502\)는 명목 \(0.05\)와 사실상 같다. 왜도 \(0.9670\)만 이론값 \(0.9428\)보다 조금 큰데, 왜도는 3차 적률이라 몬테카를로 오차가 \(0.011\)쯤이어서 그 두 배 안이다. 씨앗을 열두 개 바꿔 평균하면 \(0.9411\)로 내려앉으니 체계적 차이가 아니다. 분포 전체를 한꺼번에 묻는 KS 검정의 \(p = 0.72\)가 같은 말을 한다. \(\chi^2_9\)과 구별되지 않는다.

맨 아랫줄이 (1)의 둘째 물음에 답한다. \(\mu\)를 \(1\)에서 \(1000\)으로 옮겨도 \(E[V] = 8.9949\), \(\operatorname{sd}(V) = 4.2462\)로 조금도 달라지지 않는다. \(\sigma\) 쪽은 식에서 이미 나누어 없앴으므로 따로 확인할 것도 없다.

(1) 이론값. 둘을 갈라 보아야 한다. 먼저 중심은 멀쩡하다.

5.1절에서 증명한 \(E[S^2] = \sigma^2\)은 i.i.d.와 유한한 분산만 요구하고 정규성은 요구하지 않는다. 그러므로 지수모집단에서도 \(E[V]\)는 \(\chi^2_9\)의 평균과 정확히 같다.

어긋나는 것은 폭이다. 본문의

에 \(n = 10\), \(\beta_2 = 9\), \(\sigma^2 = 1\)을 넣으면 \(\operatorname{Var}(S^2) = \frac{1}{10}\!\left(9 - \frac79\right) = \frac{74}{90} = 0.8222\)이고

이다. \(\chi^2_9\)이 믿는 값은 \(2d = 18\)이므로 배율이

다. 표준편차로 옮기면 \(\sqrt{66.6} = 8.1609\) 대 \(\sqrt{18} = 4.2426\)으로 \(1.92\)배다. 연습문제 5가 구하는 극한 배율 \((\beta_2-1)/2 = 4\)에는 아직 못 미치는데, 그 차이는 유한한 \(n\)에서 오는 것이고 \(n\)을 키우면 \(4\)로 올라간다. 좋아지는 쪽이 아니라 나빠지는 쪽으로 올라간다는 것이 이 절의 주제다.

한 줄로 줄이면 중심은 정확히 맞고 폭만 틀리다. 불편성은 정규성을 묻지 않지만 표본분포는 묻기 때문이다.

(2) 모의실험.

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

n, n_sim = 10, 10_000
# 앞 모의실험과 **한 곳만** 다르다. 모집단을 정규에서 지수로 바꿨다.
# 나머지 코드는 그대로다.
samples = stats.expon().rvs(size=(n_sim, n))
s = samples.std(axis=1, ddof=1)
# Exp(1)의 분산이 1이므로 sigma^2으로 나눌 필요가 없다.
data = (n - 1) * s**2

# 지수분포는 오른쪽으로 크게 치우쳐 4차 적률이 크다.
# 그 결과 S^2 의 분포가 카이제곱보다 훨씬 무거운 꼬리를 갖게 되어
# 아래 그림에서 히스토그램이 빨간 곡선 밖으로 크게 벗어난다.

fig, ax = plt.subplots(figsize=(12, 3))
_, bins, _ = ax.hist(data, bins=100, density=True, alpha=0.7)
ax.plot(bins, stats.chi2(n-1).pdf(bins), '--r', lw=3, label='χ²(n-1) PDF')
ax.set_title('(n-1)S²/σ² from Exponential Population → χ² Approximation Fails')
ax.legend()
ax.spines[['top', 'right']].set_visible(False)
plt.show()

히스토그램이 빨간 곡선을 크게 벗어난다. 왼쪽으로는 더 몰려 있고 오른쪽으로는 곡선 밖까지 길게 뻗는다. 앞 절에서 표본평균이 지수 모집단에서도 얌전히 종 모양으로 수렴하던 것과 견주어 보라. 모집단을 바꾼 것은 똑같은데 한쪽은 아무 일도 없었고 다른 쪽은 무너졌다. 평균 계열과 분산 계열의 차이가 이것이다.

(1)의 두 수치를 그림에서 읽어 낼 수는 없으니 다시 씨앗을 박고 재어 본다. 이 블록도 앞 보기와 같은 이유로 되풀이를 20만 번으로 늘렸다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, B = 10, 200_000
d, beta2 = n - 1, 9.0                  # Exp(1) 의 첨도는 9 다

V = d * rng.exponential(size=(B, n)).var(axis=1, ddof=1)

# 일반 모집단의 Var(S^2) = (1/n)(beta2 - (n-3)/(n-1)) sigma^4 에서 온 값이다.
var_th = d**2 * (beta2 - (n - 3) / (n - 1)) / n
print(f"E[V]    카이제곱 {d}        이론 {d}        모의 {V.mean():.4f}")
print(f"sd(V)   카이제곱 {np.sqrt(2*d):.4f}   이론 {np.sqrt(var_th):.4f}   모의 {V.std(ddof=1):.4f}")
print(f"배율    Var(V) / {2*d} = {var_th/(2*d):.2f}   (극한 (beta2-1)/2 = {(beta2-1)/2:.1f})")
print(f"중앙값  카이제곱 {stats.chi2(d).median():.4f}             모의 {np.median(V):.4f}")

lo, hi = stats.chi2(d).ppf([0.025, 0.975])
print(f"명목 95% 구간 [{lo:.3f}, {hi:.3f}] 밖 = {np.mean((V < lo) | (V > hi)):.4f}"
      f"   (아래로 {np.mean(V < lo):.4f}, 위로 {np.mean(V > hi):.4f})")

출력:

E[V]    카이제곱 9        이론 9        모의 9.0027
sd(V)   카이제곱 4.2426   이론 8.1609   모의 8.1568
배율    Var(V) / 18 = 3.70   (극한 (beta2-1)/2 = 4.0)
중앙값  카이제곱 8.3428             모의 6.6685
명목 95% 구간 [2.700, 19.023] 밖 = 0.2322   (아래로 0.1413, 위로 0.0909)

예측한 대로 갈린다. 평균 \(9.0027\)은 \(\chi^2_9\)의 \(9\)와 맞고, 표준편차 \(8.1568\)은 \(\chi^2_9\)의 \(4.2426\)이 아니라 (1)이 계산한 \(8.1609\)와 맞는다. 중심은 같은 자리에 둔 채 폭만 \(1.92\)배 넓은 분포인 것이다.

평균이 같은데 중앙값이 \(6.67\) 대 \(8.34\)로 크게 아래에 있다는 것이 그림의 "왼쪽으로 몰려 있다"의 정체다. 평균을 제자리에 붙들어 두면서 몸통을 왼쪽으로 옮기려면 오른쪽 꼬리가 아주 멀리까지 뻗어야 하고, 그래서 두 어긋남이 한 그림에 함께 나타난다.

구간으로 옮기면 명목 5%가 실제 23.2%다. 아래로 \(14.1\%\), 위로 \(9.1\%\)로 양쪽이 고르게 어긋나지도 않는다. 연습문제 6이 같은 지수모집단·같은 \(n = 10\)에서 \(0.233\)을 얻은 것이 이 수치이고, 그 문제는 \(n\)을 키우면 이 값이 연습문제 5의 극한 \(1 - 0.673 = 0.327\)을 향해 올라간다는 것까지 보인다. 여기서 재어 둘 것은 그 출발점 하나다.

모집단이 정규이므로 표본분포는 정확하다:

구해야 할 것은 \(P(S^2 > 15) = P\!\left(\frac{19 S^2}{9} > \frac{19 \times 15}{9}\right) = P(\chi^2_{19} > 31.67)\)이다.

카이제곱 분포표나 소프트웨어에서 \(P(\chi^2_{19} > 31.67) \approx 0.034\)이다. \(\sigma^2 = 9\)일 때 표본분산이 15를 넘을 확률은 약 3.4%이다.

양변을 \(\sigma^2\)으로 나눈다. 왼쪽은 표준정규 \(n\)개의 제곱합이므로 \(\chi^2_n\)이다. 오른쪽 둘째 항은

로 표준정규 하나의 제곱, 곧 \(\chi^2_1\)이다.

정규모집단에서 \(\bar X\)와 \(S^2\)이 독립이므로 오른쪽 두 항도 독립이고, 독립인 카이제곱을 더하면 자유도가 더해진다. 왼쪽의 자유도가 \(n\)이고 둘째 항의 자유도가 1이므로 첫째 항의 자유도는 \(n-1\)일 수밖에 없다. 곧 \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)이다. \(\square\)

자유도 하나는 사라진 것이 아니라 옮겨 간 것이다. \(\mu\)를 모르는 탓에 제곱합의 한 조각이 "\(\bar X\)가 \(\mu\)에서 얼마나 떨어져 있는가"를 재는 데 쓰였고, 남은 \(n-1\)조각만 퍼짐을 재는 데 남는다. 자유도를 "제약 하나당 하나씩 준다"고 외우는 것보다 이 분해를 기억해 두는 편이 낫다.

유도가 어디에 기대고 있는지도 함께 보아 두자. 두 항의 독립성이 없으면 자유도를 더할 수 없고, 그 독립성이 바로 정규분포만의 성질이다(연습문제 7). 자유도 \(n-1\)이라는 익숙한 숫자조차 정규성 위에 서 있다.

from scipy import stats

for d in (5, 19, 50, 200, 1000):
    lo, hi = stats.chi2(d).ppf([0.025, 0.975])
    print(f"d={d:>5}   아래 {d-lo:8.3f}   위 {hi-d:8.3f}   비 {(hi-d)/(d-lo):.3f}")

출력:

d=    5   아래    4.169   위    7.833   비 1.879
d=   19   아래   10.093   위   13.852   비 1.372
d=   50   아래   17.643   위   21.420   비 1.214
d=  200   아래   37.272   위   41.058   비 1.102
d= 1000   아래   85.743   위   89.531   비 1.044

자유도 5에서 위쪽 거리가 아래쪽의 1.88배였던 것이 자유도 1000에서 1.04배까지 내려온다. 4.2절이 구한 왜도 \(\sqrt{8/d}\)가 0으로 줄어드는 것과 같은 이야기이며, 자유도가 커질수록 카이제곱이 정규에 가까워지므로 두 꼬리도 균형을 찾는다. 분산의 구간이 \(s^2\)을 중심으로 대칭이 아닌 까닭, 그리고 표본이 커지면 그 비대칭이 옅어지는 까닭이 이 표에 있다.

한 가지를 구별해 두어야 한다. 여기서 대칭에 가까워지는 것은 \(\chi^2_d\) 자체의 모양이다. \((n-1)S^2/\sigma^2\)이 그 \(\chi^2_d\)를 따르느냐는 전혀 다른 물음이고, 자유도를 키운다고 답이 달라지지 않는다. 본문의 "모양이 닮아 가는 것과 폭이 맞는 것은 다른 이야기"라는 문장이 이 구별을 가리킨다.

타당하지 않다. \(\frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}\)이라는 결과는 모집단이 정규일 때만 성립한다. 지수분포는 오른쪽으로 치우쳐 있고 첨도가 \(\beta_2 = 9\)여서 \(S^2\)의 분포가 \(\chi^2_9\) 분포보다 훨씬 두꺼운 꼬리를 갖게 된다.

구체적으로 정규가 아닌 모집단에서 \(S^2\)의 분산은 첨도에 의존한다: \(\text{Var}(S^2) \approx \frac{2\sigma^4}{n-1}\cdot\frac{\beta_2-1}{2}\). 지수분포에서는 \(\frac{2\sigma^4}{9}\cdot 4 = \frac{8\sigma^4}{9}\)가 되어 카이제곱 이론이 예측하는 값(\(\frac{2\sigma^4}{9}\))보다 4배 크다. 카이제곱 가정에 기반한 신뢰구간과 가설검정은 포함확률이 틀리고 제1종 오류율이 부풀려진다.

표본을 키우면 낫는가. 낫지 않는다. 배율 4는 \(n\)에 의존하지 않으므로 \(n = 10\)이든 \(n = 1000\)이든 그대로 남는다. 다음 문제가 그 이유를 식으로 확인한다.

배율. 두 값의 비는

이고 \(n \to \infty\)에서 \(\frac{n-1}{2n} \to \frac12\), \(\frac{n-3}{n-1} \to 1\)이므로 극한은 \(\frac{\beta_2-1}{2}\)이다. \(\square\)

정규분포는 \(\beta_2 = 3\)이라 배율이 1인데, 이 경우에는 극한이 아니라 모든 \(n\)에서 정확히 1이다. 카이제곱 결과가 근사가 아니라 등식인 것이 여기에 다시 나타난다. 균등분포는 \((1.8-1)/2 = 0.4\), 지수분포는 \((9-1)/2 = 4\)다. 표준편차로 옮기면 각각 \(\sqrt{0.4} = 0.63\)배와 \(\sqrt 4 = 2\)배다.

극한 포함률. \(n\)이 크면 \(S^2\)의 분포가 정규로 다가가고 그 표준편차는 \(\sigma^2\sqrt{(\beta_2-1)/n}\)이다. 그런데 카이제곱 구간은 표준편차를 \(\sigma^2\sqrt{2/n}\)로 믿고 그 \(\pm 1.96\)배만큼 뻗는다. 참 표준편차를 단위로 재면 구간이

만큼만 뻗는 셈이므로 극한 포함률은 \(2\Phi\!\left(1.96\sqrt{\frac{2}{\beta_2-1}}\right) - 1\)이다. 넣어 보면

• 균등: \(1.96\sqrt{2/0.8} = 3.10\), 포함률 \(2\Phi(3.10)-1 = 0.998\)
• 정규: \(1.96\sqrt{2/2} = 1.96\), 포함률 \(0.950\) (검산이 맞는다)
• 지수: \(1.96\sqrt{2/8} = 0.98\), 포함률 \(2\Phi(0.98)-1 = 0.673\)

이다. 지수 모집단에서 명목 95% 구간이 실제로는 67%만 담는다.

배율에 \(n\)이 없다는 것이 요점이다. 균등 쪽은 지나치게 넓어 보수적이고 지수 쪽은 지나치게 좁아 위험하지만, 어느 쪽이든 표본을 키운다고 0.95로 돌아오지 않는다. 첨도 하나만 알면 카이제곱 구간이 어느 방향으로 얼마나 망가지는지 미리 계산할 수 있다는 점도 함께 기억해 둘 만하다. 다음 문제가 유한한 \(n\)에서 실제 값을 재어 극한과 맞춰 보고, 5.6절이 같은 계산을 모집단별로 확인한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
B = 200_000

print(f"{'n':>6}{'(n-1)S^2 오류율':>18}{'Xbar |Z|>1.96':>16}")
for n in (10, 50, 200, 1000):
    d = n - 1
    lo, hi = stats.chi2(d).ppf([0.025, 0.975])
    X = rng.exponential(size=(B, n))        # Exp(1): mu = sigma^2 = 1
    V = d * X.var(axis=1, ddof=1)           # (n-1)S^2/sigma^2
    Z = (X.mean(axis=1) - 1) * np.sqrt(n)   # (Xbar - mu)/(sigma/sqrt n)
    print(f"{n:>6}{np.mean((V < lo) | (V > hi)):>18.3f}"
          f"{np.mean(np.abs(Z) > 1.96):>16.3f}")

출력:

     n      (n-1)S^2 오류율   Xbar |Z|>1.96
    10             0.233           0.045
    50             0.297           0.049
   200             0.317           0.050
  1000             0.324           0.050

두 열이 정반대로 움직인다. 오른쪽 열은 \(n = 10\)에서 이미 0.045이고 표본을 키우면 0.050에 더 바짝 붙는다. 중심극한정리가 약속한 대로다. 왼쪽 열은 \(n = 10\)에서 0.233이고 표본을 키우면 더 나빠져 0.324로 간다. 포함률로 옮기면 \(1 - 0.324 = 0.676\)이고, 연습문제 5가 계산한 극한 포함률 0.673에 수렴하는 중이다.

작은 표본에서 오히려 덜 나빠 보이는 것도 눈여겨볼 만하다. \(n\)이 작을 때는 \(\chi^2_{n-1}\) 자체가 넓어 어긋남을 얼마간 가려 주는데, 표본이 커지면 그 완충이 사라지고 첨도의 효과만 남는다. 자료를 더 모으면 결론이 더 믿을 만해진다는 통상의 기대가 여기서는 성립하지 않는다.

같은 모집단에서 같은 표본으로 잰 두 수치라는 점이 이 표의 힘이다. 모집단이 문제가 아니라 통계량이 몇 차 적률에 매여 있는가가 문제라는 것을 한 줄에 보여 준다.

공분산부터. \(\mu_3 = E[(X-\mu)^3]\)이 3차 중심적률일 때 \(\text{Cov}(\bar X, S^2) = \mu_3/n\)이다. 대칭인 모집단은 \(\mu_3 = 0\)이므로 두 통계량이 무상관이고, 정규분포에서는 거기에 더해 무상관이 곧 독립까지 간다. 치우친 모집단에서는 \(\mu_3 \ne 0\)이라 무상관조차 깨진다.

상관계수. \(\text{Var}(\bar X) = \sigma^2/n\)이고 \(\text{Var}(S^2) \approx (\beta_2-1)\sigma^4/n\)이므로

이다. 분자와 분모에 모두 \(1/n\)이 있어 \(n\)이 약분된다. 곧 이 상관은 모집단의 왜도를 \(\sqrt{\beta_2-1}\)로 나눈 값이며 표본크기와 무관하다. \(\text{Exp}(1)\)에서는 \(\mu_3 = 2\), \(\sigma = 1\), \(\beta_2 = 9\)이므로 \(2/\sqrt 8 = 0.707\)이다.

import numpy as np

rng = np.random.default_rng(0)
B = 200_000

# 정규모집단과 지수모집단에서 (Xbar, S^2)의 상관을 잰다.
# 이론값은 왜도 / sqrt(kappa - 1) 이며 n에 의존하지 않는다.
for name, rvs in [("Normal", lambda s: rng.normal(size=s)),
                  ("Exp(1)", lambda s: rng.exponential(size=s))]:
    for n in (10, 100):
        X = rvs((B, n))
        r = np.corrcoef(X.mean(axis=1), X.var(axis=1, ddof=1))[0, 1]
        print(f"{name:>8}  n = {n:>3}   corr(Xbar, S^2) = {r:+.3f}")

출력:

  Normal  n =  10   corr(Xbar, S^2) = +0.002
  Normal  n = 100   corr(Xbar, S^2) = +0.000
  Exp(1)  n =  10   corr(Xbar, S^2) = +0.701
  Exp(1)  n = 100   corr(Xbar, S^2) = +0.706

정규 쪽은 \(n\)이 무엇이든 0이고, 지수 쪽은 \(n\)이 열 배가 되어도 0.70에 그대로 머문다.

\(t\) 통계량에 미치는 영향. \(t = (\bar X - \mu)/(S/\sqrt n)\)에서 분자가 크면 분모도 함께 커지므로 \(t\)가 그만큼 커지지 못하고, 반대로 \(\bar X\)가 작으면 \(S\)도 작아 \(t\)가 더 음수 쪽으로 간다. 그 결과 \(t\)의 분포가 왼쪽으로 치우친다. 양측검정에서는 두 꼬리의 오차가 어느 정도 상쇄되지만 단측검정에서는 상쇄가 일어나지 않아 문제가 드러난다.

왜 이 문제를 여기에 두는가. \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)과 \(t\) 분포의 유도가 같은 독립성 하나에 함께 기대고 있기 때문이다. 정규성을 가정한다는 말이 실은 여러 겹의 성질을 한꺼번에 빌려 쓴다는 뜻임을 이 상관계수가 보여 준다.

증명. \(\chi^2_d\)의 적률생성함수는 \(M(t) = (1-2t)^{-d/2}\)(\(t < 1/2\))이다. 독립이면 합의 MGF가 곱이므로

이고, 이것이 곧 \(\chi^2_{a+b}\)의 MGF다. MGF가 분포를 유일하게 결정하므로 \(U+V \sim \chi^2_{a+b}\)이다. \(\square\)

정의로 보면 더 분명하다. \(\chi^2_a\)는 독립인 표준정규 \(a\)개의 제곱합이고 \(\chi^2_b\)는 또 다른 \(b\)개의 제곱합이니, 둘을 더하면 \(a+b\)개의 제곱합이다. 자유도는 "더해진 제곱의 개수"이므로 그냥 더해진다.

import numpy as np
from scipy import stats

a, b = 4, 7
U = stats.chi2.rvs(a, size=400_000, random_state=1)
V = stats.chi2.rvs(b, size=400_000, random_state=2)
S = U + V
print(f"  합의 평균 {S.mean():.4f}  분산 {S.var():.4f}")
print(f"  chi2({a+b}) 이론 평균 {a+b}  분산 {2*(a+b)}")
print(f"  KS 검정 p = {stats.kstest(S, 'chi2', args=(a+b,)).pvalue:.4f}")

출력:

  합의 평균 10.9949  분산 21.8936
  chi2(11) 이론 평균 11  분산 22
  KS 검정 p = 0.5618

KS 검정의 \(p = 0.56\)으로 \(\chi^2_{11}\)과 구별되지 않는다.

왜 결정적인가 — 제곱합이 쪼개지기 때문이다. 분산분석의 항등식

에서, 정규성 아래 세 항이 각각 카이제곱을 따르고 자유도도 같은 방식으로 쪼개진다(\(N-1 = (k-1) + (N-k)\)). 가법성이 없다면 "자유도를 나눠 갖는다"는 말 자체가 성립하지 않는다.

여기에 독립성이 더 필요하다. 가법성은 두 조각이 독립일 때만 쓸 수 있는데, 집단 간 제곱합과 집단 내 제곱합이 실제로 독립이라는 것이 코크런 정리의 내용이다. 연습문제 7의 \(\bar X \perp S^2\)이 \(k=1\)인 특수한 경우이며, 11장의 \(F\) 통계량이 이 두 조각의 비로 만들어진다.

뒤집으면 자유도를 세는 규칙이 된다. 제곱합을 직교하는 조각으로 나눌 때마다 자유도가 그 조각의 차원만큼 배분된다. 회귀에서 잔차 자유도가 \(n - p\)인 것도 같은 셈법이며, 여기서 잃는 자유도 하나하나가 "추정한 모수 하나"에 대응한다(연습문제 2).

import numpy as np
from scipy import stats

print(f"{'d':>5}{'정규근사 오차':>16}{'WH 오차':>12}")
for d in (5, 10, 20, 50, 100):
    q = stats.chi2.ppf([0.05, 0.95], d)                  # 참 분위수
    simp = d + stats.norm.ppf([0.05, 0.95]) * np.sqrt(2 * d)
    m, v = 1 - 2 / (9 * d), 2 / (9 * d)
    wh = d * (stats.norm.ppf([0.05, 0.95]) * np.sqrt(v) + m) ** 3
    print(f"{d:>5}{np.abs(simp-q).max():>16.4f}{np.abs(wh-q).max():>12.4f}")

출력:

    d         정규근사 오차       WH 오차
    5          1.3470      0.0266
   10          1.2963      0.0153
   20          1.2538      0.0087
   50          1.2128      0.0043
  100          1.1912      0.0025

차이가 압도적이다. \(d = 100\)에서도 단순 정규근사는 분위수를 \(1.19\)나 틀리는데, 윌슨–힐퍼티는 \(0.0025\)로 \(476\)배 정확하다.

단순 근사가 좀처럼 나아지지 않는 것에 주목하라. \(d\)를 \(5\)에서 \(100\)으로 \(20\)배 키워도 오차가 \(1.35\)에서 \(1.19\)로 거의 줄지 않는다. 카이제곱의 왜도가 \(\sqrt{8/d}\)로 천천히 줄기 때문이며, 정규근사는 왜도를 아예 \(0\)으로 놓는다.

세제곱근이 하는 일이 그 왜도를 죽이는 것이다. \(\chi^2\)는 오른쪽 꼬리가 긴데, 세제곱근은 큰 값을 강하게 눌러 분포를 대칭에 가깝게 만든다. 로그를 취하면 오히려 왼쪽으로 지나치게 당겨지고, 제곱근은 덜 당긴다. \(1/3\) 제곱이 그 중간에서 왜도를 거의 정확히 상쇄한다.

어디에 쓰이는가. 카이제곱 분위수표가 없던 시절의 계산 도구였고, 지금도 자유도가 매우 클 때 유용하다. 유전체 분석처럼 자유도가 수만인 검정에서 scipy의 chi2.sf가 수치적으로 불안정해질 수 있는데, 변환을 거치면 표준정규 꼬리 계산으로 바뀌어 안정적이다.

일반적인 교훈이 하나 있다. 왜도가 남은 분포를 정규로 근사할 때는 먼저 변환해 왜도를 없앤 뒤 근사하는 것이 훨씬 낫다. 분산의 신뢰구간을 로그 척도에서 만드는 것, 상관계수에 피셔 \(z\) 변환을 쓰는 것이 모두 같은 발상이다.

불편이 아니다. \(\sqrt{\cdot}\)는 오목함수이므로 옌센 부등식에 의해

이다. 즉 \(S\)는 \(\sigma\)를 과소추정한다. 등호는 \(S^2\)이 상수일 때만 성립하는데, 그럴 리가 없다.

정확한 상수를 구한다. \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)이므로 \(S = \sigma\sqrt{\chi^2_{n-1}/(n-1)}\)이고, \(\chi^2_d\)의 제곱근의 기댓값이

이므로 \(d = n-1\)을 넣으면

를 얻는다. \(\square\)

import numpy as np
from scipy import special

rng = np.random.default_rng(0)
print(f"{'n':>5}{'c4':>10}{'편향 %':>10}{'모의 E[S]/sigma':>18}")
for n in (2, 5, 10, 30, 100):
    c4 = np.sqrt(2 / (n - 1)) * special.gamma(n / 2) / special.gamma((n - 1) / 2)
    S = rng.normal(0, 1, (300_000, n)).std(axis=1, ddof=1)
    print(f"{n:>5}{c4:>10.6f}{(c4-1)*100:>10.3f}{S.mean():>18.6f}")

출력:

    n        c4      편향 %     모의 E[S]/sigma
    2  0.797885   -20.212          0.800341
    5  0.939986    -6.001          0.938909
   10  0.972659    -2.734          0.972315
   30  0.991418    -0.858          0.991022
  100  0.997478    -0.252          0.997573

\(n = 2\)에서 \(20\%\)나 과소추정한다. \(n = 10\)이면 \(2.7\%\), \(n = 30\)이면 \(0.9\%\)로 줄고, \(c_4 \to 1\)이므로 점근적으로는 불편이다.

"불편성은 변환을 따라가지 않는다"는 것이 핵심이다. \(S^2\)이 불편이라고 해서 \(\sqrt{S^2}\)이 불편일 이유가 없다. 일반적으로 \(g\)가 비선형이면 \(E[g(T)] \ne g(E[T])\)이며, 정규분포 쪽 연습문제 10의 델타 방법이 그 차이를 근사하는 도구였다.

실무에서 보정하는 자리가 있다. 품질관리의 관리도에서 표본별 \(s\)를 평균내어 \(\sigma\)를 추정할 때 \(\bar s / c_4\)로 나눈다. \(c_4\) 표가 관리도 교재에 실려 있는 이유가 이것이며, 부분군 크기가 \(4\sim5\)로 작아 편향이 \(6\%\) 수준이라 무시할 수 없기 때문이다.

그런데 대부분의 통계 추론에서는 보정하지 않는다. \(t\) 통계량 \((\bar X-\mu)/(S/\sqrt n)\)은 \(S\)가 편향된 채로도 정확히 \(t_{n-1}\)을 따른다. 분포를 유도할 때 편향이 이미 셈에 들어가 있기 때문이다. 불편성이 그 자체로 목표는 아니라는 점을 보여 주는 좋은 예이며, 6장에서 불편성과 평균제곱오차를 견줄 때 다시 만난다. \(\square\)