모수 대 통계량¶

(1) 해석적으로. 표본평균은 불편이므로 \(E[\bar X] = \mu = 100\)이고, 표준오차는

이다(\(n/N = 10^{-4}\)이라 유한모집단 수정은 무시해도 좋다).

되풀이 \(M = 5{,}000\)번으로 이 두 값을 추정하는 데 따르는 오차는 따로 있다.

• 표본평균들의 평균은 \(M\)개의 독립 관측을 평균한 것이므로 그 표준오차가 \(\operatorname{SE}/\sqrt M = 2.1213/\sqrt{5000} = 0.0300\)이다.
• 표본평균들의 표준편차는 정규표본의 표준편차 추정이므로 그 표준오차가 대략 \(\operatorname{SE}/\sqrt{2(M-1)} = 2.1213/99.99 = 0.0212\)다.

이 두 수가 "얼마나 가까워야 가까운 것인가"의 자다. 이 자가 없으면 \(99.969\)가 \(100\)에 가까운지 아닌지 말할 수 없다.

(2) 수치적으로. 아래 코드를 돌리면 표본평균들의 평균이 \(99.969\), 표준편차가 \(2.097\)로 나온다.

중심. \(99.969\)는 \(100\)에서 \(0.031\) 떨어져 있고 몬테카를로 오차가 \(0.0300\)이므로 \(z = -1.03\)이다. 어긋남이라 할 수 없다. 한 가지 더 정확히 말하자면, 비복원으로 뽑는 대상은 \(N(100, 15^2)\)이 아니라 생성된 \(50\)만 명이고 그들의 평균은 \(100.0258\)이다. 그것을 참값으로 쓰면 \(z = -1.9\)가 되는데, 역시 \(2\) 안이다.

퍼짐. \(2.097\)은 이론값 \(2.1213\)에서 \(0.024\) 떨어져 있고 몬테카를로 오차가 \(0.0212\)이므로 \(z = -1.15\)다. 역시 맞는다.

모집단의 표준편차 \(15\)가 표본평균의 표준편차 \(2.1\)로 줄어든 것이 \(1/\sqrt{50} = 0.1414\)배이며, 추론이 가능한 이유 자체다.

"""표본평균의 표본분포를 눈으로 확인한다."""

import matplotlib
matplotlib.use("Agg")
import numpy as np
import matplotlib.pyplot as plt

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
# 글꼴을 바꾸면 마이너스 기호가 깨지므로 unicode_minus 도 함께 꺼 준다.
plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

rng = np.random.default_rng(0)

# === 모수: 우리가 알아내려는 대상 ===
# 현실에서는 이 두 값을 절대 알 수 없다. 모의실험이므로 답을 미리 정해 둔다.
mu_true = 100
sigma_true = 15
population = rng.normal(mu_true, sigma_true, size=500_000)

# === 표본을 5000번 되풀이해 뽑고, 그때마다 표본평균을 기록한다 ===
# 이것이 이 보기의 핵심이다. 현실에서는 표본을 한 번만 뽑으므로
# 표본평균도 하나뿐이다. 여기서는 "만약 다시 뽑는다면 얼마가 나올까"를
# 5000번 되풀이해 그 분포를 직접 만들어 본다.
n = 50                # 표본 하나의 크기
num_samples = 5_000   # 되풀이 횟수
sample_means = np.array([
    rng.choice(population, n, replace=False).mean()
    for _ in range(num_samples)
])

# === 표본분포에서 두 가지를 확인한다 ===
#   (1) 중심: 표본평균들의 평균이 참 mu와 같은가?  -> 불편성
#   (2) 퍼짐: 그 표준편차가 이론값 sigma/sqrt(n)과 맞는가? -> 표준오차 공식
print(f"참 모평균 mu:      {mu_true}")
print(f"표본평균들의 평균: {sample_means.mean():.3f}")
print(f"이론 표준오차:     {sigma_true / np.sqrt(n):.3f}")
print(f"관측 표준오차:     {sample_means.std(ddof=1):.3f}")

# 5000개의 표본평균을 히스토그램으로 그린다.
# 모집단이 정규분포이므로 표본분포도 정규분포이며, 참 mu를 중심으로 대칭이다.
fig, ax = plt.subplots(figsize=(8, 3))
ax.hist(sample_means, bins=40, density=True, alpha=0.75,
        color="#1565C0", edgecolor="white")
ax.axvline(mu_true, color="#D32F2F", lw=2,
           label=r"참 모평균 $\mu = 100$")      # 참값 표시
ax.set_xlabel("표본평균")
ax.set_ylabel("밀도")
ax.set_title("표본평균의 표본분포 (n = 50)")
ax.legend()
fig.tight_layout()
fig.savefig("sampling_distribution_mean.png", dpi=170, facecolor="white",
            bbox_inches="tight")

출력:

참 모평균 mu:      100
표본평균들의 평균: 99.969
이론 표준오차:     2.121
관측 표준오차:     2.097

숫자와 그림이 같은 말을 한다.

• 중심. 표본평균들의 평균이 99.969로 참값 100과 사실상 같다. 히스토그램의 봉우리도 빨간 선(참 \(\mu\)) 위에 놓여 있다. 표본평균이 \(\mu\)의 불편추정량이라는 뜻이다.
• 퍼짐. 이론값 \(\sigma/\sqrt{n} = 15/\sqrt{50} = 2.121\)과 관측값 2.097이 거의 맞는다.

주목할 것은 모집단의 표준편차가 15인데 표본평균들의 표준편차는 2.1이라는 점이다. 개별 관측값보다 표본평균이 일곱 배쯤 덜 흔들린다. 이 축소가 \(1/\sqrt{n}\) 배이며, 통계적 추론이 가능한 이유 자체다.

(a) 모수 — 성인 모집단 전체를 기술한다. (b) 통계량 — 2,000가구 표본에서 계산했다. (c) 모수 — 생산분 전체를 기술한다. (d) 통계량 — 학생 50명에서 계산했다. (e) 모수 — 동전의 성질이다(개념적 반복의 "모집단"). (f) 통계량 — \(\hat{\beta}_1\)은 모집단 모수 \(\beta_1\)의 추정량이다.

항등식 \(\sum_i (X_i - \bar{X})^2 = \sum_i X_i^2 - n\bar{X}^2\)을 이용한다. 기댓값을 취하면

빼면

따라서 \(\mathbb{E}[S^2] = \frac{(n-1)\sigma^2}{n-1} = \sigma^2\)이다. \(\square\)

대신 \(n\)으로 나누면 \(\mathbb{E}[\tilde S^2] = \frac{n-1}{n}\sigma^2\)가 되어 체계적으로 과소추정하는 추정량이 되며, 특히 \(n\)이 작을 때 해롭다.

\(\mathrm{SE}(\bar X) = \sigma / \sqrt{n} = 10/\sqrt{n}\):

\(n\)이 네 배가 될 때마다 표준오차는 절반이 된다 — 표준적인 \(\sqrt{n}\) 속도다. 정밀도를 크게 높이려면 표본 크기를 자릿수 단위로 늘려야 하며, 이것이 통계적 표본 크기 계획의 기본 경제학이다.

빈도주의적 의미에서는 옳지 않다. 빈도주의 추론에서 \(p\)는 (미지이지만) 고정된 상수이지 확률변수가 아니므로, "\(p\)가 확률 0.95로 \((a, b)\) 안에 있다" 같은 확률 진술은 적용되지 않는다.

95% 신뢰구간의 올바른 해석은 절차에 대한 진술이다. 표집과 신뢰구간 구성을 여러 번 반복하면 그렇게 얻은 구간의 95%가 참값 \(p\)를 포함한다는 뜻이다. 특정 구간 \((0.30, 0.50)\)은 \(p\)를 포함하거나 포함하지 않거나 둘 중 하나이며, \(p\) 자체에 대한 확률 진술은 없다.

베이즈 추론에서는 모수에 대한 확률 진술을 할 수 있지만, 사전분포를 지정하고 사후분포를 보고한 뒤에야 가능하며 이는 다른 방식의 추론이다. 두 해석은 대중 보도에서 자주 혼동된다.

기댓값의 선형성: \(\mathbb{E}[\hat p] = (1/n)\sum \mathbb{E}[X_i] = p\).

독립성: \(\mathrm{Var}(\hat p) = (1/n^2)\sum \mathrm{Var}(X_i) = (1/n^2) \cdot n\, p(1-p) = p(1-p)/n\).

함수 \(p \mapsto p(1-p)\)는 위로 볼록한 포물선으로 \(p = 1/2\)에서 최대가 된다(\(p(1-p) = 0.25\)). 표본 크기를 계획할 때 \(p\)를 모르면 보통 최악의 경우인 \(p(1-p) = 0.25\)를 택하며, 95% 신뢰수준에서 목표 오차한계 ME에 대해 익숙한 공식 \(n \approx 1/(4 \mathrm{ME}^2)\)이 나온다.

정규 i.i.d. 표본에서 \(\sigma^2\)을 추정하는 경우를 생각하자. 최대가능도추정량(MLE)은 \(n - 1\)이 아니라 \(n\)으로 나눈다.

이는 편향되어 있다: \(\mathbb{E}[\tilde S^2] = \frac{n-1}{n}\sigma^2\). 그러나 분산은 불편추정량 \(S^2\)보다 작다. 계산하면

이고,

이다. 유한한 \(n\)에 대해 \(\mathrm{MSE}(\tilde S^2) < \mathrm{MSE}(S^2)\) — MLE는 편향되어 있음에도 MSE가 더 작다.

일반적인 교훈은 이것이다: 편향은 분산과 맞바꿔진다. 축소추정량(제임스–스타인, 능형회귀)은 약간의 편향을 감수하고 분산을 크게 줄임으로써 이 사실을 활용한다. 편향이 언제나 나쁜 것은 아니며, 2차원 정확도 예산에서 조절할 수 있는 하나의 손잡이일 뿐이다.

정규표본에서 \(SS = \sum(X_i-\bar X)^2\)는 \(SS/\sigma^2 \sim \chi^2_{n-1}\)을 만족하므로

이다. 따라서

이고, \(c\)로 미분해 \(0\)으로 두면

를 얻는다. 즉 \(n+1\)로 나누는 것이 MSE를 최소화한다. 이때 \(\operatorname{MSE} = 2\sigma^4/(n+1)\)이다.

import numpy as np

rng = np.random.default_rng(0)
n, B = 8, 400_000
x = rng.normal(0, 1, (B, n))
SS = ((x - x.mean(1, keepdims=True)) ** 2).sum(1)

print(f"{'분모':>6}{'평균':>10}{'편향':>10}{'분산':>10}{'MSE':>10}{'이론 MSE':>10}")
for name, d, mse_th in [("n-1", n - 1, 2 / (n - 1)),
                        ("n",   n,     (2 * n - 1) / n ** 2),
                        ("n+1", n + 1, 2 / (n + 1))]:
    e = SS / d
    print(f"{name:>6}{e.mean():>10.4f}{e.mean() - 1:>+10.4f}"
          f"{e.var():>10.4f}{((e - 1) ** 2).mean():>10.4f}{mse_th:>10.4f}")

출력:

    분모        평균        편향        분산       MSE    이론 MSE
   n-1    0.9996   -0.0004    0.2866    0.2866    0.2857
     n    0.8747   -0.1253    0.2194    0.2351    0.2344
   n+1    0.7775   -0.2225    0.1733    0.2229    0.2222

MSE가 \(0.2866 > 0.2351 > 0.2229\)로 단조 감소하며 이론값과 맞는다. 편향은 반대 방향으로 커진다.

그런데 왜 아무도 \(n+1\)을 쓰지 않는가. 이유가 여러 겹이다.

• 이 결과는 정규성에 의존한다. \(\operatorname{Var}(SS)\) 계산에 정규분포의 4차 적률이 들어갔다. 꼬리가 두꺼운 분포에서는 최적 \(c\)가 달라진다.
• MSE는 여러 기준 중 하나일 뿐이다. 손실함수를 바꾸면 최적값도 바뀐다.
• 불편성은 합쳐질 때 보존된다. 여러 집단의 분산 추정값을 합동(pooled)하거나 분산분석 표를 만들 때, 불편추정량들의 합은 여전히 불편이다. 편향된 추정량은 편향이 누적된다.
• \(S^2\)은 \(\chi^2\) 분포와 깔끔하게 맞물려 \(t\) 검정과 \(F\) 검정이 정확히 성립한다.

교훈은 MSE의 승자를 찾는 것이 아니다. "최적의 추정량"이라는 말이 반드시 어떤 기준에 대해 최적인지를 밝혀야만 뜻을 갖는다는 것이다. 기준을 바꾸면 답도 바뀐다. \(\square\)

import numpy as np
from scipy import stats

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

print(f"{'n':>6}{'통계량':>9}{'평균':>10}{'표준편차':>11}{'왜도':>9}{'첨도':>9}")
for n in (10, 100, 1000):
    x = rng.normal(0, 1, (B, n))
    for name, v in [("Xbar", x.mean(1)), ("s^2", x.var(1, ddof=1)), ("최댓값", x.max(1))]:
        print(f"{n:>6}{name:>9}{v.mean():>10.4f}{v.std():>11.4f}"
              f"{stats.skew(v):>9.3f}{stats.kurtosis(v):>9.3f}")

print("\ns^2 의 이론 왜도 = sqrt(8/(n-1))")
for n in (10, 100, 1000):
    print(f"  n = {n:>5}: {np.sqrt(8 / (n - 1)):.3f}")

출력:

     n      통계량        평균       표준편차       왜도       첨도
    10     Xbar    0.0008     0.3156   -0.006   -0.016
    10      s^2    0.9983     0.4710    0.937    1.292
    10      최댓값    1.5378     0.5874    0.413    0.345
   100     Xbar    0.0001     0.0998   -0.002    0.018
   100      s^2    0.9998     0.1423    0.280    0.108
   100      최댓값    2.5068     0.4298    0.668    0.825
  1000     Xbar    0.0001     0.0317   -0.004    0.010
  1000      s^2    1.0001     0.0447    0.090    0.009
  1000      최댓값    3.2415     0.3520    0.789    1.071

s^2 의 이론 왜도 = sqrt(8/(n-1))
  n =    10: 0.943
  n =   100: 0.284
  n =  1000: 0.089

셋이 전혀 다르게 움직인다.

\(\bar{X}\) — 어떤 \(n\)에서도 왜도와 첨도가 \(0\)이다. 모집단이 정규이므로 \(\bar{X}\)도 정확히 정규다. 표준편차는 \(1/\sqrt{n}\)대로 \(0.316 \to 0.100 \to 0.032\).

\(S^2\) — 오른쪽으로 치우쳐 있다. 왜도가 \(0.937 \to 0.280 \to 0.090\)으로 \(\sqrt{8/(n-1)}\)과 정확히 맞는다. \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)이고 카이제곱분포가 치우쳐 있기 때문이다. 다만 치우침이 \(0\)으로 가므로 \(S^2\)도 결국 정규에 가까워진다. 이것이 중심극한정리의 확장판이 적용되는 경우다.

최댓값 — 완전히 다르다. 왜도가 줄기는커녕 \(0.413 \to 0.668 \to 0.789\)로 커진다. 표준편차도 \(0.587 \to 0.430 \to 0.352\)로 느릿느릿 줄 뿐이다(\(\bar{X}\)처럼 \(1/\sqrt{n}\)이 아니라 대략 \(1/\log n\)이다).

왜 최댓값에는 중심극한정리가 통하지 않는가. 중심극한정리는 많은 항의 합에 관한 정리다. 각 관측이 조금씩 기여해 개별 기여가 씻겨 나가면서 정규성이 나온다. 최댓값은 합이 아니라 하나의 관측이며, 나머지 \(n-1\)개는 그 값을 고르는 데만 쓰인다. 그래서 극한분포가 정규가 아니라 검벨분포(왜도 \(1.1395\))이며, 위 수치가 그쪽으로 올라가고 있다.

이것이 극단값 이론의 출발점이다. 홍수 수위, 보험의 대형 손실, 금융의 최대 낙폭처럼 최댓값이 관심사인 문제에서는 정규분포에 기댄 도구를 쓰면 안 된다. 검벨·프레셰·바이불의 세 극한분포가 그 자리를 대신한다.

(a) 불편이지만 일치성 없음: 첫 관측값만 쓰기. \(\hat\mu = X_1\)로 두면 \(\mathbb{E}[X_1] = \mu\)이므로 불편이다. 그러나 \(\operatorname{Var}(X_1) = \sigma^2\)이 \(n\)과 무관하므로 자료를 아무리 모아도 정밀해지지 않는다.

(b) 편향이지만 일치성 있음: 최대가능도 분산. \(\tilde S^2 = \frac{1}{n}\sum(X_i-\bar X)^2\)는 편향이 \(-\sigma^2/n\)이지만, 이 편향이 \(n \to \infty\)에서 \(0\)으로 가고 분산도 \(0\)으로 가므로 일치성을 갖는다.

import numpy as np

rng = np.random.default_rng(0)

print("(a) X_1 은 불편이지만 일치성이 없다")
print(f"{'n':>7}{'E[X_1]':>10}{'SD[X_1]':>10}{'E[Xbar]':>10}{'SD[Xbar]':>10}")
for n in (10, 100, 1000, 10_000):
    x = rng.normal(5, 2, (20_000, n))
    print(f"{n:>7}{x[:, 0].mean():>10.4f}{x[:, 0].std():>10.4f}"
          f"{x.mean(1).mean():>10.4f}{x.mean(1).std():>10.4f}")

print("\n(b) MLE 분산은 편향되어 있지만 일치성이 있다")
print(f"{'n':>7}{'평균':>10}{'편향':>10}{'이론 편향':>11}{'표준편차':>11}")
for n in (5, 20, 100, 1000):
    v = rng.normal(0, 1, (50_000, n)).var(1, ddof=0)
    print(f"{n:>7}{v.mean():>10.5f}{v.mean() - 1:>+10.5f}"
          f"{-1 / n:>+11.5f}{v.std():>11.4f}")

출력:

(a) X_1 은 불편이지만 일치성이 없다
      n    E[X_1]   SD[X_1]   E[Xbar]  SD[Xbar]
     10    4.9880    2.0113    5.0003    0.6308
    100    5.0252    1.9937    5.0011    0.2005
   1000    5.0100    2.0053    5.0001    0.0631
  10000    4.9877    2.0242    5.0000    0.0202

(b) MLE 분산은 편향되어 있지만 일치성이 있다
      n        평균        편향      이론 편향       표준편차
      5   0.80233  -0.19767   -0.20000     0.5714
     20   0.94886  -0.05114   -0.05000     0.3067
    100   0.98978  -0.01022   -0.01000     0.1403
   1000   0.99885  -0.00115   -0.00100     0.0449

(a)에서 E[X_1]은 언제나 \(5\) 근처지만 SD[X_1]은 \(2\)에 붙박여 있다. 반면 SD[Xbar]는 \(0.632 \to 0.020\)으로 떨어진다. \(X_1\)은 표적을 평균적으로 맞히지만 결코 표적에 수렴하지 않는다.

(b)에서 편향이 \(-0.197 \to -0.0012\)로 이론값 \(-1/n\)을 따라 사라지고 표준편차도 함께 줄어든다.

두 개념의 관계.

어느 쪽이 더 중요한가. 거의 언제나 일치성이다. 일치성이 없다는 것은 자료를 무한히 모아도 답을 알 수 없다는 뜻이라 치명적이다. 불편성은 "평균적으로 맞다"일 뿐인데, 우리는 보통 자료를 한 번 얻으므로 여러 표본에 걸친 평균이 맞다는 것만으로는 위안이 되지 않는다.

가장 널리 쓰이는 추정량 중 하나인 최대가능도추정량은 일반적으로 편향되어 있고, 일치성과 점근정규성 때문에 쓰인다. \(\square\)

\(g(x) = 1/x\)는 볼록함수이므로 옌센 부등식에 의해 \(\mathbb{E}[1/\bar{X}] > 1/\mathbb{E}[\bar{X}] = 1/\mu\)이다. 즉 언제나 과대추정한다.

크기를 알려면 \(\mu\) 주위에서 2차까지 전개한다.

기댓값을 취하면 1차항이 사라지고 \(\mathbb{E}[(\bar{X}-\mu)^2] = \sigma^2/n\)이므로

이다. \(g(x)=1/x\)이면 \(g''(x) = 2/x^3\)이라

를 얻는다.

import numpy as np

rng = np.random.default_rng(0)
mu, sigma = 5., 2.

print(f"{'n':>6}{'E[1/Xbar]':>12}{'1/mu':>10}{'델타 근사':>12}{'실제 편향':>12}")
for n in (5, 20, 100, 1000):
    xbar = rng.normal(mu, sigma, (400_000, n)).mean(1)
    emp = np.mean(1 / xbar)
    approx = 1 / mu + sigma ** 2 / (n * mu ** 3)
    print(f"{n:>6}{emp:>12.6f}{1 / mu:>10.6f}{approx:>12.6f}{emp - 1 / mu:>+12.6f}")

출력:

     n   E[1/Xbar]      1/mu       델타 근사       실제 편향
     5    0.207046  0.200000    0.206400   +0.007046
    20    0.201690  0.200000    0.201600   +0.001690
   100    0.200307  0.200000    0.200320   +0.000307
  1000    0.200029  0.200000    0.200032   +0.000029

편향이 언제나 양수이고, \(n \ge 20\)부터는 델타 근사가 소수점 다섯째 자리까지 맞는다. 편향이 \(1/n\)의 속도로 줄어드는 것도 확인된다.

일반 원리. \(\hat\theta\)이 \(\theta\)의 불편추정량이어도 \(g(\hat\theta)\)은 \(g\)가 비선형인 한 \(g(\theta)\)의 불편추정량이 아니다. 편향의 부호는 \(g\)의 볼록성이 정한다. \(g\)가 볼록이면 위로, 오목이면 아래로 치우친다(앞 절 모집단과 표본의 연습문제 10에 나온 \(\mathbb{E}[S] < \sigma\)가 오목한 경우다).

어디서 문제가 되는가.

• 비율과 역수: 평균 대기시간에서 처리율(\(1/\bar{T}\))을 구하거나, 평균 비용에서 효율을 구할 때
• 오즈비와 상대위험도: \(\log\)를 취해 다루는 것이 표준인 이유가 대칭성뿐 아니라 편향 때문이기도 하다
• 로그정규 자료의 역변환: \(\exp(\overline{\log X})\)는 산술평균이 아니라 기하평균을 추정한다. 로그척도에서 회귀한 뒤 지수를 취하면 평균을 과소추정한다

처방. (1) \(n\)이 크면 편향이 \(1/n\)이라 대개 무시할 만하다. (2) 델타 방법으로 크기를 가늠해 표준오차와 비교하라. 편향이 표준오차보다 훨씬 작으면 신경 쓰지 않아도 된다. (3) 그렇지 않으면 부트스트랩으로 편향을 추정해 빼거나, 애초에 관심 있는 척도에서 직접 모형을 세워라. \(\square\)