유도와 분포이론¶

1단계. 검정통계량을 계산한다.

2단계. \(\alpha = 0.05\)(양측)에서 \(\chi^2_{20}\)의 임계값을 구한다.

• 하단 임계값: \(\chi^2_{0.975, 20} = 9.591\)
• 상단 임계값: \(\chi^2_{0.025, 20} = 34.170\)

3단계. \(9.591 < 24.667 < 34.170\)이므로 검정통계량이 채택역 안에 있다. \(H_0\)을 기각하지 못한다.

유의수준 5%에서 모분산이 15와 다르다고 결론지을 만한 증거가 충분하지 않다(\(p = 0.429\)).

(1) 임계값을 \(S^2\) 의 눈금으로 되돌린다. 채택역은 \(\chi^2_{\alpha/2,\,n-1} \le (n-1)S^2/\sigma_0^2 \le \chi^2_{1-\alpha/2,\,n-1}\) 이므로, 양변에 \(\sigma_0^2/(n-1)\) 을 곱하면

이다. 채택역이 \(\sigma_0^2\) 의 상수배 구간이고 그 두 배수 \(\chi^2_{\alpha/2}/(n-1)\), \(\chi^2_{1-\alpha/2}/(n-1)\) 은 자유도와 \(\alpha\) 만으로 정해진다. \(n = 21\), \(\alpha = 0.05\) 에서 \(\chi^2_{0.025,\,20} = 9.5908\), \(\chi^2_{0.975,\,20} = 34.1696\) 이므로 배수가 \(0.47954\) 와 \(1.70848\) 이고

이면 기각하지 못한다. 관측된 \(18.5\) 는 이 안쪽이고 위쪽 경계까지 \(7.1272\) 가 남았다. 표본분산이 \(38\%\) 더 커졌어야 기각했다는 뜻이다(\(25.6272/18.5 = 1.385\)).

신뢰구간과의 맞물림. 앞 절의 쌍대성이 여기서 그대로 보인다. \(S^2 = 18.5\) 의 \(95\%\) 신뢰구간은

이고 \(\sigma_0^2 = 15\) 가 그 안에 있다. 두 진술은 같은 부등식을 어느 쪽 변수로 풀었느냐의 차이일 뿐이다. 채택역은 \(\sigma_0^2\) 을 고정하고 \(S^2\) 에 대해 푼 것이고, 신뢰구간은 \(S^2\) 을 고정하고 \(\sigma^2\) 에 대해 푼 것이다.

확인한다.

import numpy as np
from scipy import stats

n = 21
s_squared = 18.5
sigma_0_squared = 15
alpha = 0.05

# (n-1)S^2/sigma_0^2 은 귀무가설이 참일 때 자유도 n-1 인 카이제곱을 따른다.
# 이 검정은 정규성에 매우 민감하다 — 평균에 대한 t 검정과 달리 중심극한정리의
# 보호를 받지 못하기 때문이다.
chi2_stat = (n - 1) * s_squared / sigma_0_squared

# 좌우가 대칭이 아니므로 양쪽 기각값을 따로 구한다.
df = n - 1
chi2_lower = stats.chi2.ppf(alpha / 2, df)
chi2_upper = stats.chi2.ppf(1 - alpha / 2, df)

# 양측 p-값은 두 꼬리 넓이 중 작은 쪽을 두 배 한다.
p_value = 2 * min(stats.chi2.cdf(chi2_stat, df), stats.chi2.sf(chi2_stat, df))

print(f"Test statistic: {chi2_stat:.3f}")
print(f"Critical values: [{chi2_lower:.3f}, {chi2_upper:.3f}]")
print(f"P-value: {p_value:.4f}")

if p_value < alpha:
    print("Reject H0: variance differs from the hypothesized value.")
else:
    print("Fail to reject H0: insufficient evidence of a difference.")

# (1) 어떤 S^2 이었다면 기각했을까. 임계값을 S^2 의 눈금으로 되돌린다.
s2_lo = sigma_0_squared * chi2_lower / df
s2_hi = sigma_0_squared * chi2_upper / df
print(f"\nS^2 의 채택역 = [{s2_lo:.4f}, {s2_hi:.4f}]")
print(f"  관측된 S^2 = {s_squared}  ->  위쪽 경계까지 {s2_hi - s_squared:.4f} 남았다")
print(f"  배수로 보면 sigma_0^2 x [{chi2_lower / df:.5f}, {chi2_upper / df:.5f}]")

# 쌍대 관계: sigma_0^2 = 15 가 S^2 = 18.5 의 신뢰구간 안에 있는가.
ci_lo = df * s_squared / chi2_upper
ci_hi = df * s_squared / chi2_lower
print(f"\nS^2 = {s_squared} 의 95% 신뢰구간 = ({ci_lo:.4f}, {ci_hi:.4f})")
print(f"  sigma_0^2 = {sigma_0_squared} 이 구간 안에 있는가: {ci_lo <= sigma_0_squared <= ci_hi}")

출력:

Test statistic: 24.667
Critical values: [9.591, 34.170]
P-value: 0.4290
Fail to reject H0: insufficient evidence of a difference.

S^2 의 채택역 = [7.1931, 25.6272]
  관측된 S^2 = 18.5  ->  위쪽 경계까지 7.1272 남았다
  배수로 보면 sigma_0^2 x [0.47954, 1.70848]

S^2 = 18.5 의 95% 신뢰구간 = (10.8283, 38.5787)
  sigma_0^2 = 15 이 구간 안에 있는가: True

유도한 채택역 \([7.1931, 25.6272]\) 와 신뢰구간 \((10.8283, 38.5787)\) 이 코드와 맞고, \(\sigma_0^2 = 15\) 가 구간 안이라는 판정이 \(p = 0.4290 > 0.05\) 와 일치한다.

(2) 이 모든 수가 한 줄에 매달려 있다. 위에서 쓴 것은 \(\chi^2_{20}\) 분포뿐이고, 그것은 이 쪽이 유도한 \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\) 에서 나왔으며, 그 유도는 첫 단계부터 정규성을 쓴다. \(\sum Z_i^2\) 이 카이제곱이 되려면 \(Z_i\) 가 정규여야 하기 때문이다. 정규성이 깨지면 임계값 \(9.591\) 과 \(34.170\) 이 엉뚱한 자리가 된다.

얼마나 엉뚱해지는지는 미리 계산할 수 있다. 5.3절이 유도한 결과를 빌린다. \(n\) 이 크면 \(\chi^2_{n-1}\) 의 양측 임계값이 \((n-1)\bigl(1 \pm 1.96\sqrt{2/(n-1)}\bigr)\) 이므로, 검정은 \(S^2/\sigma^2\) 의 표준편차를 \(\sqrt{2/n}\) 로 믿고 기각역을 긋는다. 그런데 참 표준편차는 \(\operatorname{Var}(S^2) \approx (\beta_2-1)\sigma^4/n\) 에서 \(\sqrt{(\beta_2-1)/n}\) 이다. 기각역이 참 표준편차 단위로는

만큼만 뻗으므로

가 된다. \(\sqrt n\) 이 분자와 분모에서 함께 약분되어 \(n\) 이 사라졌다. 임계값도 \(1/\sqrt n\) 로 줄고 참 산포도 \(1/\sqrt n\) 로 줄기 때문이다. 그러므로 이 어긋남은 표본을 키워도 낫지 않는다. 남는 것은 모집단 첨도 \(\beta_2\) 하나뿐이다.

모의실험으로 확인한다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(15)

pops = [
    ("정규", lambda s: rng.normal(size=s), 1.0, 3.0),
    ("균등", lambda s: rng.uniform(size=s), 1 / 12, 1.8),
    ("지수", lambda s: rng.exponential(size=s), 1.0, 9.0),
    ("로그정규", lambda s: rng.lognormal(0, 1, size=s),
     (np.exp(1) - 1) * np.exp(1), np.exp(4) + 2 * np.exp(3) + 3 * np.exp(2) - 3),
]

ns = [21, 100, 1000, 10_000]
hdr = "".join(f"{'n=' + str(n):>10}" for n in ns)
print(f"{'모집단':>8}{'beta2':>10}{'극한':>9}{hdr}")
for name, gen, var0, beta2 in pops:
    limit = 2 * (1 - stats.norm.cdf(1.96 * np.sqrt(2 / (beta2 - 1))))
    row = ""
    for n in ns:
        lo, hi = stats.chi2.ppf([0.025, 0.975], n - 1)
        reps = 40_000 if n <= 1000 else 4_000
        rej = 0
        done = 0
        while done < reps:
            blk = min(2000, reps - done)
            T = (n - 1) * gen((blk, n)).var(axis=1, ddof=1) / var0
            rej += np.sum((T < lo) | (T > hi))
            done += blk
        row += f"{rej / reps:>10.4f}"
    print(f"{name:>6}{beta2:>11.3f}{limit:>9.4f}{row}")
print("\n반복수: n <= 1000 은 40,000 벌, n = 10,000 은 4,000 벌")
print(f"  0.05 근처 몬테카를로 표준오차:  40,000 벌 {np.sqrt(0.05 * 0.95 / 40_000):.4f}"
      f"   4,000 벌 {np.sqrt(0.05 * 0.95 / 4_000):.4f}")

출력:

     모집단     beta2       극한      n=21     n=100    n=1000   n=10000
    정규      3.000   0.0500    0.0505    0.0485    0.0499    0.0548
    균등      1.800   0.0019    0.0040    0.0020    0.0020    0.0013
    지수      9.000   0.3271    0.2718    0.3087    0.3212    0.3222
  로그정규    113.936   0.7942    0.5763    0.6553    0.7208    0.7560

반복수: n <= 1000 은 40,000 벌, n = 10,000 은 4,000 벌
  0.05 근처 몬테카를로 표준오차:  40,000 벌 0.0011   4,000 벌 0.0034

정규 줄만 수평이다. \(0.0505\), \(0.0485\), \(0.0499\), \(0.0548\) 로 \(n\) 이 \(500\) 배 늘어나는 동안 \(0.05\) 에 머문다. 마지막 칸이 조금 높아 보이지만 반복이 \(4{,}000\) 벌뿐이라 몬테카를로 표준오차가 \(0.0034\) 이고, \(0.0548\) 은 \(0.05\) 에서 \(1.4\) 표준오차 떨어진 값이다. 우연으로 설명된다.

나머지 셋은 각자의 상수로 간다. 지수분포는 \(0.2718 \to 0.3087 \to 0.3212 \to 0.3222\) 로 극한값 \(0.3271\) 에 붙어 간다. 명목 \(5\%\) 검정이 실제로는 세 번에 한 번 꼴로 기각한다. 균등분포는 반대 방향으로 \(0.002\) 에 머물러 지나치게 보수적이다(\(\beta_2 < 3\) 이라 참 산포가 믿는 것보다 작다). 로그정규는 \(0.5763 \to 0.6553 \to 0.7208 \to 0.7560\) 으로 \(0.7942\) 를 향해 계속 올라간다.

로그정규가 아직 극한에 닿지 않은 것은 흠이 아니라 예상된 일이다. 위 유도는 \(S^2\) 의 정규근사에 기대는데, \(\beta_2 = 113.9\) 인 모집단에서 그 근사가 쓸 만해지려면 \(n\) 이 아주 커야 한다. 중요한 것은 방향이다. 네 줄 가운데 \(n\) 과 함께 \(0.05\) 쪽으로 돌아오는 줄은 하나도 없다. 정규는 처음부터 맞고, 나머지 셋은 \(n\) 을 키울수록 명목값에서 더 멀어진다.

이것이 분산 추론과 평균 추론의 갈림길이다. 평균에 대한 \(t\) 검정은 중심극한정리가 비정규성을 씻어 주므로 \(n\) 을 키우면 명목 수준에 수렴한다. 분산에는 그런 보호가 없다. 유도의 첫 단계가 정규성이고, 그 단계가 틀리면 표본크기로는 고칠 수 없다. 연습문제 4 가 "평균은 맞지만 산포가 틀린다"고 적은 것의 수치적 실물이 이 표다.

\(X_i - \mu = (X_i - \bar{X}) + (\bar{X} - \mu)\)로 쓰고 제곱하여 합한다.

두 번째 항에서 \((\bar{X} - \mu)\)는 \(i\)에 의존하지 않으므로 합 밖으로 빼냈다. 그런데

이므로 교차항이 사라진다. 세 번째 항은 \(i\)에 무관한 값을 \(n\)번 더한 것이므로 \(n(\bar{X}-\mu)^2\)이다. 따라서

이 항등식은 평행축 정리 또는 편차제곱합 분해로 불린다. 교차항이 0이 되는 것은 \(\bar{X}\)가 \(\sum_i (X_i - c)^2\)을 최소화하는 \(c\)라는 사실의 다른 표현이며, 이 때문에 \(\bar{X}\)가 최소제곱 추정량이 된다.

부수적 결과로 \(\mu \neq \bar{X}\)인 한 \(\sum_i (X_i - \bar{X})^2 < \sum_i (X_i - \mu)^2\)이다. 표본평균으로 계산한 제곱합이 참 평균으로 계산한 것보다 항상 작으며, 이것이 \(n\)이 아니라 \(n-1\)로 나누어 편향을 보정하는 이유이다. \(\square\)

정규성 아래의 증명 개요. \(Z_i = (X_i - \mu)/\sigma\)라 하고 벡터 \(\mathbf{Z} = (Z_1,\ldots,Z_n)^\top \sim \mathcal{N}(\mathbf{0}, I_n)\)을 생각한다. 첫 행이 \(\frac{1}{\sqrt{n}}(1,1,\ldots,1)\)인 직교행렬 \(Q\)를 잡고 \(\mathbf{Y} = Q\mathbf{Z}\)라 하자.

직교변환은 표준정규 벡터의 분포를 보존하므로 \(\mathbf{Y} \sim \mathcal{N}(\mathbf{0}, I_n)\)이고 성분들이 서로 독립이다. 그런데

\(Y_1\)과 \((Y_2,\ldots,Y_n)\)이 독립이므로 \(\bar{X}\)와 \(S^2\)도 독립이다. 동시에 \(\sum_{j=2}^n Y_j^2\)은 독립인 표준정규 제곱 \(n-1\)개의 합이므로 \(\chi^2_{n-1}\)을 따른다. 이것이 본문 결과의 엄밀한 증명이다.

왜 정규분포에만 특유한가. 증명에서 결정적으로 쓴 성질은 표준정규 벡터의 분포가 직교변환에 불변이라는 것이다. 이 회전불변성을 갖는 유일한 독립 성분 분포가 정규분포이다(Maxwell의 정리).

실제로 Lukacs의 정리는 더 강한 결과를 준다. \(X_1,\ldots,X_n\)이 i.i.d.일 때 \(\bar{X}\)와 \(S^2\)이 독립일 필요충분조건이 바탕 분포가 정규라는 것이다.

실무적 함의. 비정규 자료에서는 \(\bar{X}\)와 \(S^2\)이 상관되어 있다. 예컨대 오른쪽으로 치우친 분포에서는 큰 \(\bar{X}\)가 큰 \(S^2\)과 함께 나타나는 경향이 있다. 이 상관이 \(t\) 통계량의 분포를 왜곡하는 한 원인이며, 분산 검정에서는 카이제곱 근사를 무너뜨린다. \(\square\)

비대칭성의 근원. \(\chi^2_\nu\)의 왜도는 \(\sqrt{8/\nu}\)이다.

작은 자유도에서 분포가 강하게 오른쪽으로 치우쳐 있으므로 상단 분위수가 평균 \(\nu\)에서 하단 분위수보다 훨씬 멀리 떨어져 있다. 신뢰구간을 만들 때 \(\sigma^2\)을 이 분위수로 나누므로 결과 구간의 상한이 하한보다 점추정값에서 멀어진다.

\(\nu \to \infty\)이면 왜도가 \(0\)으로 가고 분포가 정규에 접근하여 구간이 대칭에 가까워진다.

근사의 정확도. 두 가지 근사를 비교하자.

• 단순 정규근사: \(\chi^2_\nu \approx \nu + z\sqrt{2\nu}\)
• Wilson-Hilferty 근사: \(\chi^2_\nu \approx \nu\left(1 - \frac{2}{9\nu} + z\sqrt{\frac{2}{9\nu}}\right)^3\)

import numpy as np
from scipy import stats

z = 1.96
print(f"{'df':>5} {'exact':>10} {'normal':>10} {'Wilson-H':>10}")
for df in [5, 20, 50, 200]:
    exact = stats.chi2.ppf(0.975, df)
    norm_ap = df + z * np.sqrt(2 * df)
    wh = df * (1 - 2 / (9 * df) + z * np.sqrt(2 / (9 * df)))**3
    print(f"{df:>5} {exact:>10.3f} {norm_ap:>10.3f} {wh:>10.3f}")

출력:

   df      exact     normal   Wilson-H
    5     12.833     11.198     12.822
   20     34.170     32.396     34.172
   50     71.420     69.600     71.424
  200    241.058    239.200    241.061

단순 정규근사는 \(\nu = 5\)에서 12.8을 11.2로 추정하여 13% 오차를 낸다. \(\nu = 200\)에서도 여전히 0.8% 부족하다. 근사가 치우침을 무시하므로 상단 분위수를 체계적으로 과소추정한다.

Wilson-Hilferty 근사는 \(\nu = 5\)에서 이미 소수 둘째 자리까지 맞고 \(\nu \geq 20\)에서는 사실상 정확하다. \((\chi^2_\nu/\nu)^{1/3}\)이 거의 정규라는 사실에 기반하며, 세제곱근 변환이 치우침을 제거한다. Box-Cox 계열의 변환이 왜 유용한지를 보여주는 좋은 예이다. \(\square\)

평균의 계산(정규성 불필요). \(X_i\)가 i.i.d.이고 \(\mathbb{E}[X_i] = \mu\), \(\operatorname{Var}(X_i) = \sigma^2\)이라고만 가정하자. 연습문제 1의 항등식에 기댓값을 취하면

따라서

이며, 이는 \(\chi^2_{n-1}\)의 평균과 일치한다. 어떤 분포에서든 성립한다. 동시에 \(\mathbb{E}[S^2] = \sigma^2\), 곧 \(S^2\)이 \(\sigma^2\)의 불편추정량임을 보여준다.

평균만으로는 부족한 이유. 검정을 하려면 통계량의 전체 분포가 필요하지 평균만으로는 안 된다. 두 분포가 같은 평균을 가져도 임계값이 전혀 다를 수 있다.

구체적으로 15.1절에서 보았듯

인데 \(\chi^2_{n-1}\)의 분산은 \(2(n-1)\)이다. 두 값이 일치하는 것은 \(\gamma_2 = 0\), 곧 정규분포일 때뿐이다.

곧 비정규 자료에서도 검정통계량의 평균은 맞지만 산포가 틀린다. 기준분포의 중심은 제자리인데 폭이 좁으므로, 실제 통계량이 임계값 밖으로 자주 벗어나 제1종 오류가 팽창한다.

이것이 Cochran 정리가 필요한 이유이다. 정규성 아래에서만 1차 적률뿐 아니라 분포 전체가 정확히 카이제곱이 된다. \(\square\)