X̄의 표본분포 (Exponential)¶

(1) 이론값. \(\text{Exp}(1)\)은 \(\mu = 1\), \(\sigma = 1\)이므로 중심과 폭은 치우침과 무관하게 곧바로 나온다.

모양도 정확히 안다. \(\bar X \sim \text{Gamma}(\text{형상}=5,\ \text{비율}=5)\)이고 \(\text{Gamma}(k)\)의 왜도는 \(2/\sqrt{k}\), 초과첨도는 \(6/k\)이므로

이다. 일반 공식 \(\gamma_1/\sqrt{n}\), \((\beta_2-3)/n\)에 \(\text{Exp}(1)\)의 \(\gamma_1 = 2\), \(\beta_2 = 9\)를 넣은 것과 같은 값이다. 왜도가 \(0.894\)라는 것은 \(n = 5\)에서 치우침이 거의 그대로 남아 있다는 뜻이다. 균등 모집단에서는 같은 \(n\)에서 왜도가 처음부터 \(0\)이었다.

꼬리확률의 오차는 감마와 정규를 각각 계산해 견주면 된다. \(z = (2-1)/0.44721 = 2.23607\)이므로

정규근사가 꼬리확률을 \(2.3\)배 작게 본다. 중심과 폭은 소수점까지 맞게 주면서도 꼬리는 절반 이하로 틀리는 것이며, \(\sigma\)를 아는데도 그렇다. 덧붙여 정규근사는 \(P(\bar X \le 0) = \Phi(-2.23607) = 0.012674\)라는 값을 내놓는데 \(\bar X > 0\)은 확률 \(1\)이다. 두 수가 똑같은 것은 \(0\)과 \(2\)가 평균 \(1\)에서 같은 거리에 있기 때문이고, 정규근사가 없는 대칭을 강제한 흔적이 바로 그것이다.

(2) 모의실험.

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

np.random.seed(1)

sample_size = 5
n_samples = 10_000
n_population = 10_000

# Exp(1)에서 큰 모집단을 만든다. 치우친 모집단이다.
population = np.random.exponential(size=(n_population,))

# 표본을 딱 하나 뽑는다. 현실에서 우리가 실제로 갖게 되는 것이 이것뿐이다.
# 아래 가운데 패널에 점 몇 개로 그려진다.
single_sample = np.random.choice(population, size=sample_size, replace=False)

# 표본을 되풀이해 뽑으며 표본평균을 기록한다. 이 값들의 분포가 표집분포다.
sample_means = [
    np.mean(np.random.choice(population, size=sample_size, replace=False))
    for _ in range(n_samples)
]
sample_means = np.asarray(sample_means)

# 중심과 폭: Exp(1) 이론값, 실현 모집단, 모의값을 나란히 둔다.
mu_pop, sd_pop = population.mean(), population.std()
print(f"Exp(1) 이론        평균 = 1.0000,  sigma/sqrt(n) = {1 / np.sqrt(sample_size):.4f}")
print(f"실현 모집단 1만 개  평균 = {mu_pop:.4f},  s/sqrt(n)     = {sd_pop / np.sqrt(sample_size):.4f}")
print(f"모의 1만 회        평균 = {sample_means.mean():.4f},  표준편차      = {sample_means.std(ddof=1):.4f}")

# 모양. X-bar ~ Gamma(n, rate=n) 이므로 Exp(1) 이론은 왜도 2/sqrt(n), 초과첨도 6/n 이다.
# 그런데 실현된 1만 개는 Exp(1) 의 꼬리를 다 담지 못하므로 겨냥값이 따로 있다.
print(f"\nExp(1) 이론        왜도 = {2 / np.sqrt(sample_size):+.4f},  초과첨도 = {6 / sample_size:+.4f}")
print(f"실현 모집단 1만 개  왜도 = {stats.skew(population):+.4f},  초과첨도 = {stats.kurtosis(population):+.4f}  (최대값 {population.max():.2f})")
print(f"  거기서 나오는 X-bar 겨냥값   = {stats.skew(population) / np.sqrt(sample_size):+.4f},  {stats.kurtosis(population) / sample_size:+.4f}")
print(f"모의 1만 회        왜도 = {stats.skew(sample_means):+.4f},  초과첨도 = {stats.kurtosis(sample_means):+.4f}")

# 모양 통계량의 몬테카를로 요동은 정규 공식 sqrt(6/B), sqrt(24/B) 를 쓸 수 없다.
# 그 공식은 정규모집단에서만 맞다. 그래서 실제로 400 번 되풀이해 재어 본다.
rng = np.random.default_rng(11)
rep = np.array([
    [stats.skew(m), stats.kurtosis(m), m.std(ddof=1)]
    for m in (population[rng.integers(0, n_population, size=(n_samples, sample_size))].mean(axis=1)
              for _ in range(400))
])
print(f"  같은 일을 400 번 되풀이한 요동 = {rep[:, 0].std():.4f},  {rep[:, 1].std():.4f}  (표준편차는 {rep[:, 2].std():.4f})")

# 정확분포(감마)와 정규근사를 견준다.
exact = stats.gamma(a=sample_size, scale=1 / sample_size)
approx = stats.norm(loc=1, scale=1 / np.sqrt(sample_size))
grid = np.linspace(0, 5, 200_001)
gap = np.abs(exact.cdf(grid) - approx.cdf(grid))
print(f"\n정확분포 대 정규근사  최대 누적분포 차 = {gap.max():.5f}  (x = {grid[gap.argmax()]:.4f})")
for x in (1.5, 2.0):
    print(f"P(X-bar > {x})  정확 = {exact.sf(x):.6f},  정규근사 = {approx.sf(x):.6f},  모의 = {(sample_means > x).mean():.6f}")
print(f"P(X-bar <= 0)  정확 = 0.000000,  정규근사 = {approx.cdf(0):.6f}")

# 모집단과 표집분포를 나란히 그린다.
# 세 패널을 sharex=True 로 묶는 것이 이 그림의 핵심 장치다.
# 가로 눈금이 같아야 세 분포의 **퍼짐**을 직접 견줄 수 있다.
#   위   모집단      : 가장 넓다
#   가운데 표본 하나  : 모집단에서 뽑은 점 몇 개
#   아래  표집분포    : 눈에 띄게 좁다. 이 좁아짐이 sigma/sqrt(n) 이다.
fig, (ax0, ax1, ax2) = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

_, bins, _ = ax0.hist(population, bins=100)
ax0.set_title("Population Distribution (Exponential)")

ax1.scatter(single_sample, np.zeros_like(single_sample), s=100)
ax1.set_title(f"Sample Distribution (n = {sample_size})")

ax2.hist(sample_means, bins=bins)
ax2.set_title("Sampling Distribution of X-bar")

plt.tight_layout()
plt.show()

출력:

Exp(1) 이론        평균 = 1.0000,  sigma/sqrt(n) = 0.4472
실현 모집단 1만 개  평균 = 0.9885,  s/sqrt(n)     = 0.4387
모의 1만 회        평균 = 0.9886,  표준편차      = 0.4312

Exp(1) 이론        왜도 = +0.8944,  초과첨도 = +1.2000
실현 모집단 1만 개  왜도 = +1.9196,  초과첨도 = +5.1090  (최대값 8.99)
  거기서 나오는 X-bar 겨냥값   = +0.8585,  +1.0218
모의 1만 회        왜도 = +0.8353,  초과첨도 = +0.9186
  같은 일을 400 번 되풀이한 요동 = 0.0340,  0.1625  (표준편차는 0.0036)

정확분포 대 정규근사  최대 누적분포 차 = 0.05963  (x = 0.9835)
P(X-bar > 1.5)  정확 = 0.132062,  정규근사 = 0.131776,  모의 = 0.121500
P(X-bar > 2.0)  정확 = 0.029253,  정규근사 = 0.012674,  모의 = 0.023900
P(X-bar <= 0)  정확 = 0.000000,  정규근사 = 0.012674

중심과 폭은 맞는다. 비복원으로 뽑으므로 겨냥값은 실현된 1만 개의 평균 \(0.9885\)이고, 모의 평균 \(0.9886\)이 몬테카를로 오차 \(0.0044\)의 \(0.02\)배만큼 떨어져 있다. 모의 표준편차 \(0.4312\)는 겨냥값 \(0.4387\)보다 \(1.7\%\) 작은데, 400 번 되풀이해 재어 본 요동이 \(0.0036\)이므로 2배쯤 어긋난 셈이다.

모양에서 재미있는 일이 생긴다. 모의 왜도 \(0.8353\)은 \(\text{Exp}(1)\) 이론값 \(0.8944\)보다 작고, 초과첨도 \(0.9186\)은 이론값 \(1.2\)보다 한참 작다. 요동이 각각 \(0.0340\)과 \(0.1625\)이니 \(1.7\)배씩 어긋난 것으로, 못 봐줄 정도는 아니지만 둘이 같은 방향으로 치우쳐 있는 것이 수상하다.

범인은 몬테카를로가 아니라 실현된 모집단이다. \(\text{Exp}(1)\)의 왜도는 \(2\), 초과첨도는 \(6\)인데, 거기서 뽑힌 1만 개는 왜도 \(1.9196\), 초과첨도 \(5.1090\)밖에 되지 않는다. 최대값이 \(8.99\)에서 끊겼기 때문이다. \(\text{Exp}(1)\)에서 \(9\)를 넘는 값은 확률 \(e^{-9} = 0.00012\)로 1만 개에 한 개꼴이니 꼬리가 통째로 잘려 나간 셈이고, 적률의 차수가 높을수록 그 손실이 크다. 평균은 \(1.1\%\), 분산은 \(3.8\%\) 어긋나는데 4차 중심적률은 \(9\)여야 할 것이 \(7.507\)로 \(17\%\)가 사라졌다.

잘린 모집단에서 나오는 겨냥값은 \(0.8585\)와 \(1.0218\)이다. 이것과 견주면 모의값 \(0.8353\), \(0.9186\)이 각각 요동의 \(0.7\)배, \(0.6\)배만큼 떨어져 있어 깔끔하게 맞는다. 어긋난 것은 모의실험이 아니라 "모집단 1만 개를 \(\text{Exp}(1)\) 자체로 여긴" 가정이었다.

꼬리확률. 정확분포와 정규근사의 누적분포 차는 최대 \(0.0596\)으로 균등의 \(0.0057\)보다 열 배 크다. \(x = 1.5\)에서는 둘이 \(0.1321\)과 \(0.1318\)로 거의 같은데 이것은 우연이며, \(x = 2\)로 가면 \(0.0293\) 대 \(0.0127\)로 벌어진다. 모의값 \(0.1215\)와 \(0.0239\)는 잘린 모집단에서 뽑은 것이라 겨냥값이 각각 \(0.12557\)과 \(0.02627\)로 내려가 있고(같은 모집단에서 200만 번 다시 뽑아 확인했다), 거기서 비율의 몬테카를로 오차 \(0.0034\)와 \(0.0016\)의 \(1.2\)배, \(1.5\)배만큼 떨어져 있다.

요점은 셋째 줄이다. 폭은 \(\sigma/\sqrt n\)으로 정확히 주면서 꼬리는 \(2.3\)배 틀리는 것이 \(n = 5\)의 지수 모집단이고, 이것이 "\(n \ge 30\)이면 된다"는 어림짐작이 어디서 깨지는지 알려 준다.

적률생성함수는 \(t < \lambda\)에서 \(M_X(t) = \frac{\lambda}{\lambda - t}\)이다.

1차 적률:

2차 적률:

분산:

\(\square\)

\(X_i \sim \text{Exp}(\lambda)\)의 MGF는 \(M_{X_i}(t) = \frac{\lambda}{\lambda - t}\)이다.

독립성에 의해 합의 MGF는:

이는 \(\text{Gamma}(n, \lambda)\) 분포(형상 \(n\), 비율 \(\lambda\))의 MGF이다. MGF가 분포를 유일하게 결정하므로 \(S_n \sim \text{Gamma}(n, \lambda)\)이다. \(\square\)

\(S_5 \sim \text{Gamma}(5, 1)\)일 때 \(\bar{X} = S_5 / 5\)이므로 \(P(\bar{X} > 2) = P(S_5 > 10)\)이다.

Python으로:

from scipy import stats
# 정확한 값: 지수표본의 합은 감마분포를 따른다.
p_exact = 1 - stats.gamma.cdf(10, a=5, scale=1)
# 정규근사. 평균 1, 표준오차 1/sqrt(5) 이다.
p_normal = 1 - stats.norm.cdf(2, loc=1, scale=1/5**0.5)
print(f"Exact (Gamma):       {p_exact:.6f}")
print(f"Normal approximation: {p_normal:.6f}")

출력:

Exact (Gamma):       0.029253
Normal approximation: 0.012674

정확한 값은 약 0.0293이고 정규근사는 약 0.0127이다. \(n = 5\)는 중심극한정리가 지수분포의 치우침을 온전히 보정하기에 너무 작아, 정규근사가 오른쪽 꼬리 확률을 크게 과소추정한다. \(\square\)

왜도가 \(\gamma_1\)인 i.i.d. 확률변수에 대해 \(\bar{X} = \frac{1}{n}\sum X_i\)의 왜도는:

\(\bar{X}\)의 3차 중심적률이 (독립인 복사본 \(n\)개를 더하고 \(n\)으로 나누므로) \(\mu_3 / n^2\)이고 \((\text{Var}(\bar{X}))^{3/2} = (\sigma^2/n)^{3/2}\)이므로:

\(\text{Exp}(1)\)에서 \(\gamma_1 = 2\)이므로 \(\gamma_1(\bar{X}) = 2/\sqrt{n}\)이다.

\(2/\sqrt{n} < 0.5\)로 두면:

따라서 왜도가 0.5 아래로 떨어지려면 \(n \ge 17\)이 필요하다. \(\square\)

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

np.random.seed(1)
population = np.random.exponential(size=10_000)

sample_means = [
    np.mean(np.random.choice(population, size=50, replace=False))
    for _ in range(10_000)
]

fig, ax = plt.subplots(figsize=(10, 4))
ax.hist(sample_means, bins=60, density=True, alpha=0.5, label="Simulated")

x = np.linspace(min(sample_means), max(sample_means), 200)
ax.plot(x, stats.norm.pdf(x, loc=1, scale=1/np.sqrt(50)), "r--", lw=2,
        label="N(1, 1/50)")
ax.legend()
ax.set_title("Sampling Distribution of X-bar, n=50")
plt.show()

\(n = 50\)에서는 표본평균의 히스토그램이 거의 대칭이고 \(N(1, 1/50)\) 밀도를 바짝 따라간다. \(\bar{X}\)의 왜도가 \(2/\sqrt{50} \approx 0.28\)로 충분히 작아 정규근사가 아주 잘 맞는다. \(\square\)

최대가능도. 로그가능도가 \(\ell(\lambda) = n\ln\lambda - \lambda\sum x_i\)이므로

이다.

기대값. 연습문제 2에서 \(S_n = \sum X_i \sim \text{Gamma}(n,\lambda)\)이므로

이다(\(n \ge 2\)에서만 수렴한다). 따라서

이다. \(\square\)

항상 과대추정한다. \(1/x\)가 볼록함수이므로 옌센 부등식에서 예상되는 방향이다. 편향의 크기는 \(\lambda/(n-1)\)로, \(n=5\)면 25%나 된다.

불편 보정.

로 두면 불편이 된다. \(n=1\)이면 기대값이 존재하지 않으므로 보정 자체가 불가능하다.

이 예가 보여 주는 것은 최대가능도추정이 불변성은 가져도 불편성은 갖지 않는다는 점이다. \(\mu = 1/\lambda\)의 MLE는 \(\bar X\)이고 이것은 불편인데, 역수를 취하는 순간 편향이 생긴다. 어느 모수화로 문제를 적느냐에 따라 불편성이 나타났다 사라진다.

표본평균. \(\operatorname{Var}(X) = 1/\lambda^2\)이므로

중앙값 기반. 표본중앙값의 점근분산은 일반적으로

이다. 지수분포에서 중앙값은 \(m = \ln2/\lambda\)이고 그 점의 밀도는

이므로 \(\operatorname{Var}(\tilde X) \approx 1/(n\lambda^2)\)이다. 평균의 추정량은 \(\tilde X/\ln2\)이므로

비교. 상대효율이

로, 중앙값 기반 추정량의 분산이 두 배 넘게 크다. 같은 정밀도를 얻으려면 표본이 2.08배 필요하다.

뜻. 모형이 정말 지수분포라면 표본평균을 써야 한다. \(\bar X\)가 충분통계량의 함수이자 최소분산불편추정량이기 때문이다.

그래도 중앙값을 쓸 이유가 있다면 모형이 틀렸을 가능성이다. 자료에 오염이 섞이거나 꼬리가 지수보다 두꺼우면 \(\bar X\)는 몇 개의 큰 값에 끌려가지만 중앙값은 버틴다. 효율 52%를 보험료로 내고 강건성을 사는 거래이며, 고장시간 자료처럼 극단값의 신뢰도가 낮은 경우에는 합리적인 선택일 수 있다.

\(\hat\lambda = 1/2.5 = 0.4\)명/분, \(\hat\mu = 1/2.0 = 0.5\)명/분이므로

이다. 서비스 자체는 2분인데 대기까지 합치면 10분이다. 이용률 \(\rho = \lambda/\mu = 0.8\)의 효과다.

문제점 1 — 극도로 불안정하다. \(\hat W\)는 두 추정값의 차이의 역수다. \(\hat\mu - \hat\lambda = 0.1\)인데 각 추정값의 표준오차가

수준이다. 분모가 자기 표준오차보다 조금 큰 정도에 지나지 않는다. 델타 방법으로 계산하면

으로 추정값 자체보다 크다. 사실상 아무것도 말하지 못한다.

문제점 2 — 분모가 음수가 될 수 있다. \(\hat\mu < \hat\lambda\)가 나오면 \(\hat W\)가 음수가 되고, 이는 "대기행렬이 발산한다"는 뜻이다. 표본이 작으면 실제로 \(\rho<1\)인데도 이런 일이 일어난다.

문제점 3 — 비선형 변환의 편향. \(1/(\mu-\lambda)\)가 볼록이므로 \(\hat W\)는 \(W\)를 체계적으로 과대추정한다.

교훈. 혼잡한 시스템(\(\rho \to 1\))의 성능 지표는 추정하기 대단히 어렵다. 작은 추정오차가 결과를 폭발시킨다. 실무에서는 \(\rho\)를 직접 추정해 신뢰구간을 만들고, \(W\)는 부트스트랩이나 베이즈 방법으로 비대칭 구간을 보고하는 편이 낫다.

(가) 정확한 구간. \(2n\lambda\bar X = 2\lambda S_n \sim \chi^2_{2n}\)이 추축량이다. \(2n = 10\)이므로

에서 \(\mu = 1/\lambda\)에 대해 풀면

이다. 포함확률이 정확히 0.95이며, \(\bar x\)를 중심으로 극도로 비대칭이다. 위로 3배까지 뻗는다.

(나) 정규근사 구간. \(\operatorname{SE}(\bar X) = \mu/\sqrt n\)을 \(\bar x/\sqrt5\)로 추정하면

이다. 실제 포함확률은 \(P(0.124\bar X \le 1 \le 1.877\bar X)\), 즉 \(P(0.533 \le \bar X \le 8.06)\)인데 \(\bar X \sim \text{Gamma}(5, \text{척도}=1/5)\)이므로 계산하면 약 0.868이다.

비교.

정규근사는 명목 95%에서 실제로는 87%밖에 담지 못한다. 게다가 어긋남이 한쪽으로 몰려 있다. 위쪽 한계가 너무 낮아 참값이 위로 빠져나가는 경우가 대부분이다.

원인은 두 가지다. 첫째, \(n=5\)에서 \(\bar X\)의 왜도가 \(2/\sqrt5 = 0.894\)로 여전히 크다. 둘째, 표준오차를 \(\mu\)가 아니라 \(\bar x\)로 추정하면서 \(\bar x\)가 작게 나온 표본에서 구간이 함께 좁아지는 이중의 문제가 생긴다.

결론. 정확한 추축량을 알 수 있으면 반드시 그것을 쓴다. 지수·감마·포아송처럼 카이제곱 관계가 있는 분포에서는 정확한 구간이 공짜로 얻어진다.

최솟값. \(P(\min > t) = \{e^{-\lambda t}\}^n = e^{-n\lambda t}\)이므로

최댓값. \(P(\max \le t) = (1-e^{-\lambda t})^n\)이고, 순서통계량의 기대값 공식에서

이다(\(\gamma \approx 0.5772\)는 오일러 상수). 이는 최댓값을 \(\min\)부터 차례로 쌓아 올리는 분해에서 나오며, 무기억성 덕분에 \(k\)번째 구간이 \(\text{Exp}((n-k+1)\lambda)\)를 따른다.

세 통계량의 움직임.

셋의 성격이 완전히 다르다.

• \(\min\)은 \(1/n\)의 속도로 0에 붙는다. 산포도 같은 속도로 줄어 상대적 불확실성은 그대로다.
• \(\bar X\)는 큰수의 법칙대로 모평균에 수렴하고 산포가 \(1/\sqrt n\)로 준다. 유일하게 안정된 것이다.
• \(\max\)는 \(\ln n\)으로 천천히 발산하지만 산포가 줄지 않는다. 실제로 \(\lambda\max - \ln n\)이 굼벨분포로 수렴한다. 표본을 아무리 늘려도 최댓값의 불확실성은 그대로라는 뜻이다.

극단값을 다룰 때 표본평균의 직관이 통하지 않는 이유가 여기 있다. 평균은 자료가 쌓이면 안정되지만 최댓값은 그렇지 않다. 100년 홍수위나 최대 손실액 추정이 어려운 근본적인 까닭이며, 그래서 극단값 이론이라는 별도의 분야가 있다.