평균의 표본분포¶

크기 25인 표본을 반복해서 뽑으면 표본평균들이 100 주위에 모이며 전형적인 편차는 0.8이다.

from scipy import stats
print(f"P(X_bar > 155) = {stats.norm.sf(1.25):.4f}")

출력:

P(X_bar > 155) = 0.1056

from scipy import stats
print(f"P(6.8 < X_bar < 7.2) = {stats.norm.cdf(0.93) - stats.norm.cdf(-0.93):.4f}")

출력:

P(6.8 < X_bar < 7.2) = 0.6476

모집단이 정규이므로 \(n = 5\)에서도 결과가 정확하다:

from scipy import stats
print(f"P(X_bar > 72) = {stats.norm.sf(0.447):.4f}")

출력:

P(X_bar > 72) = 0.3274

물이 떨어진다는 것은 \(\bar{X} > 110/50 = 2.2\)를 뜻한다:

\(n = 5\)이고 정규성 가정이 없으면 중심극한정리를 믿고 적용할 수 없다. 모집단 모양에 관한 추가 정보 없이는 이 확률을 구할 수 없다.

모집단이 치우쳐 있지만 \(n = 100\)이면 중심극한정리를 쓰기에 충분히 크다:

from scipy import stats
print(f"P(X_bar > 2100) = {stats.norm.sf(2):.4f}")

출력:

P(X_bar > 2100) = 0.0228

(1) 이론값. \(\text{Uniform}(0,1)\)은 \(\mu = 1/2\), \(\sigma = 1/\sqrt{12} = 0.288675\)이므로

이다. 이 두 값에는 근사가 없다. 모집단이 균등이든 무엇이든 정확하다.

여기서 한 걸음 더 나가는 것이 요령이다. 모의실험이 내놓는 두 숫자 자체가 확률변수이므로 그것들의 흔들림도 미리 계산할 수 있다. 되풀이 횟수를 \(B = 10{,}000\)이라 하면

이고, 표준편차 추정값의 흔들림은 대략

이다. 그러므로 \(0.5000\)과 \(0.1291\)이 그대로 나올 것을 기대하면 안 된다. 소수 셋째 자리에서 흔들리는 것이 정상이며, 어긋남이 위 두 수의 두세 배 안에 들면 공식이 맞은 것이다.

점선이 가두는 확률도 계산해 둔다. \(n = 5\)에서 \(\bar X\)의 정확한 분포는 어윈-홀 분포를 \(5\)로 나눈 것이므로

를 얻는다. 정규분포라면 \(2\Phi(1) - 1 = 0.6827\)이니 \(0.0103\) 작다. 균등의 \(\bar X\)는 꼬리가 정규보다 짧아(\(n = 5\)에서 초과첨도 \(-0.24\)) 가운데 한 표준편차 안에 든 질량이 오히려 적다.

(2) 모의실험.

import matplotlib.pyplot as plt
import numpy as np
from math import comb, factorial
from scipy import stats

np.random.seed(0)

# 크기 5짜리 균등표본 U(0,1)을 1만 번 뽑아 그때마다 표본평균을 기록한다.
X_bar = []
for _ in range(10_000):
    x = np.random.uniform(size=(5,))
    X_bar.append(x.mean())

# 1만 개 표본평균의 평균과 표준편차.
#   평균     -> 참 평균 0.5 에 가까워야 한다
#   표준편차 -> 이것이 표준오차다. 이론값은 sigma/sqrt(n) 이며
#              U(0,1)의 sigma = 1/sqrt(12) 이므로
#              (1/sqrt(12))/sqrt(5) = 0.1291 이 나와야 한다.
average = np.array(X_bar).mean()
standard_error = np.array(X_bar).std()

print(f"Estimated Mean of X_bar:  {average:.4f}")
print(f"Standard Error of X_bar:  {standard_error:.4f}")

# 이론값과 1만 번 되풀이에서 기대할 요동을 함께 적는다.
n, B = 5, 10_000
mu_true, se_true = 0.5, 1 / np.sqrt(12 * n)
print(f"이론값  mu = {mu_true:.4f},  sigma/sqrt(n) = {se_true:.4f}")
print(f"몬테카를로 오차  평균 {se_true / np.sqrt(B):.4f},  표준오차 {se_true / np.sqrt(2 * B):.4f}")
print(f"어긋남을 오차 단위로  평균 {abs(average - mu_true) / (se_true / np.sqrt(B)):.2f}배,  "
      f"표준오차 {abs(standard_error - se_true) / (se_true / np.sqrt(2 * B)):.2f}배")

# 점선 두 개가 가두는 확률. n=5 의 정확한 분포는 어윈-홀을 5로 나눈 것이다.
def cdf_irwin_hall(x, n):
    x = np.asarray(x, dtype=float)
    total = np.zeros_like(x)
    for k in range(n + 1):
        total += (-1) ** k * comb(n, k) * np.clip(x - k, 0, None) ** n
    return np.clip(total / factorial(n), 0, 1)

p_exact = float(cdf_irwin_hall(n * (mu_true + se_true), n) - cdf_irwin_hall(n * (mu_true - se_true), n))
p_normal = float(stats.norm.cdf(1) - stats.norm.cdf(-1))
p_sim = float(np.mean(np.abs(np.array(X_bar) - average) <= standard_error))
print(f"점선 사이의 확률  정확 {p_exact:.4f},  정규근사 {p_normal:.4f},  모의 {p_sim:.4f}")

# 히스토그램에 이론값을 겹쳐 그린다.
fig, ax = plt.subplots(figsize=(12, 3))
ax.set_title("Sampling Distribution of X-bar")
ax.hist(X_bar, bins=100, density=True, alpha=0.3)
ax.vlines(average, ymin=0, ymax=5, color="k", lw=5, label="Mean")
ax.vlines(average + standard_error, ymin=0, ymax=5,
          color="k", ls="--", label="Mean +/- SE")
ax.vlines(average - standard_error, ymin=0, ymax=5,
          color="k", ls="--")
ax.legend()
plt.show()

출력:

Estimated Mean of X_bar:  0.4981
Standard Error of X_bar:  0.1287
이론값  mu = 0.5000,  sigma/sqrt(n) = 0.1291
몬테카를로 오차  평균 0.0013,  표준오차 0.0009
어긋남을 오차 단위로  평균 1.47배,  표준오차 0.42배
점선 사이의 확률  정확 0.6724,  정규근사 0.6827,  모의 0.6685

두 숫자 모두 예측한 띠 안에 든다. 모의 평균 \(0.4981\)이 \(0.5\)에서 \(0.0019\) 떨어져 있어 몬테카를로 오차 \(0.0013\)의 \(1.47\)배, 모의 표준오차 \(0.1287\)이 \(0.1291\)에서 \(0.0004\) 떨어져 있어 \(0.0009\)의 \(0.42\)배다. "가깝다"가 아니라 "오차의 \(1.5\)배와 \(0.4\)배"라고 말할 수 있다는 것이 이 보기의 쓸모다. 둘 다 \(2\)배 안쪽이므로 공식과 어긋난다고 볼 근거가 없다.

띠의 폭을 알면 반대 방향으로도 쓸 수 있다. 모의 표준오차를 소수 넷째 자리까지 믿고 싶다면 \(\mathrm{SE}/\sqrt{2B} < 0.00005\), 곧 \(B > 3{,}300{,}000\)번을 되풀이해야 한다. 1만 번으로 셋째 자리까지가 한계다.

점선 두 개. 모의 비율 \(0.6685\)가 정확값 \(0.6724\)에서 \(0.0039\) 떨어져 있고, 비율의 몬테카를로 오차가 \(\sqrt{0.672 \times 0.328/10^4} = 0.0047\)이니 \(0.8\)배다. 반면 정규근사 \(0.6827\)과는 \(0.0142\), 곧 오차의 \(3.0\)배나 떨어져 있다. 1만 번의 모의실험이 "정확한 어윈-홀"과 "정규근사"를 구별해 낼 만큼 촘촘하다는 뜻이고, 그래서 \(n = 5\)에서 \(\bar X\)가 아직 정규가 아니라는 것을 모의실험만으로도 보일 수 있다.

다만 공식 \(\sigma/\sqrt n\) 쪽에는 그런 어긋남이 없다. 중심과 폭은 모집단의 모양과 무관하게 정확하고, 모양만 근사다. 이 쪽이 처음부터 말해 온 구분이 숫자로 갈라져 나온 자리다.

(a) \(\mathbb{E}[\bar X] = \mu = 75\). 표본평균은 \(n\)과 무관하게 불편이다.

(b) \(\mathrm{SE}(\bar X) = \sigma/\sqrt n = 18/3 = 6\). \(n\)을 네 배인 36으로 늘리면 표준오차가 절반인 3이 된다.

\(\mathrm{SE} = 20/\sqrt{25} = 4\). \(Z = (155 - 150)/4 = 1.25\).

\(P(\bar X > 155) = 1 - \Phi(1.25) = 0.1056\)으로 약 10.6%이다.

참고: 모집단의 정규성이 주어졌으므로 \(\bar X\)는 어떤 \(n\)에 대해서도 정확히 정규분포이다. 중심극한정리가 필요 없다.

\(\mathrm{SE} = 100/\sqrt 5 \approx 44.7\). 그대로 계산하면 \(Z = 10/44.7 \approx 0.224\)이므로 \(P \approx 0.41\)이다.

그러나: \(n = 5\)는 중심극한정리를 쓰기에 너무 작고 모집단의 모양도 알려져 있지 않다. \(\bar X\)의 표본분포는 심하게 치우쳤거나(예: 수명이 지수분포일 때) 꼬리가 두꺼울 수 있다. 정규근사가 크게 빗나갈 수 있다.

Chebyshev 부등식을 이용한 보수적 한계: \(P(|\bar X - 800| \ge 10) \le \sigma^2/(n \cdot 10^2) = 10000/(5 \cdot 100) = 20\)으로 아무 정보도 주지 못한다.

실용적 결론: 분포 가정이 없거나 \(n\)이 더 크지 않으면 점 추정값을 포기하고 보수적 한계만 보고해야 한다. 실제 전구 수명은 흔히 와이불분포를 따르며, 이를 가정하면 정확한 계산이 가능하다.

\(\mathrm{SE} = 500/10 = 50\). \(Z = (2100 - 2000)/50 = 2\).

\(P(\bar X > 2100) = 1 - \Phi(2) = 0.0228\)로 약 2.3%이다.

왜 타당한가: \(n = 100\)은 크므로 중심극한정리에 의해 모집단 모양과 무관하게 \(\bar X\)가 근사적으로 정규분포이다. 모집단의 왜도 \(\gamma_1\)은 \(\bar X\)에서 \(\gamma_1/\sqrt{100} = \gamma_1/10\)로 줄어들므로, 치우침이 온건한 수준이라면 남는 몫이 무시할 만하다.

주의: 모집단의 치우침이 심하면(예: \(\sigma_{\log} > 1\)인 로그정규) \(n = 100\)에서도 \(\bar X\)에 치우침이 남을 수 있다. Berry-Esseen 한계에 따르면 정규근사의 오차는 \(\sim \mathbb{E}|X|^3/(\sigma^3 \sqrt n)\)이다. 약간 치우친 매출 자료라면 \(n = 100\)으로 충분하지만, 오른쪽으로 심하게 치우친 자료에서는 더 큰 \(n\)이나 붓스트랩이 안전하다.

다음이 필요하다:

최소 표본크기 \(n = 256\)이 필요하다. \(\square\)

\(n > 0\)에서 \(f(n) = \sigma / \sqrt{n} = \sigma \cdot n^{-1/2}\)이라 하자.

1계도함수:

따라서 \(f\)는 순감소한다. 표본이 클수록 언제나 표준오차가 작아진다.

2계도함수:

따라서 \(f\)는 순볼록이다. 볼록성은 \(n\)이 커질수록 표준오차의 감소 속도가 느려짐을 뜻한다. 실용적으로 말하면, 관측값을 하나 더 얻을 때마다 표준오차가 줄어드는 폭이 직전보다 작아진다. 표본크기를 늘리는 데에는 한계 수확 체감이 있다. \(\square\)

모든 쌍의 상관이 \(\rho\)로 같다고 하자(등상관 구조). 그러면

이므로

이다. \(\rho = 0\)이면 익숙한 \(\sigma^2/n\)으로 돌아온다.

\(\rho > 0\)이고 \(n\to\infty\)이면

이다. 표본을 아무리 늘려도 표준오차가 0으로 가지 않는다. 관측값들이 같은 정보를 되풀이해 담고 있어 새 관측이 새 정보를 주지 못하기 때문이다.

괄호 안의 \(1+(n-1)\rho\)를 설계효과(DEFF)라 하고, \(n/\text{DEFF}\)를 유효 표본크기라 한다. \(\rho=0.05\)이고 \(n=20\)인 군집이 50개면 관측이 1000개지만 유효 표본크기는

으로 절반에 지나지 않는다.

현실에서 이런 상황이 흔하다. 같은 학교 학생들, 같은 환자에게서 반복 측정한 값, 시계열의 이웃한 관측값이 모두 그렇다. 독립을 가정하고 \(\sigma/\sqrt n\)을 쓰면 표준오차를 심하게 과소평가하고, 그 결과 신뢰구간이 좁아지고 있지도 않은 유의성이 나온다. 군집 표준오차, 혼합효과 모형, 일반화추정방정식이 모두 이 문제를 다루는 도구다.

부트스트랩을 쓴다.

rng = np.random.default_rng(0)
x = rng.uniform(0, 1, 5)              # 우리가 가진 표본은 이것뿐
boot = [rng.choice(x, size=5, replace=True).mean() for _ in range(10_000)]
print(f"부트스트랩 SE = {np.std(boot, ddof=1):.4f}")

출력:

부트스트랩 SE = 0.1421

원리. 참 모집단에서 반복 표집하는 대신 관측된 표본 자체를 모집단으로 삼아 복원추출한다. 표본이 모집단을 잘 대신한다면 재표집으로 얻은 평균들의 흩어짐이 참 표집분포의 흩어짐에 가깝다.

보기 7과의 관계. 보기 7의 모의실험은 참 모집단 \(\text{Uniform}(0,1)\)에서 1만 번 표집했다. 그것이 가능한 것은 모집단을 알기 때문이며, 현실에서는 불가능하다. 부트스트랩은 같은 일을 표본 하나만 가지고 흉내 낸다.

\(n=5\)에서는 결과가 그리 좋지 않다. 얻은 값 0.1421을 이론값 \(\sigma/\sqrt5 = 0.1291\)과 곧바로 견주면 안 된다. 부트스트랩이 흉내 내는 것은 참 모집단이 아니라 손에 쥔 표본이므로, 비교 대상은 그 표본의 \(s/\sqrt n\)이다. 이 표본은 \(s = 0.357\)이라 \(s/\sqrt5 = 0.1597\)이었다.

그 값과 견주면 부트스트랩 SE가 체계적으로 작다. 복원추출이 \(s^2\)이 아니라 \(\frac{n-1}{n}s^2\)을 재현하기 때문이며, 비율이 \(\sqrt{4/5} = 0.894\)다. 실제로 \(0.894 \times 0.1597 = 0.1428\)로 관측값 0.1421과 맞는다.

이론값 0.1291에서 벗어난 몫의 대부분은 이 편향이 아니라 표본 하나의 운이다. \(n=5\)짜리 표본에서 \(s\) 자체가 크게 흔들리기 때문이며, 다른 시드를 쓰면 이론값보다 작게 나오기도 한다. 두 종류의 오차를 구별해 두어야 한다. 편향은 \(n\)이 크면 사라지고, 표본의 운은 \(n\)이 커야 잦아든다.

부트스트랩을 쓸 이유. 평균이라면 \(s/\sqrt n\)이라는 공식이 있으니 굳이 필요 없다. 부트스트랩이 빛나는 것은 공식이 없는 통계량이다. 중앙값, 절사평균, 사분위수 범위, 두 추정값의 비, 상관계수, 회귀에서 파생된 복잡한 양의 표준오차를 모두 같은 절차로 얻을 수 있다.

표본중앙값 \(\tilde X\)의 점근분포는

이다(\(m\)은 모중앙값, \(f\)는 밀도). 정규분포에서는 \(m=\mu\)이고

이므로

이다. \(\square\)

비교. \(\operatorname{SE}(\bar X) = \sigma/\sqrt n\)이므로

로, 중앙값의 분산이 평균의 1.57배다. 정규모집단에서 중앙값을 쓰면 표본의 36%를 버리는 셈이다. 표본 100개로 중앙값을 내는 것이 표본 64개로 평균을 내는 것과 같다.

그래도 중앙값을 쓸 때.

꼬리가 두꺼우면 순서가 뒤집힌다는 점이 중요하다. \(t_3\) 모집단에서는 중앙값의 점근효율이 평균의 1.6배이고, 코시에서는 평균이 아예 수렴하지 않는 반면 중앙값은 제대로 작동한다.

타협안으로 절사평균이 있다. 양끝 10~20%를 버리고 평균을 내면 정규분포에서 효율을 거의 잃지 않으면서(95% 이상) 이상치에 대한 저항성을 얻는다. 후버 M-추정량도 같은 절충이다.

\(\hat\mu = \sum_i a_iX_i\)라 하면 불편성은 \(E[\hat\mu] = \mu\sum_i a_i = \mu\), 즉 \(\sum_i a_i = 1\)과 같다. 독립이므로

이다. 따라서 \(\sum a_i = 1\) 아래에서 \(\sum a_i^2\)을 최소로 하면 된다.

코시-슈바르츠로. \((\sum_i a_i\cdot1)^2 \le (\sum_i a_i^2)(\sum_i 1^2)\)이므로

이고 등호는 \(a_i\)가 모두 같을 때, 즉 \(a_i = 1/n\)일 때만 성립한다. 그때 \(\operatorname{Var} = \sigma^2/n\)이다. \(\square\)

다른 증명. \(a_i = 1/n + d_i\)로 두면 \(\sum d_i = 0\)이고

으로 곧바로 나온다.

정규성은 필요 없다. 쓴 것은 \(E[X_i]=\mu\), \(\operatorname{Var}(X_i)=\sigma^2\), 그리고 무상관성(공분산 0)뿐이다. 독립일 필요조차 없고 무상관이면 충분하다. 이 점에서 가우스-마르코프 정리와 같은 성격이며, 실제로 이 결과는 가우스-마르코프의 가장 단순한 경우다.

범위의 한계. "선형"과 "불편"이라는 두 제약을 걸었다는 점이 중요하다.

• 제약을 풀면 더 나은 추정량이 있을 수 있다. 꼬리가 두꺼운 모집단에서 절사평균이나 중앙값이 \(\bar X\)보다 분산이 작다(비선형이라 이 정리를 비껴간다).
• 불편성을 포기하면 평균제곱오차가 더 작은 추정량이 있다. 3차원 이상에서 여러 평균을 동시에 추정할 때 제임스-스타인 추정량이 \(\bar X\)를 일률적으로 이긴다는 사실이 유명하다.
• 모집단이 정규라면 \(\bar X\)는 선형이 아닌 추정량까지 통틀어 최소분산불편추정량이다(완비충분통계량의 함수이므로). 정규성은 결론을 "선형" 제약 없이 확장하는 데 쓰인다.