S²의 표본분포 (Normal)¶

(1) 이론값. \(E[S^2] = \sigma^2\)은 정규성 없이 성립하므로 넷 다 모집단 분산을 그대로 겨냥한다. 틀릴 수 있는 것은 폭뿐이고, 그것은 위의 식

에 \(n = 100\)을 넣어 \((\beta_2 - 97/99)\sigma^4/100\)으로 계산한다. \(\chi^2(2) = \text{Exp}(1/2)\)라 지수와 모양이 같아 첨도도 \(9\)로 같고, 분산만 \(4\)배다.

분산 배율은 \(\frac{n-1}{2n}\left(\beta_2 - \frac{n-3}{n-1}\right)\)이고 \(n = 100\)에서 \(0.495\,(\beta_2 - 0.9798)\)이다. 첫 줄이 정확히 \(1\)이 되는 것이 정리 1이고, 나머지 세 줄이 \(1\)에서 벗어난 양이 이 쪽의 주제다. \(n = 100\)인데도 지수·카이제곱에서 폭이 두 배 가까이 넓고 균등에서 \(0.64\)배로 좁다.

(2) 모의실험.

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

np.random.seed(1)

n_population = 10_000
n_sample = 100
n_sim = 1_000

# 모집단 넷을 준비한다. 정규 하나와 정규가 아닌 셋이다.
# 카이제곱 결과는 정규모집단에서만 정확하므로,
# 나머지 셋에서 어긋나는 모습을 보는 것이 이 그림의 목적이다.
# random_state를 각각 다르게 주어 네 모집단이 서로 무관하게 만든다.
populations = {
    "Normal(0,1)":  stats.norm().rvs(n_population, random_state=1),
    "Exp(1)":       stats.expon().rvs(n_population, random_state=2),
    "Chi-sq(2)":    stats.chi2(df=2).rvs(n_population, random_state=3),
    "Uniform(0,1)": stats.uniform().rvs(n_population, random_state=4),
}

# 명목 모집단의 분산과 첨도. 실현된 1만 개는 이 값에서 조금 비껴난다.
nominal = {"Normal(0,1)": (1.0, 3.0), "Exp(1)": (1.0, 9.0),
           "Chi-sq(2)": (4.0, 9.0), "Uniform(0,1)": (1 / 12, 1.8)}

def sd_s2(sigma2, beta2, n):
    """Var(S^2) = (beta2 - (n-3)/(n-1)) sigma^4 / n 의 제곱근."""
    return np.sqrt((beta2 - (n - 3) / (n - 1)) * sigma2 ** 2 / n)

fig, axes = plt.subplots(1, len(populations), figsize=(16, 3.5))

print("모집단          beta2 명목/실현   E[S^2] 실현/모의    sd(S^2) 명목/실현/모의   배율 (모의/카이제곱)^2  이론 배율")
for ax, (name, population) in zip(axes, populations.items()):
    # 표본분산 S^2 의 표집분포를 모의실험으로 얻는다.
    s2_sims = np.array([
        np.random.choice(population, size=n_sample, replace=False).var(ddof=1)
        for _ in range(n_sim)
    ])

    # 명목 모집단, 실현된 1만 개, 모의실험의 세 층을 나란히 적는다.
    s2_nom, b2_nom = nominal[name]
    s2_real, b2_real = population.var(), stats.kurtosis(population) + 3
    # 카이제곱이 예측하는 폭. 정규모집단에서만 맞는 값이다.
    chi2_sd = np.sqrt(2 / (n_sample - 1)) * s2_real
    ratio_sim = (s2_sims.std(ddof=1) / chi2_sd) ** 2
    ratio_thm = (n_sample - 1) / (2 * n_sample) * (b2_real - (n_sample - 3) / (n_sample - 1))
    print(f"{name:14s} {b2_nom:5.2f} / {b2_real:5.3f}   {s2_real:7.5f} / {s2_sims.mean():7.5f}   "
          f"{sd_s2(s2_nom, b2_nom, n_sample):.5f} / {sd_s2(s2_real, b2_real, n_sample):.5f} / {s2_sims.std(ddof=1):.5f}"
          f"        {ratio_sim:5.3f}            {ratio_thm:5.3f}")

    # 모의실험으로 얻은 값들의 히스토그램.
    _, bins, _ = ax.hist(s2_sims, density=True, bins=30,
                         alpha=0.5, edgecolor="white",
                         label=r"simulated $S^2$")

    # 이론적 카이제곱 밀도를 S^2 의 눈금으로 옮겨 그린다.
    # 정리는 (n-1)S^2/sigma^2 ~ chi^2(n-1) 이므로
    # S^2 = (sigma^2/(n-1)) * chi^2 = X/c  (단, c = (n-1)/sigma^2) 이다.
    #
    # 변수변환 Y = X/c 의 밀도는 f_Y(y) = f_X(cy) * c 다.
    # 마지막에 곱하는 c 가 그 야코비안이며, 이것을 빠뜨리면
    # 곡선의 넓이가 1이 되지 않아 히스토그램과 눈금이 어긋난다.
    df = n_sample - 1
    sigma2 = population.var()
    c = df / sigma2
    x_grid = np.linspace(bins[0], bins[-1], 300)
    pdf = stats.chi2(df).pdf(x_grid * c) * c
    ax.plot(x_grid, pdf, "--r", lw=2, alpha=0.7,
            label=r"$\chi^2$-based PDF")

    ax.set_title(name)
    ax.set_xlabel(r"$S^2$")

axes[0].set_ylabel("Density")
axes[-1].legend(fontsize=8)
plt.tight_layout()
plt.show()

출력:

모집단          beta2 명목/실현   E[S^2] 실현/모의    sd(S^2) 명목/실현/모의   배율 (모의/카이제곱)^2  이론 배율
Normal(0,1)     3.00 / 3.029   0.99757 / 1.00183   0.14213 / 0.14279 / 0.14371        1.027            1.014
Exp(1)          9.00 / 8.166   0.98539 / 0.97830   0.28320 / 0.26415 / 0.26405        3.554            3.557
Chi-sq(2)       9.00 / 8.242   3.86995 / 3.86573   1.13280 / 1.04290 / 1.04117        3.583            3.595
Uniform(0,1)    1.80 / 1.785   0.08472 / 0.08441   0.00755 / 0.00760 / 0.00776        0.415            0.398

불편성은 넷 다 맞는다. 모의 \(E[S^2]\)이 실현 모집단의 분산에서 떨어진 거리가 몬테카를로 오차(\(\operatorname{sd}(S^2)/\sqrt{1000}\)) 단위로 각각 \(0.9\)배, \(0.9\)배, \(0.1\)배, \(1.3\)배다. 정규든 아니든, 첨도가 \(1.8\)이든 \(9\)이든 중심은 흔들리지 않는다.

폭은 실현 이론값과 맞는다. 셋째 열의 둘째 수(실현)와 셋째 수(모의)를 견주면 상대차가 \(+0.6\%\), \(-0.04\%\), \(-0.17\%\), \(+2.1\%\)다. 되풀이 1천 번에서 표준편차의 몬테카를로 요동이 \(1/\sqrt{2000} = 2.2\%\)쯤이므로 넷 다 그 안에 든다. 마지막 두 열의 배율도 \(1.027\) 대 \(1.014\), \(3.554\) 대 \(3.557\), \(3.583\) 대 \(3.595\), \(0.415\) 대 \(0.398\)로 맞아 떨어진다. 배율은 표준편차의 제곱이라 상대오차가 두 배로 보이는 것을 감안해야 한다.

그런데 명목 이론값과는 어긋난다. 지수에서 (1)이 예측한 \(\operatorname{sd}(S^2) = 0.28320\)과 배율 \(3.970\)이 모의에서는 \(0.26405\)와 \(3.554\)로 나왔다. 둘째 열이 범인을 가리킨다. \(\text{Exp}(1)\)의 첨도는 \(9\)인데 거기서 뽑힌 1만 개의 첨도는 \(8.166\)밖에 되지 않는다. 카이제곱 모집단에서도 \(9\)가 \(8.242\)로 내려갔다. 첨도는 4차 적률이라 유한한 모집단이 꼬리를 다 담지 못하면 곧바로 작아진다. 실현된 \(\beta_2\)를 식에 넣으면 \(0.26415\)와 \(3.557\)이 나와 모의값과 소수 셋째 자리까지 맞는다.

꼬리가 가벼운 쪽에서는 이 문제가 없다. 균등의 첨도는 \(1.8\)에서 \(1.785\)로 거의 그대로고, 정규도 \(3.0\)에서 \(3.029\)다. 어긋남이 생기는 곳은 꼬리가 무거운 모집단뿐이며, 그것은 식이 틀린 것이 아니라 "1만 개를 \(\text{Exp}(1)\) 자체로 여긴" 가정이 틀린 것이다.

그림의 네 패널을 보면 히스토그램은 넷 다 종 모양에 가깝다. \(S^2\)에도 중심극한정리가 작용하기 때문이다. 그런데 붉은 카이제곱 곡선과 겹치는 것은 첫째 패널뿐이고, 둘째·셋째는 곡선보다 낮고 넓으며 넷째는 높고 좁다. 모양이 닮아 가는 일과 폭이 맞는 일이 따로 논다는 것이 숫자로도, 그림으로도 확인된다.

여기서 \(\sigma^2 = 4\), \(n = 10\)이므로 \(\frac{(n-1)S^2}{\sigma^2} = \frac{9S^2}{4} \sim \chi^2(9)\)이다.

구해야 할 것은:

Python으로:

from scipy import stats
p = 1 - stats.chi2.cdf(13.5, df=9)
print(f"P(S^2 > 6) = {p:.4f}")

출력:

P(S^2 > 6) = 0.1413

결과는 \(P(\chi^2(9) > 13.5) \approx 0.1413\)이다. \(\square\)

정의에서 출발한다:

합을 전개하면:

기댓값을 취하면:

따라서:

\(n - 1\)로 나누면:

\(\square\)

정리에 의해 \(\frac{(n-1)S^2}{\sigma^2} \sim \chi^2(n-1)\)이다.

\(Q = \frac{(n-1)S^2}{\sigma^2}\)라 하면 \(S^2 = \frac{\sigma^2 Q}{n-1}\)이므로:

\(\square\)

\(\text{Exp}(1)\)에서 \(\sigma^2 = 1\), \(\beta_2 = 9\)이다.

\(S^2\)의 분산에 대한 일반 공식은:

\(n = 100\)이면:

정규이론 값은:

지수 모집단에서 \(S^2\)의 분산이 정규이론 값보다 약 4배 크다. 지수분포의 두꺼운 꼬리(초과첨도 \(= 6\)) 때문에 \(S^2\)의 변동성이 커지기 때문이다. \(\square\)

\((n-1)S^2/\sigma^2 \sim \chi^2(n-1)\)이므로:

\(\sigma^2\)에 대해 정리하면:

카이제곱 분포표(또는 Python)에서:

from scipy import stats

# 자유도 24인 카이제곱 분포의 양쪽 2.5% 지점.
# 정규분포와 달리 두 값이 0을 중심으로 대칭이 아니다.
lower_chi2 = stats.chi2.ppf(0.025, df=24)
upper_chi2 = stats.chi2.ppf(0.975, df=24)
print(f"chi2_0.025(24) = {lower_chi2:.4f}")
print(f"chi2_0.975(24) = {upper_chi2:.4f}")

# 신뢰구간. 큰 분위수가 **아래쪽** 한계를 만든다는 점에 주의하라.
# sigma^2 이 분모에 있으므로 부등식을 뒤집으면 순서가 바뀐다.
S2, n = 12, 25
print(f"95% CI = ({(n-1)*S2/upper_chi2:.2f}, {(n-1)*S2/lower_chi2:.2f})")

출력:

chi2_0.025(24) = 12.4012
chi2_0.975(24) = 39.3641
95% CI = (7.32, 23.22)

\(\sigma^2\)에 대한 95% 신뢰구간은 약 \((7.32, 23.22)\)이다. 카이제곱 분포가 치우쳐 있으므로 이 구간이 \(S^2 = 12\)를 중심으로 대칭이 아님에 유의하라. \(\square\)

검정. \(\dfrac{S_1^2/\sigma_1^2}{S_2^2/\sigma_2^2} \sim F_{n_1-1,\,n_2-1}\)이므로 \(H_0:\sigma_1^2=\sigma_2^2\) 아래에서

양측 5%의 기각역은 \(F < F_{11,14,0.025} = 0.298\) 또는 \(F > F_{11,14,0.975} = 3.095\)이다. \(1.546\)은 그 사이에 있으므로 기각하지 못한다. p-값은 0.438이다.

신뢰구간. 위 추축량을 \(\sigma_1^2/\sigma_2^2\)에 대해 풀면

이다.

구간이 대단히 넓다. 참 분산비가 0.5일 수도 5일 수도 있다. \(n\)이 각각 12, 15로 작아 분산 추정 자체가 부정확하기 때문이다(앞서 본 대로 분산을 10% 이내로 알려면 200개가 필요하다).

실무적 경고. 이 \(F\) 검정은 정규성에 극도로 민감하다. 모집단이 조금만 꼬리가 두꺼워도 실제 오류율이 명목 5%를 크게 웃돈다. 게다가 "합동 \(t\) 검정을 쓸지 정하려고 먼저 등분산을 검정한다"는 이단계 절차는 권장되지 않는다. 검정력이 낮아 실제 차이를 놓치기 쉽고, 두 검정이 겹쳐 전체 오류율이 통제되지 않기 때문이다. 처음부터 웰치 검정을 쓰는 것이 낫다.

착상. 두 검정 모두 산포 자체를 다루는 대신 산포를 나타내는 새 변수를 만들어 평균을 비교한다.

각 관측값을 중심으로부터의 절대편차로 바꾼다.

그다음 이 \(z\) 값들에 대해 보통의 분산분석(또는 \(t\) 검정) 을 한다. 집단마다 산포가 다르면 \(z\)의 평균이 달라지므로 검출된다.

• 르빈 검정: \(c_i\)로 집단 평균 \(\bar x_i\)를 쓴다.
• 브라운-포사이드 검정: \(c_i\)로 집단 중앙값 \(\tilde x_i\)를 쓴다.

왜 강건한가. \(F\) 검정은 4차 적률(첨도)에 의존하는 통계량을 쓰지만, 르빈류는 \(z\)의 평균을 비교하므로 중심극한정리의 보호를 받는다. 정규성이 깨져도 오류율이 크게 어긋나지 않는다.

비교.

중앙값을 쓰는 브라운-포사이드가 가장 강건해 기본값으로 권장된다. SciPy의 stats.levene(*samples, center='median')이 기본값으로 중앙값을 쓰는 것이 이 때문이다. 바틀렛 검정도 있지만 \(F\) 검정처럼 정규성에 민감하다.

둘 다 쓰되 Q-Q 그림에 더 기대야 한다.

정규성 검정의 구조적 문제.

• 귀무가설이 "정규다"이다. 기각하지 못한다고 해서 정규라는 증거가 아니다. "증거 없음"과 "없다는 증거"는 다르다.
• \(n\)이 작으면 검정력이 없다. \(n=15\)쯤이면 꽤 심하게 치우친 모집단도 잡아내지 못한다. 그런데 정규성이 가장 중요한 때가 바로 \(n\)이 작을 때다.
• \(n\)이 크면 지나치게 예민하다. \(n=5000\)이면 실무적으로 무시할 만한 미세한 어긋남도 \(p<0.001\)로 기각된다. 그런데 \(n\)이 크면 중심극한정리 덕분에 정규성이 덜 중요해진다.

두 문제가 정확히 반대 방향이다. 검정이 답을 주는 영역과 답이 필요한 영역이 어긋나 있다.

Q-Q 그림의 장점. 어긋남의 크기와 방향을 보여 준다. "오른쪽 꼬리가 약간 두껍다", "관측값 두 개만 벗어나 있다", "전체가 휘어 있다"를 구별할 수 있다. 어긋남이 실무적으로 문제가 되는 종류인지 판단할 근거를 준다.

덧붙일 점.

• 분산 추론은 정규성에 극도로 민감하므로 진단이 특히 중요하다. 평균 추론은 중심극한정리가 어느 정도 구해 준다.
• "검정 → 통과 → 카이제곱 방법 사용"이라는 이단계 절차 자체가 전체 오류율을 흐트러뜨린다.
• 정규성이 의심스러우면 검정으로 씨름하기보다 애초에 정규성을 덜 요구하는 방법(부트스트랩, 르빈 검정, 순열검정)으로 옮겨 가는 편이 낫다.

코크런 정리. \(\mathbf{Z} \sim N_n(\mathbf{0}, I)\)이고 이차형식들의 합이

이며 \(\sum_k r_k = n\)이면, 각 \(\mathbf{Z}^\top A_k\mathbf{Z} \sim \chi^2_{r_k}\)이고 서로 독립이다.

핵심은 자유도의 합이 \(n\)과 정확히 같다는 조건이다. 그때만 \(A_k\)들이 서로 직교하는 사영행렬이 되어 독립성이 따라 나온다.

표본분산에의 적용. \(Z_i = (X_i-\mu)/\sigma\)로 두면

이고 \(1 + (n-1) = n\)이다. 코크런 정리에 따라

• \(n\bar Z^2 = \{\sqrt n(\bar X-\mu)/\sigma\}^2 \sim \chi^2_1\),
• \(\sum_i(Z_i-\bar Z)^2 = (n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\),
• 둘은 독립

이다. 앞의 것이 \(\bar X\)만의 함수이고 뒤의 것이 \(S^2\)이므로, \(\bar X\)와 \(S^2\)의 독립성이 따라 나온다. 두 결과가 사실 하나의 정리다.

분산분석에서의 쓰임. \(k\)개 집단, 전체 \(N\)개 관측에서 제곱합이

로 쪼개지고 자유도가

로 정확히 맞아떨어진다. 코크런 정리에 따라 \(H_0\) 아래에서

이고, 따라서 그 비인

이다. \(F\) 검정이 정확히 \(F\) 분포를 따르는 근거가 코크런 정리이며, 분자와 분모의 독립성이 그 핵심이다.

이원배치, 요인설계, 회귀의 제곱합 분해도 모두 같은 논리를 따른다. 자유도가 맞아떨어지는지 확인하는 것이 분산분석표를 검산하는 방법인 이유이기도 하다.

\(t_{0.975,24} = 2.064\)이다.

(가) 신뢰구간 — 모평균 \(\mu\)를 겨냥한다.

(나) 예측구간 — 다음 한 관측값 \(X_{n+1}\)을 겨냥한다. 예측오차의 분산이 \(\sigma^2(1+1/n)\)이므로

(다) 허용구간 — 모집단의 95%를 95% 신뢰수준으로 담는 구간이다. \(\bar x \pm k\,s\) 꼴이며, \(n=25\), 내용비율 0.95, 신뢰수준 0.95에서 \(k \approx 2.63\)이다.

폭의 비교.

셋을 혼동하는 것이 가장 흔한 오류 중 하나다.

• "평균 수명의 95% 신뢰구간이 (98, 102)시간"은 개별 제품의 95%가 그 사이에 있다는 뜻이 전혀 아니다.
• 예측구간과 허용구간은 \(n\to\infty\)에서도 폭이 0으로 가지 않는다. 개별 값의 변동은 자료를 아무리 모아도 사라지지 않기 때문이다.
• 예측구간과 허용구간의 차이는 미묘하다. 예측구간은 "다음 한 개"를, 허용구간은 "모집단의 95%"를 담는다. 후자가 더 강한 주장이라 언제나 더 넓다.

품질관리 규격이나 환경 기준은 대개 허용구간으로 정한다. "제품의 99%가 규격 안에 든다"를 보증해야 하기 때문이다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, B = 25, 100_000
lo, hi = stats.chi2.ppf([0.025, 0.975], n - 1)

print(f"n = {n}, 목표 신뢰수준 95%")
print(f"{'모집단':>12}{'초과첨도':>12}{'실제 포함률':>14}")
cases = [("정규", stats.norm), ("균등", stats.uniform), ("지수", stats.expon),
         ("t(5)", stats.t(5)), ("로그정규", stats.lognorm(1.0))]
for name, d in cases:
    X = d.rvs(size=(B, n), random_state=rng)
    s2, v = X.var(1, ddof=1), d.var()
    cover = np.mean(((n - 1) * s2 / hi <= v) & (v <= (n - 1) * s2 / lo))
    print(f"{name:>12}{float(d.stats(moments='k')):>12.2f}{cover:>14.4f}")

출력:

n = 25, 목표 신뢰수준 95%
         모집단        초과첨도        실제 포함률
          정규        0.00        0.9495
          균등       -1.20        0.9963
          지수        6.00        0.7208
        t(5)        6.00        0.8338
        로그정규      110.94        0.4165

(scipy의 stats(moments='k')는 초과첨도 \(\gamma_2\)를 돌려준다. 이 장에서 쓰는 첨도는 \(\beta_2 = \gamma_2 + 3\)이므로 차례로 \(\beta_2 = 3,\ 1.8,\ 9,\ 9,\ 113.94\)다.)

결과가 참담하다. 로그정규 자료에서 명목 95% 구간이 실제로는 42%만 포함한다. 절반 이상의 경우에 참 분산이 구간 밖에 있다.

본문의 극한값과 견주어 보라. 균등은 극한 0.998에 유한표본 0.9963, 지수는 극한 0.673에 유한표본 0.7208이다. \(n = 25\)에서는 카이제곱 분포 자체가 아직 넓어 어긋남이 조금 가려져 있고, \(n\)을 키우면 극한값 쪽으로 내려간다. 표본을 늘리는 것이 사태를 개선하지 않는다는 뜻이다.

첨도만으로는 다 설명되지 않는다. 지수와 \(t_5\)는 첨도가 \(\beta_2 = 9\)로 똑같은데 포함률이 0.7208과 0.8338로 다르다. 차이를 만드는 것은 치우침이다. 지수는 왜도가 2인 비대칭 분포이고 \(t_5\)는 대칭이다. 배율 \((\beta_2-1)/2\)는 점근적인 폭만 말해 줄 뿐, 유한표본에서 구간이 어느 쪽으로 빗나가는지까지 정하지는 못한다.

대안. 부트스트랩으로 \(S^2\)의 분포를 재표집해 직접 추정하거나, 점근분산에 표본첨도를 넣어 보정하거나, 곱셈적 자료라면 로그를 취한 뒤 분석한다. 분산이라는 척도 자체를 포기하고 사분위범위나 MAD로 옮겨 가는 것도 방법이다. 어느 쪽이든 먼저 첨도를 재는 것이 순서다. \(\hat\beta_2\)가 4를 넘으면 카이제곱 구간을 신뢰하지 않는 편이 안전하다.