S²의 표본분포 (Exponential)¶

(1) 이론값. \(\text{Exp}(1)\)은 \(\sigma^2 = 1\), \(\mu_3 = 2\), \(\beta_2 = 9\)다. 불편성은 치우침과 무관하므로 두 경우 모두 \(E[S^2] = 1\)이다. 폭은

에 넣는다.

카이제곱 예측은 \(\sqrt{2/(n-1)}\)이므로 \(0.471405\)와 \(0.142134\)다. 분산 배율은

이고 극한값 \(4\)에 아래에서 다가간다. 표준편차로는 두 배 가까이 넓다.

상관계수는 \(\operatorname{Cov}(\bar X, S^2) = \mu_3/n\), \(\operatorname{Var}(\bar X) = \sigma^2/n\)을 쓰면 \(n\)이 약분되어

가 된다. \(n = 10\)에서 \(2/\sqrt{8.22222} = 0.69749\), \(n = 100\)에서 \(2/\sqrt{8.02020} = 0.70621\)이고 극한이 \(2/\sqrt{8} = 0.70711\)이다. \(n\)이 들어 있던 자리가 모두 지워진 것이 요점이다.

(2) 모의실험.

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

rng = np.random.default_rng(1)

beta2, mu3, sigma2_true = 9.0, 2.0, 1.0

fig, axes = plt.subplots(1, 2, figsize=(12, 3.5))
for ax, n in zip(axes, (10, 100)):
    samples = rng.exponential(size=(100_000, n))
    s2 = samples.var(axis=1, ddof=1)
    xbar = samples.mean(axis=1)

    # 이론값과 모의값을 나란히 적는다. 여기서 MC오차는 정규가정 아래의 어림값이다.
    sd_true = np.sqrt((beta2 - (n - 3) / (n - 1)) * sigma2_true ** 2 / n)
    sd_chi2 = np.sqrt(2 / (n - 1)) * sigma2_true
    ratio_true = (n - 1) / (2 * n) * (beta2 - (n - 3) / (n - 1))
    corr_true = mu3 / (np.sqrt(sigma2_true) ** 3 * np.sqrt(beta2 - (n - 3) / (n - 1)))
    print(f"n = {n}")
    print(f"  E[S^2]   이론 {sigma2_true:.6f}   모의 {s2.mean():.6f}   (MC오차 {sd_true / np.sqrt(len(s2)):.6f})")
    print(f"  sd(S^2)  이론 {sd_true:.6f}   모의 {s2.std(ddof=1):.6f}")
    print(f"  카이제곱이 예측하는 sd = {sd_chi2:.6f}   분산 배율  이론 {ratio_true:.4f}   모의 {(s2.std(ddof=1) / sd_chi2) ** 2:.4f}")
    print(f"  corr(X-bar, S^2)  이론 {corr_true:.4f}   모의 {np.corrcoef(xbar, s2)[0, 1]:.4f}")

    # 꼬리가 아주 길다. n=10 에서는 S^2 이 18을 넘는 표본도 나오므로
    # 99.5백분위에서 자르지 않으면 가운데가 뭉개져 보이지 않는다.
    hi = np.percentile(s2, 99.5)
    _, bins, _ = ax.hist(s2, bins=60, range=(0, hi), density=True,
                         alpha=0.5, edgecolor="white", label=r"simulated $S^2$")

    df, sigma2 = n - 1, 1.0
    c = df / sigma2
    g = np.linspace(1e-6, hi, 300)
    ax.plot(g, stats.chi2(df).pdf(g * c) * c, "--r", lw=2, label=r"$\chi^2$-based PDF")

    ax.set_title(f"Exp(1),  n = {n}")
    ax.set_xlabel(r"$S^2$")
    ax.set_xlim(0, hi)

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

출력:

n = 10
  E[S^2]   이론 1.000000   모의 0.991682   (MC오차 0.002867)
  sd(S^2)  이론 0.906765   모의 0.898319
  카이제곱이 예측하는 sd = 0.471405   분산 배율  이론 3.7000   모의 3.6314
  corr(X-bar, S^2)  이론 0.6975   모의 0.6957
n = 100
  E[S^2]   이론 1.000000   모의 1.000781   (MC오차 0.000896)
  sd(S^2)  이론 0.283200   모의 0.283307
  카이제곱이 예측하는 sd = 0.142134   분산 배율  이론 3.9700   모의 3.9730
  corr(X-bar, S^2)  이론 0.7062   모의 0.7056

\(n = 100\)은 깔끔하게 맞는다. \(E[S^2] = 1.000781\)이 몬테카를로 오차 \(0.000896\)의 \(0.9\)배, \(\operatorname{sd}(S^2) = 0.283307\)이 이론값과 \(0.04\%\) 차이, 배율 \(3.9730\)이 이론 \(3.9700\)과 맞고, 상관 \(0.7056\)이 이론 \(0.7062\)와 맞는다.

\(n = 10\)에서는 몬테카를로 오차를 조심해야 한다. 모의 \(E[S^2] = 0.991682\)가 \(1\)에서 \(0.0083\) 떨어져 있어 위에 적힌 오차 \(0.002867\)의 \(2.9\)배다. 같은 일을 120 번 되풀이해 재어 보면 평균이 \(0.999753\)으로 맞고 요동의 실제 폭이 \(0.00273\)이니, 이 한 번이 운 나쁜 쪽에 떨어진 것이다.

\(\operatorname{sd}(S^2)\)에서는 한 가지가 더 있다. 모의 \(0.898319\)가 이론 \(0.906765\)보다 \(0.93\%\) 작은데, 표준편차의 몬테카를로 요동을 정규가정 공식 \(\mathrm{sd}/\sqrt{2B} = 0.00203\)으로 어림하면 \(4\)배 어긋난 것으로 보인다. 그러나 그 공식은 재는 대상이 정규일 때만 맞다. \(n = 10\)의 지수모집단에서 \(S^2\)은 첨도가 아주 큰 분포이므로 요동이 실제로 \(0.00604\)로 세 배 크고, 그 단위로는 \(1.4\)배다. 120 번의 평균은 \(0.906358\)로 이론값 \(0.906765\)와 맞는다. 정규용 오차 공식을 비정규 대상에 쓰면 어긋남을 과장해 읽게 된다는 것이 이 보기가 덤으로 주는 교훈이다.

배율과 상관계수는 두 \(n\)에서 모두 맞는다. 상관이 \(0.6957\)과 \(0.7056\)으로 \(n\)이 열 배 달라져도 거의 그대로인 것이 눈에 띈다. 이것이 "\(n\)을 키워도 고쳐지지 않는다"의 가장 간결한 증거다.

(1) 극한 포함률. 구간은 \(W = (n-1)S^2/\sigma^2\)이 \(\chi^2(n-1)\)의 \(2.5\)-\(97.5\) 백분위 \([\ell, u]\)에 들 때 참값을 덮는다. \(W\)의 평균은 정확히 \(n-1\)이지만 분산은

로 카이제곱이 믿는 \(2(n-1)\)의 네 배다. \(n\)이 크면 \(W\)가 정규에 가까워지고 \(\ell, u \approx (n-1) \pm 1.96\sqrt{2(n-1)}\)이므로, \(W\)를 자기 표준편차 \(\sqrt{2(n-1)r}\)로 재면 그 경계가 \(\pm 1.96/\sqrt{r}\)에 놓인다. 따라서

명목 \(95\%\)가 실제로는 \(67\%\)다. 균등모집단에서는 \(0.998\)로 지나치게 넓었으나 여기서는 반대쪽으로, 그리고 훨씬 더 심하게 어긋난다.

(2) 모의실험.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
sigma2 = 1.0             # Exp(1)의 참 분산
beta2 = 9.0              # Exp(1)의 첨도

print("명목 신뢰수준 95%")
for n in (10, 30, 100, 1000):
    s2 = rng.exponential(size=(100_000, n)).var(axis=1, ddof=1)
    lo, hi = stats.chi2(n - 1).ppf([0.025, 0.975])
    cover = np.mean(((n - 1) * s2 / hi <= sigma2) & (sigma2 <= (n - 1) * s2 / lo))

    # W = (n-1)S^2/sigma^2 의 평균·표준편차·왜도. 카이제곱이라면 왜도가 sqrt(8/(n-1)) 여야 한다.
    w = (n - 1) * s2 / sigma2
    r = (n - 1) / (2 * n) * (beta2 - (n - 3) / (n - 1))
    sd_w = np.sqrt(2 * (n - 1) * r)

    # 예측 두 가지. 정규 근사는 적률 둘만, 이동감마 근사는 왜도까지 맞춘다.
    z_lo, z_hi = (lo - (n - 1)) / sd_w, (hi - (n - 1)) / sd_w
    pred_normal = stats.norm.cdf(z_hi) - stats.norm.cdf(z_lo)
    k = 4 / stats.skew(w) ** 2
    b = sd_w / np.sqrt(k)
    a = (n - 1) - b * k
    pred_gamma = stats.gamma(k).cdf((hi - a) / b) - stats.gamma(k).cdf(max((lo - a) / b, 0))

    print(f"n = {n:>4}:  실제 포함률 {cover:.3f}   배율 r = {r:.4f}   "
          f"왜도(W) 모의 {stats.skew(w):.3f} 카이제곱 {np.sqrt(8 / (n - 1)):.3f}   "
          f"예측 정규 {pred_normal:.3f} 이동감마 {pred_gamma:.3f}")
print(f"극한: 2*Phi(1.96/sqrt(4)) - 1 = {2 * stats.norm.cdf(stats.norm.ppf(0.975) / 2) - 1:.4f}")

출력:

명목 신뢰수준 95%
n =   10:  실제 포함률 0.764   배율 r = 3.7000   왜도(W) 모의 3.332 카이제곱 0.943   예측 정규 0.670 이동감마 0.908
n =   30:  실제 포함률 0.716   배율 r = 3.9000   왜도(W) 모의 1.768 카이제곱 0.525   예측 정규 0.672 이동감마 0.715
n =  100:  실제 포함률 0.690   배율 r = 3.9700   왜도(W) 모의 0.969 카이제곱 0.284   예측 정규 0.673 이동감마 0.681
n = 1000:  실제 포함률 0.673   배율 r = 3.9970   왜도(W) 모의 0.296 카이제곱 0.089   예측 정규 0.673 이동감마 0.674
극한: 2*Phi(1.96/sqrt(4)) - 1 = 0.6729

극한은 맞는다. \(n = 1000\)의 실제 포함률 \(0.673\)이 (1)에서 구한 \(0.6729\)와 소수 셋째 자리까지 일치한다. 비율의 몬테카를로 오차가 \(\sqrt{0.673 \times 0.327/10^5} = 0.0015\)이니 들어맞는다고 할 수 있다.

작은 \(n\)에서는 정규 예측이 맞지 않는다. 둘째 열의 배율 \(r\)은 \(n = 10\)에서 이미 \(3.70\)으로 극한 \(4\)에 가까운데, 정규 예측은 네 \(n\)에서 \(0.670\)–\(0.673\)으로 거의 움직이지 않는다. 실제는 \(0.764 \to 0.716 \to 0.690 \to 0.673\)으로 내려간다. 빠진 것은 \(r\)이 아니라 모양이다.

셋째 열이 그것을 가리킨다. \(n = 10\)에서 \(W\)의 왜도가 \(3.332\)인데 카이제곱이라면 \(\sqrt{8/9} = 0.943\)이어야 한다. 세 배 반이나 더 치우쳐 있다. 오른쪽으로 길게 늘어진 분포는 질량 대부분이 평균 왼쪽에 몰려 있고, 구간 \([\ell, u]\)는 평균을 기준으로 왼쪽 \(-0.77\) 표준편차, 오른쪽 \(+1.23\) 표준편차에 걸쳐 있다. 왼쪽 경계가 짧아 거기서 조금 잃지만, 질량이 몰려 있는 중앙 왼쪽을 구간이 품고 있어 전체로는 정규가 예측한 것보다 더 많이 덮는다.

왜도를 맞추면 따라잡히는지 확인하려고 마지막 열에 이동감마 근사를 두었다. \(W \approx a + b\,\text{Gamma}(k)\)로 놓고 평균·분산·왜도 셋을 맞춘 것이다. \(n = 30\)에서 \(0.715\)(실제 \(0.716\)), \(n = 100\)에서 \(0.681\)(실제 \(0.690\)), \(n = 1000\)에서 \(0.674\)(실제 \(0.673\))로 정규 예측보다 뚜렷이 낫다. 왜도를 넣은 것이 \(n\) 의존성의 대부분을 설명한다.

다만 \(n = 10\)에서는 이동감마가 \(0.908\)을 주어 실제 \(0.764\)를 크게 넘어선다. 적률 셋으로도 모자란다는 뜻이고, 그 자리에서는 왜도 \(3.33\)에 대응하는 감마 형상모수가 \(k = 4/3.33^2 = 0.36\)으로 극단적이어서 근사 자체가 믿을 만하지 않다. 세 적률을 맞춘 근사가 두 적률보다 낫지만 \(n = 10\)에서는 아직 모자라다는 것을 숨기지 않고 적어 둔다.

어느 쪽이든 결론은 같다. \(0.764\)에서 \(0.673\)으로 가는 내림세는 고쳐지는 방향이 아니다. 카이제곱이 작은 \(n\)에서 우연히 덜 틀려 보이던 것이 \(n\)이 커지며 벗겨지는 것이고, 바닥은 \(0.673\)이다.

\(\mu = 1\)이므로 중심적률을 원점적률로 전개한다. \(E[X^k] = k!\)이므로 \(E[X]=1\), \(E[X^2]=2\), \(E[X^3]=6\), \(E[X^4]=24\)이고

이다. \(\sigma^2 = 1\)이므로 \(\beta_2 = 9\), 왜도는 \(\mu_3/\sigma^3 = 2\)다.

\(\lambda\)에 의존하지 않는다. 왜도와 첨도는 척도에 불변이고 \(\text{Exp}(\lambda)\)는 \(\text{Exp}(1)\)을 \(1/\lambda\)배 한 것이기 때문이다. 어떤 비율모수를 쓰든 배율은 언제나 4다.

(a) 참값. \(\beta_2 = 9\), \(\sigma^4 = 1\), \(n = 100\)이므로

(b) 카이제곱 예측.

비가 \(0.283/0.142 = 1.99\)로 두 배다. 모의실험에서 얻은 \(\text{Var}(S^2) = 0.0805\)와도 맞는다.

구간 폭. 카이제곱 구간은 실제 필요한 폭의 절반이다. 거꾸로 말하면 제대로 된 구간은 카이제곱 구간보다 두 배 넓어야 한다. 포함률이 0.95에서 0.69로 떨어지는 크기가 이 정도 왜곡에서 나온다.

\(n\)이 크면 두 분포 모두 정규로 근사된다. 참 분포는

이고, 카이제곱 구간은 \(S^2\)의 표준편차를 \(\text{SD}_{\chi^2} = \sigma^2\sqrt{2/n}\)로 믿고 \(\pm 1.96\,\text{SD}_{\chi^2}\)만큼 뻗는다. 따라서 구간이 참 표준편차 단위로는

만큼만 뻗는다. 그러므로 포함률의 극한은

이다. 보기 2의 \(n = 1000\) 값 0.673과 소수 셋째 자리까지 맞는다. \(\square\)

이 계산은 일반적인 공식을 준다. 명목 신뢰수준 \(1-\alpha\)에서 극한 포함률은

이다. 균등분포(\(\beta_2 = 1.8\))를 넣으면 \(2\Phi(1.96\sqrt{2.5}) - 1 = 2\Phi(3.10) - 1 = 0.998\)로 앞 페이지의 값이 나온다. 첨도 하나만 알면 카이제곱 구간이 얼마나 망가지는지 미리 계산할 수 있다.

\(\text{Cov}(\bar X, S^2) = \mu_3/n\)이고 \(\text{Var}(\bar X) = \sigma^2/n\), \(\text{Var}(S^2) \approx (\beta_2-1)\sigma^4/n\)이므로

이다. 분자와 분모가 모두 \(1/n\)이라 \(n\)이 약분된다. 지수분포는 \(\mu_3 = 2\), \(\sigma = 1\), \(\beta_2 = 9\)이므로 \(2/\sqrt 8 = 0.707\)이다. \(\square\)

분자의 \(\mu_3/\sigma^3\)은 왜도이므로, 이 상관계수는 왜도를 \(\sqrt{\beta_2-1}\)로 나눈 것이다. 모집단이 치우친 만큼 \(\bar X\)와 \(S^2\)이 얽힌다.

\(t\) 통계량에 미치는 영향. \(t = (\bar X - \mu)/(S/\sqrt n)\)에서 분자가 크면 분모도 함께 커진다. 오른쪽으로 치우친 모집단에서 \(\bar X\)가 우연히 크게 나온 표본은 \(S\)도 크므로 \(t\)가 그만큼 커지지 못하고, 반대로 \(\bar X\)가 작으면 \(S\)도 작아 \(t\)가 더 음수 쪽으로 간다. 그 결과 \(t\)의 분포가 왼쪽으로 치우친다. 단측검정의 실제 유의수준이 한쪽은 명목보다 높고 다른 쪽은 낮아진다.

치우친 자료에서 \(t\) 검정의 양측 오류율이 그럭저럭 맞는 것도 이 때문이다. 두 꼬리의 오차가 서로 상쇄되기 때문이며, 단측검정에서는 상쇄가 일어나지 않아 문제가 드러난다.

(1) 모집단이 지수라는 것을 안다면 정확한 방법이 있다. 4장 지수분포 연습문제에서 본 \(2\lambda\sum_i X_i \sim \chi^2_{2n}\)을 쓰면 \(\lambda\)의 정확한 구간을 얻고, 지수분포에서는 \(\sigma^2 = 1/\lambda^2\)이므로 변환해서 분산의 구간으로 옮길 수 있다. 전제가 강하다(분포족을 알아야 한다).

(2) 첨도를 추정해 정규근사를 쓴다.

전제는 \(n\)이 충분히 커서 \(S^2\)의 정규근사가 통하고 \(\hat\beta_2\)가 안정적이라는 것이다. 문제는 \(\hat\beta_2\) 자체가 4차 적률 추정이라 매우 불안정하다는 점이다. 꼬리가 두꺼운 자료에서 표본첨도는 표본크기에 따라 계속 커지기도 한다.

(3) 부트스트랩. \(S^2\)의 표집분포를 재표집으로 직접 추정한다. 분포족도 첨도 공식도 필요 없고, 전제는 표본이 모집단을 대표한다는 것과 4차 적률이 유한하다는 것뿐이다. 5장의 부트스트랩 표준오차 페이지에서 다룬다. 실무의 기본 선택이다.

한 가지 덧붙이면, 애초에 분산 자체가 관심사인지 되물을 필요가 있다. 치우친 자료에서는 분산보다 사분위범위나 로그 변환 후의 산포가 더 뜻 있는 요약일 때가 많다.

"\(n \ge 30\)" 규칙은 \(\bar X\)에만 적용된다. 이 규칙은 중심극한정리의 수렴 속도에 관한 경험칙이고, \(S^2\)의 카이제곱 관계는 수렴의 문제가 아니라 전제의 문제다. 정규성이 깨지면 \(n\)을 늘려도 다른 분포로 수렴하는 것이 아니라 애초에 틀린 분포를 쓰고 있는 것이다.

한 걸음 더 나아가면, 규칙 자체도 \(\bar X\)에 대해서조차 무비판적으로 쓸 수 없다. 4장에서 본 것처럼 지수모집단의 합은 왜도가 \(2/\sqrt n\)로 줄어들어 \(n = 30\)에서 0.37이 남고, 꼬리 확률을 다루는 경우에는 그 정도로도 부족하다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
for n in (10, 30, 100):
    S2 = rng.exponential(1, (200_000, n)).var(1, ddof=1)
    print(f"  n={n:>4}: S^2 왜도 {stats.skew(S2):>8.4f}   "
          f"log S^2 왜도 {stats.skew(np.log(S2)):>8.4f}")

출력:

  n=  10: S^2 왜도   2.9497   log S^2 왜도  -0.1872
  n=  30: S^2 왜도   1.7490   log S^2 왜도  -0.0274
  n= 100: S^2 왜도   0.9636   log S^2 왜도   0.0366

로그가 왜도를 거의 완전히 없앤다. \(n=30\)에서 \(1.75\)가 \(-0.027\)로 줄었다. \(n=10\)에서도 \(2.95 \to -0.19\)로, 원래 척도의 \(n=100\)(\(0.96\))보다 훨씬 대칭이다.

왜 이렇게 잘 듣는가. \(S^2\)은 본질적으로 제곱합이고 제곱합은 곱셈적으로 변동한다. "참값의 \(1.5\)배"와 "참값의 \(1/1.5\)배"가 대칭적인 사건인데, 원래 척도에서는 \(+0.5\sigma^2\)과 \(-0.33\sigma^2\)로 비대칭하게 보인다. 로그를 취하면 \(\pm0.405\)로 대칭이 된다. 분산은 덧셈이 아니라 곱셈의 눈금 위에 있는 양이다.

그래서 구간은 로그 척도에서 만든다.

를 만든 뒤 지수를 취해 되돌린다. 지수함수가 단조이므로 포함확률이 보존되고, 되돌린 구간은 자동으로 양수이며 비대칭이 된다.

이 방법의 장점 셋.

세 번째가 특히 유용하다. \(\operatorname{sd}(\log S^2)\)가 \(\sigma\)에 의존하지 않으므로, 표준오차를 구하는 데 \(\sigma\)의 추정값이 필요 없다. 앞 절 \(F\) 분포 연습문제 3에서 \(\operatorname{Var}(\log S^2) \approx 2/(n-1)\)을 쓴 것도 같은 이유였다.

다만 \(\beta_2\)는 여전히 필요하다. 지수모집단이면 \(\beta_2 = 9\)이므로 \(\operatorname{sd}(\log S^2) \approx \sqrt{8/n}\)이고, 정규라면 \(\sqrt{2/n}\)이다. 로그 변환이 정규성 가정을 없애 주지는 않는다. 왜도 문제만 고칠 뿐이며, 첨도를 모르면 여전히 연습문제 5의 방법(부트스트랩 등)이 필요하다.

import numpy as np

rng = np.random.default_rng(0)
print(f"{'n':>5}{'S^2 의 표준오차':>16}{'Xbar^2 의 표준오차':>20}{'효율비':>10}")
for n in (10, 30, 100):
    X = rng.exponential(1, (200_000, n))
    a, b = X.var(1, ddof=1), X.mean(1) ** 2
    print(f"{n:>5}{a.std():>16.5f}{b.std():>20.5f}{(a.std()/b.std())**2:>10.3f}")

출력:

    n      S^2 의 표준오차       Xbar^2 의 표준오차       효율비
   10         0.90631             0.71051     1.627
   30         0.52044             0.38089     1.867
  100         0.28363             0.20299     1.952

\(\bar X^2\)이 훨씬 정밀하다. 효율비가 \(n\)이 커지면서 \(2\)에 수렴한다. 같은 정밀도를 얻는 데 \(S^2\)은 표본이 두 배 필요하다는 뜻이다.

왜 이득이 생기는가. \(S^2\)은 모집단이 무엇인지 모른다고 가정하고 2차 적률만으로 분산을 추정한다. 반면 \(\bar X^2\)은 "지수분포다"라는 정보를 쓴다. 지수분포에서는 평균 하나가 분포 전체를 결정하므로, 평균을 잘 추정하면 분산도 공짜로 잘 추정된다. 1차 적률이 2차 적률보다 훨씬 안정적이라는 것이 이 절의 주제였고, 그 안정성을 분산 추정에 끌어온 셈이다.

대가는 모형 오설정 위험이다.

두 번째 줄이 핵심이다. 모집단이 지수가 아니라 감마\((2, \theta)\)라면 \(\sigma^2 = 2\theta^2\)이지만 \(\mu^2 = 4\theta^2\)이므로, \(\bar X^2\)은 참 분산의 두 배로 수렴한다. 표본을 아무리 키워도 틀린 값에 수렴하며, 이것은 편향이 아니라 비일치성이라 회복할 방법이 없다.

이것이 모수적 방법과 비모수적 방법의 일반적인 맞바꿈이다. 모형을 가정하면 정보가 늘어 효율이 오르지만, 가정이 틀리면 답 자체가 틀린다. \(S^2\)은 효율이 낮은 대신 어떤 모집단에서도 \(\sigma^2\)을 향해 간다.

실무에서는 진단을 먼저 한다. 지수성을 믿을 근거가 있는지(변동계수가 \(1\)에 가까운지, Q-Q 그림이 직선인지) 확인하고, 확신이 서면 모수적 추정을, 아니면 \(S^2\)을 쓴다. 절충안으로 감마분포처럼 한 단계 넓은 족을 가정하는 것도 있으며, 지수를 특수한 경우로 포함하므로 위험이 줄어든다.

이론값. \(\text{Exp}(1)\)의 중앙값은 \(\ln 2 = 0.6931\)이다. MAD는 \(|X - \ln 2|\)의 중앙값이므로 \(P(|X-\ln 2| \le m) = 0.5\)를 푼다. 사분위수는 \(Q_1 = -\ln(0.75) = 0.2877\), \(Q_3 = -\ln(0.25) = 1.3863\)이다.

import numpy as np
from scipy import stats
from scipy.optimize import brentq

rng = np.random.default_rng(0)
med = np.log(2)
f = lambda m: (stats.expon.cdf(med + m) - stats.expon.cdf(max(med - m, 0))) - 0.5
mad_th = brentq(f, 1e-6, 5)

X = rng.exponential(1, (100_000, 50))
mad = np.median(np.abs(X - np.median(X, axis=1, keepdims=True)), axis=1)
iqr = np.percentile(X, 75, axis=1) - np.percentile(X, 25, axis=1)

print(f"  이론 MAD = {mad_th:.4f},  모의 {mad.mean():.4f}")
print(f"  이론 IQR = {stats.expon.ppf(.75) - stats.expon.ppf(.25):.4f},  "
      f"모의 {iqr.mean():.4f}")
print(f"  1.4826 * MAD = {1.4826 * mad_th:.4f}   (참 sigma = 1)")
print(f"  IQR / 1.349  = {(stats.expon.ppf(.75)-stats.expon.ppf(.25))/1.349:.4f}")

출력:

  이론 MAD = 0.4812,  모의 0.4738
  이론 IQR = 1.0986,  모의 1.0710
  1.4826 * MAD = 0.7134   (참 sigma = 1)
  IQR / 1.349  = 0.8144

환산상수를 그대로 쓰면 \(\sigma\)를 \(29\%\) 과소추정한다. \(1.4826 \times \text{MAD} = 0.713\)인데 참값은 \(1\)이다. IQR 쪽도 \(0.814\)로 \(19\%\) 낮다.

이유는 상수가 정규분포 전용이기 때문이다. \(1.4826 = 1/\Phi^{-1}(0.75)\)는 "정규분포에서 MAD를 \(\sigma\)로 바꾸는 환산율"이다. 지수분포는 꼬리가 훨씬 두꺼워 같은 \(\sigma\)에 대해 중앙 부분이 더 좁으므로, 정규용 상수로는 모자란다.

지수분포용 상수는 \(1/0.4812 = 2.078\)이다. 모집단을 알면 상수를 다시 구하면 된다. 그런데 모집단을 알면 애초에 강건 측도를 쓸 이유가 줄어든다(연습문제 8). 여기에 강건 측도의 근본적인 난점이 있다.

그래서 강건 측도는 "\(\sigma\)의 추정"보다 "산포 자체의 요약"으로 쓰는 것이 정직하다. 두 집단의 MAD를 견주는 것은 뜻이 있지만, MAD에서 \(\sigma\)를 역산해 정규 이론에 넣는 것은 가정을 몰래 들여오는 일이다.

덤으로 하나. 지수분포의 변동계수는 정확히 \(1\)이다(\(\sigma = \mu\)). 위 모의에서 \(S/\bar X\)가 \(0.98\)로 나오는데, 이 값이 \(1\)에서 크게 벗어나면 지수 가정을 의심할 수 있다. 상수 환산이 필요 없는 진단이라 실무에서 유용하다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, B, REP = 30, 999, 2000
cover_chi = cover_boot = 0
width_chi = width_boot = 0.0

for _ in range(REP):
    x = rng.exponential(1, n)
    s2 = x.var(ddof=1)
    # (1) 카이제곱 구간
    lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
    hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)
    cover_chi += (lo <= 1 <= hi); width_chi += hi - lo
    # (2) 부트스트랩 백분위 구간
    idx = rng.integers(0, n, (B, n))
    bs = x[idx].var(axis=1, ddof=1)
    blo, bhi = np.percentile(bs, [2.5, 97.5])
    cover_boot += (blo <= 1 <= bhi); width_boot += bhi - blo

print(f"  n={n}, 지수모집단, 참 sigma^2 = 1")
print(f"  카이제곱  포함률 {cover_chi/REP:.3f}   평균 폭 {width_chi/REP:.3f}")
print(f"  부트스트랩 포함률 {cover_boot/REP:.3f}   평균 폭 {width_boot/REP:.3f}")

출력:

  n=30, 지수모집단, 참 sigma^2 = 1
  카이제곱  포함률 0.719   평균 폭 1.166
  부트스트랩 포함률 0.744   평균 폭 1.381

부트스트랩이 거의 도움이 되지 않는다. 포함률이 \(0.719\)에서 \(0.744\)로 겨우 \(0.025\) 오르는 데 그치고, 둘 다 명목값 \(0.95\)에서 한참 멀다. 게다가 구간은 오히려 \(18\%\) 넓어졌다(\(1.166 \to 1.381\)). 폭을 더 쓰고도 포함률을 거의 못 산 셈이다.

이 결과는 예상보다 나쁘다. 부트스트랩은 분포 가정을 없애 주므로 크게 나아질 것 같은데 그렇지 않다.

왜 부트스트랩도 부족한가. 부트스트랩은 모집단의 모양을 표본에서 읽어 오므로 원리적으로는 첨도가 \(9\)라는 사실을 반영할 수 있다. 카이제곱이 \(\beta_2 = 3\)을 강요하는 것과 달리 이 부분은 해결된다. 그러나

• \(n = 30\)의 표본이 지수분포의 꼬리를 제대로 담지 못한다. 4차 적률을 추정해야 하는데, 꼬리가 두꺼운 분포에서 4차 적률은 관측 몇 개가 좌우한다. 재표본은 원표본에 없는 큰 값을 만들어 내지 못하므로, 원표본이 꼬리를 과소대표하면 부트스트랩도 그대로 과소대표한다.
• 백분위 구간은 \(S^2\) 분포의 치우침을 보정하지 않는다. 연습문제 7에서 본 대로 \(S^2\)은 왜도 \(1.75\)로 심하게 치우쳐 있는데, 백분위 구간은 재표본 분포의 꼬리를 그대로 잘라 쓸 뿐이라 이 치우침이 구간을 엉뚱한 쪽으로 밀어낸다. 폭이 넓어졌는데도 포함률이 오르지 않은 것이 그 증상이다. 구간이 참값을 가운데 두지 못하고 한쪽으로 치우쳐 있다는 뜻이기 때문이다.

개선 방향 둘.

"부트스트랩은 만능이 아니다"가 요점이다. 분포 가정을 없애 주지만 표본이 모집단을 대표한다는 가정은 그대로 남는다. 그리고 4차 적률처럼 추정이 어려운 양에 기대는 문제에서는 그 가정이 특히 위태롭다.

일반 규칙 하나. 부트스트랩은 평균 같은 안정된 통계량에서 매우 잘 듣고, 분산·첨도·극값 같은 고차 또는 꼬리 의존 통계량에서는 덜 듣는다. 이 절이 보인 "\(\bar X\)는 괜찮고 \(S^2\)은 아니다"라는 구도가 부트스트랩에도 그대로 이어진다. 17장에서 부트스트랩의 작동 조건을 다룰 때 이 예가 다시 나온다. \(\square\)