분산 신뢰구간의 포함확률 모의실험¶

(1) 포함률은 잴 것이 못 된다. 추축량 \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\) 때문에 포함률은 정확히 \(0.95\)이고 이것은 \(n\)과 \(\mu\), \(\sigma\)에 전혀 의존하지 않는다. 모의실험이 \(0.95\)를 재현해도 새로 알게 되는 것이 없고, \(0.93\)이 나오면 그것은 되풀이 \(100\)회의 몬테카를로 잡음일 뿐이다.

잴 것은 폭이다. 포함률이 고정되어 있으므로 이 구간의 쓸모를 가르는 것은 얼마나 좁은가 하나뿐이다. 폭은

로 \(S^2\)에 정비례한다(\(\ell = \chi^2_{0.025,\,n-1}\), \(u = \chi^2_{0.975,\,n-1}\)). \(S^2\)이 \(\sigma^2\)의 불편추정량이므로 \(E[S^2] = \sigma^2\)이고, 따라서 기대 폭은 닫힌 꼴로 나온다.

\(n = 12\)이면 \(\ell = 3.81575\), \(u = 21.92005\) 이므로

이다. 추정하려는 값 \(\sigma^2 = 4\)의 \(2.38\)배다. 폭의 표준편차도 같은 식에서 나온다. \(\operatorname{SD}(S^2) = \sigma^2\sqrt{2/(n-1)}\) 이므로

로, 폭 자체가 크게 흔들린다.

(2) 모의실험.

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

np.random.seed(42)          # 아래 출력과 그림을 재현하려면 고정한다

n_simulations = 100
n_samples = 12
mu, sigma = 0.0, 2.0
alpha = 0.05

true_var = sigma**2
df = n_samples - 1
# 임계값은 표본과 무관하므로 반복문 밖에서 한 번만 구한다.
# 자유도만으로 정해지는 값이기 때문이다.
chi2_lo = chi2.ppf(alpha / 2, df=df)
chi2_hi = chi2.ppf(1 - alpha / 2, df=df)

lowers = np.empty(n_simulations)
uppers = np.empty(n_simulations)
centers = np.empty(n_simulations)

for i in range(n_simulations):
    x = np.random.normal(loc=mu, scale=sigma, size=n_samples)
    s2 = x.var(ddof=1)
    lowers[i] = df * s2 / chi2_hi
    uppers[i] = df * s2 / chi2_lo
    centers[i] = s2

covered = (lowers <= true_var) & (true_var <= uppers)
n_fail = int((~covered).sum())
coverage_pct = 100.0 * covered.mean()
print(f"Coverage: {coverage_pct:.1f}%, Failures: {n_fail}")

# (1) 잴 것은 폭이다. 기대 폭은 E[S^2] = sigma^2 에서 닫힌 꼴로 나온다.
widths = uppers - lowers
exp_w = df * true_var * (1 / chi2_lo - 1 / chi2_hi)
print(f"기대 폭 (이론)  {exp_w:.5f}    표준편차 (이론) {exp_w * np.sqrt(2 / df):.5f}")
print(f"평균 폭 (모의)  {widths.mean():.5f}    표준편차 (모의) {widths.std(ddof=1):.5f}")
print(f"폭 최소 {widths.min():.4f}  중앙값 {np.median(widths):.4f}  최대 {widths.max():.4f}")
print(f"참값 4 와 8 을 동시에 담는 구간 {int(((lowers <= 4) & (8 <= uppers)).sum())} 개")
print(f"참값 2 와 8 을 동시에 담는 구간 {int(((lowers <= 2) & (8 <= uppers)).sum())} 개")

출력:

Coverage: 95.0%, Failures: 5
기대 폭 (이론)  9.52386    표준편차 (이론) 4.06099
평균 폭 (모의)  9.35875    표준편차 (모의) 4.08747
폭 최소 2.5049  중앙값 9.0140  최대 24.7756
참값 4 와 8 을 동시에 담는 구간 69 개
참값 2 와 8 을 동시에 담는 구간 26 개

기대 폭이 맞는다. 이론값 \(9.52386\) 과 모의 평균 \(9.35875\) 의 차이는 \(0.165\)인데, 폭의 표준편차가 \(4.06\) 이므로 \(100\)개 평균의 몬테카를로 표준오차는 \(4.06/\sqrt{100} = 0.406\) 이다. \(0.41\) 표준오차 안이다. 폭의 표준편차도 이론 \(4.06099\) 대 모의 \(4.08747\) 로 맞는다.

포함률은 \(95.0\%\)가 나왔지만 이것은 확인이 아니다. (1)에서 말한 대로 참값이 \(0.95\)임을 이미 알고 있고, \(100\)회의 표준오차가 \(0.0218\) 이라 \(93\)개나 \(97\)개가 나왔어도 똑같이 "맞는다"고 했을 것이다.

폭이 말해 주는 것이 훨씬 아프다. 평균 폭 \(9.36\) 은 추정하려는 \(\sigma^2 = 4\) 의 \(2.3\)배이고, 가장 넓은 구간은 \(24.78\) 로 참값의 여섯 배가 넘는다. 구체적으로 보면

• \(100\)개 중 \(69\)개가 \(4\)와 \(8\)을 동시에 담는다. 곧 열에 일곱은 "분산이 두 배인가"조차 가르지 못한다.
• \(100\)개 중 \(26\)개가 \(2\)와 \(8\)을 동시에 담는다. 넷에 하나는 네 배 차이도 가르지 못한다.

포함률이 정확히 \(0.95\)라는 것과 구간이 쓸모 있다는 것은 전혀 다른 말이다. \(n = 12\)로 분산을 추정하는 일은 이 정도로 막연하며, 그 사실은 포함률이 아니라 폭에만 나타난다.

(1) 그림.

# 구간 하나를 가로선 하나로 그린다. 참값을 담은 구간은 검정, 놓친 구간은
# 빨강이다. 세로 점선이 참값이고, 빨간 선이 몇 개인지 세는 것이 곧 포함확률을
# 재는 일이다. 구간마다 길이가 다른 까닭은 표본마다 s 가 다르기 때문이다.
fig, ax = plt.subplots(figsize=(12, 12))
for i in range(n_simulations):
    color = "k" if covered[i] else "r"
    ax.plot([lowers[i], uppers[i]], [i, i], lw=2, color=color)
    ax.plot(centers[i], i, marker="o", ms=3, color=color)

ax.axvline(true_var, linestyle="--", linewidth=1.5, color="r")
ax.set_title(f"{n_simulations} Chi-square Variance CIs | n={n_samples}, CL=95%")
ax.set_yticks([])
ax.set_xlabel("Variance")
plt.tight_layout()
plt.show()

(1) 점이 왼쪽에 치우친다. 점은 \(s^2\)이고 구간의 두 끝은 \(s^2\)에 각각 \((n-1)/u\)와 \((n-1)/\ell\)을 곱한 값이다. \(n = 12\)에서 그 두 배율이

이다. 점에서 왼쪽 끝까지는 \(s^2\)의 \(1 - 0.50182 = 0.49818\)배, 오른쪽 끝까지는 \(2.88279 - 1 = 1.88279\)배다. 오른쪽이 \(3.78\)배 멀다. 어느 표본에서든 같은 비이므로 \(100\)개 선이 모두 똑같이 치우쳐 있다.

까닭은 \(\chi^2_{11}\)이 오른쪽으로 늘어져 있다는 것 하나다. \(S^2\)이 참값보다 크게 나오는 쪽의 폭이 넓으므로, 그것을 되돌려 \(\sigma^2\)을 가두려면 오른쪽에 자리를 더 내주어야 한다.

평균의 \(t\)-구간 그림과 두 가지가 다르다.

둘째 줄이 그림에서 눈에 띄는 모습이다. 오른쪽으로 멀리 간 선은 반드시 길고, 왼쪽에 몰린 선은 반드시 짧다. 평균의 그림에서는 길고 짧은 선이 참값 양쪽에 고르게 섞여 있었다.

(2) \(\sigma\) 척도에서는 비대칭이 줄되 사라지지 않는다. 양 끝에 제곱근을 취하면 배율도 제곱근이 되어

이다. 점(\(s\))에서 왼쪽까지가 \(s\)의 \(0.29160\)배, 오른쪽까지가 \(0.69788\)배로 오른쪽이 \(2.39\)배 멀다. \(\sigma^2\) 척도의 \(3.78\)배보다 완화되었지만 여전히 \(1\)이 아니다.

제곱근이 비대칭을 줄이는 것은 로그가 아니기 때문이다. 배율의 비는 \(\sigma^2\) 척도에서 \(2.88279/0.50182 = 5.745\), \(\sigma\) 척도에서 그 제곱근인 \(2.397\) 이다. 어떤 거듭제곱을 취해도 이 비는 \(1\)이 되지 않고 \(1\)에 가까워질 뿐이다. 비대칭은 척도의 문제가 아니라 \(\chi^2_{11}\) 자체가 치우쳐 있다는 사실의 문제다. \(n\)을 키워 \(\chi^2_{n-1}\)이 대칭에 가까워져야 비로소 사라진다.

\(\text{df} = 19\), \(\alpha = 0.05\)일 때:

\(\sigma^2\)의 95% 신뢰구간은

제곱근을 취하면 \(\sigma\)의 95% 신뢰구간은

\(\square\)

\(Z_i = (X_i - \mu)/\sigma\)라 하면 \(Z_i \overset{\text{iid}}{\sim} N(0,1)\)이다. 그러면

항등식 \(X_i - \mu = (X_i - \bar{X}) + (\bar{X} - \mu)\)로 분해하면:

\(\sigma^2\)으로 나누면:

둘째 항은 \(Z = (\bar{X}-\mu)/(\sigma/\sqrt{n}) \sim N(0,1)\)인 \(Z^2\)이므로 \(\chi^2_1\)을 따른다. (정규성 아래 \(\bar{X}\)와 \(S^2\)의 독립성, 즉 Cochran 정리에 의해) 두 항은 독립이다. 따라서 독립인 카이제곱 확률변수의 가법성에 의해

이다. \(\square\)

from scipy.stats import chi2
np.random.seed(0)
n, n_sim, true_var = 20, 10_000, 1.0
df = n - 1
chi2_lo = chi2.ppf(0.025, df)
chi2_hi = chi2.ppf(0.975, df)
covers = 0
for _ in range(n_sim):
    x = np.random.exponential(1.0, size=n)
    s2 = x.var(ddof=1)
    lo = df * s2 / chi2_hi
    hi = df * s2 / chi2_lo
    if lo <= true_var <= hi:
        covers += 1
print(f"Coverage: {100 * covers / n_sim:.1f}%")

출력:

Coverage: 72.7%

95%를 훨씬 밑돈다. 스무 번에 한 번 놓치리라 믿었는데 네 번에 한 번 놓친다. 지수분포는 오른쪽으로 치우쳐 있고 초과첨도가 6이므로 \((n-1)S^2/\sigma^2\)이 \(\chi^2_{n-1}\)을 따르지 않는다.

앞 페이지의 평균 신뢰구간과 비교해 보라. 같은 지수자료, 비슷한 표본크기에서 \(t\)-구간의 포함확률은 90% 근처였고 \(n\)을 키우면 95%로 올라왔다. 여기서는 그렇지 않다. 이 코드에서 \(n\)만 바꿔 보면

로 오히려 나빠지면서 67% 근처로 수렴한다.

이유는 분명하다. 평균의 구간은 중심극한정리가 받쳐 주므로 \(n\)이 커지면 어떤 분포에서든 맞아 들어간다. 반면 카이제곱 추축량은 정규성 자체에 의존하며, 표본을 아무리 늘려도 모집단이 정규가 되지는 않는다. \(\text{Var}(S^2)\)은 첨도에 따라 달라지는데 카이제곱 구간은 정규분포의 첨도를 가정해 너비를 정하므로, 첨도가 큰 자료에서는 구간이 처음부터 너무 좁다. \(n\)을 키우면 그 "너무 좁음"이 더 또렷해질 뿐이다. 비정규 자료의 분산 구간에는 붓스트랩을 쓰라. \(\square\)

하한은 \(L = (n-1)S^2 / \chi^2_{1-\alpha/2}\)이고 상한은 \(U = (n-1)S^2 / \chi^2_{\alpha/2}\)이다. \(S^2\)으로부터의 거리는:

\(\chi^2_{n-1}\) 분포가 오른쪽으로 치우쳐 있으므로 \(\chi^2_{\alpha/2} < n-1 < \chi^2_{1-\alpha/2}\)이다(\(\chi^2_{n-1}\)의 평균 \(n-1\)이 위쪽 분위수에 더 가깝다). 따라서 \((n-1)/\chi^2_{\alpha/2} > 1\)이고 \((n-1)/\chi^2_{1-\alpha/2} < 1\)인데, 전자가 1을 넘는 양이 후자가 1에 못 미치는 양보다 커서 \(U - S^2 > S^2 - L\)이 된다. 예컨대 \(\text{df} = 19\)에서는 \(19/8.907 - 1 = 1.13\) 대 \(1 - 19/32.852 = 0.42\)이다. 카이제곱분포의 오른쪽 꼬리가 길어 분모의 작은 카이제곱 값이 큰 \(\sigma^2\) 값을 만들기 때문에 신뢰구간의 위쪽 꼬리가 더 멀리 뻗는다. \(\square\)

\(\sigma\)의 95% 신뢰구간은 \((\sqrt{9.25}, \sqrt{34.13}) = (3.04, 5.84)\)이다. 구간 전체가 6보다 아래이므로 95% 신뢰수준에서 자료는 \(\sigma < 6\)과 부합한다. 동등하게, 값 \(\sigma = 6\)(즉 \(\sigma^2 = 36\))은 \(\sigma^2\)의 95% 신뢰구간 \((9.25, 34.13)\) 밖에 있다. 따라서 5% 유의수준에서 \(\sigma < 6\)이라고 결론지을 충분한 증거가 있다. \(\square\)

이유. \(g(t)=\sqrt t\)는 \((0,\infty)\)에서 순증가하는 함수다. 따라서

로 두 사건이 논리적으로 동일하다. 확률도 같을 수밖에 없다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n, M = 20, 40_000
x = rng.normal(0, 2, (M, n))                  # σ = 2, σ² = 4
s2 = x.var(1, ddof=1)
lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)

print(f"σ² = 4 포함확률  {np.mean((lo <= 4) & (4 <= hi)):.6f}")
print(f"σ  = 2 포함확률  "
      f"{np.mean((np.sqrt(lo) <= 2) & (2 <= np.sqrt(hi))):.6f}")

σ² = 4 포함확률  0.948075
σ  = 2 포함확률  0.948075

소수 여섯째 자리까지 같다. 우연이 아니라 같은 사건이기 때문이다.

일반 원리. 모수 \(\theta\)의 구간에 단조함수 \(g\)를 적용하면 \(g(\theta)\)의 구간이 되고, 신뢰수준이 그대로 보존된다. 이것을 변환 등변성이라 한다.

주의할 점 세 가지.

1. 감소함수면 양끝이 뒤바뀐다. \(\theta\)의 구간 \((L,U)\)에 \(g(t)=1/t\)를 적용하면 \((1/U,\ 1/L)\)이다.

2. 점추정값은 옮겨 가지 않는다. \(E[S^2]=\sigma^2\)이지만 \(E[S]\ne\sigma\)다. 앞서 본 \(c_4\) 보정이 필요한 이유다. 구간은 옮겨 가고 불편성은 옮겨 가지 않는다.

3. 왈드 구간은 등변이 아니다. \(\hat\theta\pm z\widehat{\operatorname{SE}}\)에 \(g\)를 적용한 것과, \(g(\hat\theta)\pm z\widehat{\operatorname{SE}}_g\)는 다르다. 델타법으로 표준오차를 다시 계산하면 결과가 달라진다. 추축량이나 우도비로 만든 구간만이 등변이다.

폭은 \(S^2\)에 비례한다.

따라서 폭의 분포는 \(S^2\)의 분포, 즉 스케일된 카이제곱이다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
n, M = 20, 40_000
x = rng.normal(0, 2, (M, n))                  # σ² = 4
s2 = x.var(1, ddof=1)
lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)
w = hi - lo

q = np.percentile(w, [10, 25, 50, 75, 90])
print(f"폭의 십분위 {q[0]:.3f}  사분위 {q[1]:.3f}  중앙값 {q[2]:.3f}  "
      f"사분위 {q[3]:.3f}  구분위 {q[4]:.3f}")
print(f"평균 {w.mean():.3f}   90%/10% 비 {q[4] / q[0]:.2f}")

폭의 십분위 3.822  사분위 4.764  중앙값 6.015  사분위 7.448  구분위 8.927
평균 6.231   90%/10% 비 2.34

폭이 3.8에서 8.9까지 흔들린다. 열 번 중 한 번은 폭이 8.9를 넘고, 또 한 번은 3.8보다 좁다. 같은 \(n\), 같은 모집단인데 구간의 길이가 두 배 넘게 차이 난다.

왜 중요한가.

1. 표본크기 계획이 어렵다. "폭을 2 이하로"라는 목표를 세워도, \(S^2\)이 크게 나오면 달성하지 못한다. 확률적 보장("90% 확률로 폭이 2 이하")으로 계획을 세워야 하며, 이를 위한 공식이 따로 있다.

2. 좁은 구간이 좋은 소식이 아니다. 폭이 작다는 것은 \(S^2\)이 작았다는 뜻이고, 그것은 모분산을 과소추정했을 가능성이 높다는 뜻이다. 좁은 구간과 낮은 하한이 함께 나타난다.

3. 평균 폭만 보고하면 오해를 부른다. 방법 비교에서 "평균 폭 6.24"만 적으면 이 변동이 숨는다. 사분위나 백분위를 함께 적는다.

평균과의 대비. \(t\) 구간의 폭도 \(S\)에 비례하므로 확률변수지만, \(S\)는 \(S^2\)보다 변동이 작다(제곱근이 변동을 압축한다). 그래서 분산 구간의 폭이 평균 구간의 폭보다 훨씬 불안정하다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
M = 20_000
print(f"{'n':>6s} {'포함확률':>9s}")
for n in [20, 50, 100, 500, 2000]:
    y = rng.exponential(1.0, (M, n))          # 분산 1, 초과첨도 6
    s2 = y.var(1, ddof=1)
    lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
    hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)
    print(f"{n:6d} {np.mean((lo <= 1) & (1 <= hi)):9.4f}")

     n     포함확률
    20    0.7296
    50    0.6993
   100    0.6970
   500    0.6772
  2000    0.6747

회복되기는커녕 오히려 나빠진다. \(n=20\)에서 0.73, \(n=2000\)에서 0.67이다.

왜 그런가. 평균의 \(t\) 구간은 중심극한정리 덕에 \(n\)이 커지면 회복된다. 분산은 다르다.

로 점근분산이 4차 적률에 의존한다. 정규에서는 \(\mu_4-\sigma^4=2\sigma^4\)인데, 지수분포에서는 \(\mu_4=9\sigma^4\)이라 \(8\sigma^4\)로 네 배다.

한편 카이제곱 구간은 \((n-1)S^2/\sigma^2\sim\chi^2_{n-1}\)을 쓰는데, 이 분포의 분산이 \(2(n-1)\)이므로 \(2\sigma^4/n\)짜리 변동만 반영한다. 실제 변동이 \(8\sigma^4/n\)이므로 구간의 폭이 필요한 것의 절반이다.

핵심. \(n\to\infty\)에서 이 비율은 변하지 않는다. \(\sqrt{2/8}=0.5\)가 그대로 남는다. 카이제곱 구간은 비정규 아래에서 점근적으로도 타당하지 않다. 극한 포함확률을 직접 계산하면

from scipy import stats
import numpy as np

for kappa, name in [(0, "정규"), (6, "지수"), (3, "t(6)쯤")]:
    r = np.sqrt(2 / (kappa + 2))              # 반영된 폭 / 필요한 폭
    print(f"{name:6s} 초과첨도 {kappa}  극한 포함확률 "
          f"{2 * stats.norm.cdf(1.96 * r) - 1:.4f}")

정규     초과첨도 0  극한 포함확률 0.9500
지수     초과첨도 6  극한 포함확률 0.6729
t(6)쯤  초과첨도 3  극한 포함확률 0.7849

극한값 0.673이 모의실험의 0.675와 잘 맞는다.

교훈. "\(n\)이 크면 괜찮다"는 직관이 모든 절차에 적용되지 않는다. 중심극한정리가 보호하는 것은 평균이지 분산이 아니다. 분산·첨도·상관계수처럼 고차 적률에 의존하는 양은 분포가정을 따로 확인해야 한다.

import numpy as np
from scipy import stats

def bonett(x, alpha=0.05):
    n = len(x)
    z = stats.norm.ppf(1 - alpha / 2)
    trim = 1 / (2 * np.sqrt(max(n - 4, 1)))
    m = stats.trim_mean(x, min(trim, 0.49))
    xb = x.mean()
    g4 = n * np.sum((x - m)**4) / np.sum((x - xb)**2)**2
    se = np.sqrt((g4 - (n - 3) / n) / (n - 1))
    c = n / (n - z)
    s2 = x.var(ddof=1)
    return c * s2 * np.exp(-z * se), c * s2 * np.exp(z * se)

rng = np.random.default_rng(13)
n, M, B = 50, 3_000, 999
hc = hb = hp = 0
for _ in range(M):
    x = rng.exponential(1.0, n)               # 참 분산 1
    s2 = x.var(ddof=1)
    lo = (n - 1) * s2 / stats.chi2.ppf(0.975, n - 1)
    hi = (n - 1) * s2 / stats.chi2.ppf(0.025, n - 1)
    hc += (lo <= 1 <= hi)
    lb, hb_ = bonett(x)
    hb += (lb <= 1 <= hb_)
    idx = rng.integers(0, n, (B, n))
    star = x[idx].var(1, ddof=1)
    l, h = np.percentile(star, [2.5, 97.5])
    hp += (l <= 1 <= h)
print(f"카이제곱 {hc / M:.4f}   보넷 {hb / M:.4f}   "
      f"부트스트랩 백분위 {hp / M:.4f}")

카이제곱 0.7113   보넷 0.8813   부트스트랩 백분위 0.8073

순위. 보넷(0.881) > 부트스트랩 백분위(0.807) > 카이제곱(0.711).

부트스트랩이 기대만 못한 이유. 백분위 부트스트랩은 일차정확(\(O(n^{-1/2})\) 오차)이라 치우친 통계량에서 수렴이 느리다. \(S^2\)의 표집분포는 지수분포 아래에서 크게 오른쪽으로 치우쳐 있고, 백분위법은 그 치우침을 반대 방향으로 반영하는 결함이 있다.

개선책. BCa나 부트스트랩-\(t\)를 쓰면 이차정확(\(O(n^{-1})\))이 되어 훨씬 낫다. 가속상수 \(a\)가 왜도를 보정한다.

권고.

더 근본적인 질문. 지수분포 자료에서 "분산"이 정말 관심사인가? 지수분포는 \(\sigma=\mu\)이므로 평균만 추정하면 분산이 따라온다. 분포를 모형화하는 것이 요약통계를 강건하게 추정하는 것보다 나은 경우가 많다.

분포 목록. 첨도를 넓게 훑되 참 분산을 아는 것들로 고른다.

import numpy as np
from scipy import stats

rows = [
    ("균등(0,1)",        1 / 12,           -1.2),
    ("정규(0,1)",        1.0,               0.0),
    ("로지스틱(0,1)",     np.pi**2 / 3,      1.2),
    ("t(9) / √(9/7)",   1.0,               1.2),
    ("지수(1)",          1.0,               6.0),
    ("t(5) / √(5/3)",   1.0,               6.0),
    ("감마(2,1)",        2.0,               3.0),
    ("로그정규(0,1)",     np.exp(2) - np.exp(1), 110.9),
    ("오염 0.9N(0,1)+0.1N(0,9)", 0.9 + 0.1 * 9, 5.3),
]
print(f"{'분포':26s} {'참 분산':>10s} {'초과첨도':>9s}")
for name, var, kur in rows:
    print(f"{name:26s} {var:10.4f} {kur:9.1f}")

분포                               참 분산      초과첨도
균등(0,1)                        0.0833      -1.2
정규(0,1)                        1.0000       0.0
로지스틱(0,1)                      3.2899       1.2
t(9) / √(9/7)                  1.0000       1.2
지수(1)                          1.0000       6.0
t(5) / √(5/3)                  1.0000       6.0
감마(2,1)                        2.0000       3.0
로그정규(0,1)                      4.6708     110.9
오염 0.9N(0,1)+0.1N(0,9)         1.8000       5.3

설계상 주의점.

1. 분산을 1로 정규화한다. 분포마다 참 분산이 다르면 폭을 비교할 수 없다. \(X/\sigma\)로 표준화하면 포함확률은 그대로, 폭은 비교 가능해진다.

2. 첨도가 지표다. 위 표에서 보듯 초과첨도가 같으면 포함확률도 거의 같다. 지수와 \(t_5\)가 둘 다 6이고, 실제로 카이제곱 구간의 포함확률이 비슷하다. 왜도는 상대적으로 덜 중요하다.

3. 첨도가 존재하지 않는 분포도 넣는다. \(t_4\) 이하는 4차 적률이 무한대라 점근이론이 아예 성립하지 않는다. 이런 극단도 한 번 확인해 둘 가치가 있다.

4. 로그정규는 극단적이다. 초과첨도 110.9로, 어떤 방법도 살아남지 못한다. "이 방법은 로그정규에서도 잘 된다"는 주장은 의심해야 한다.

5. 오염 모형을 반드시 넣는다. 실무의 비정규성은 매끄러운 두꺼운 꼬리보다 소수의 이상치 형태인 경우가 많다. 혼합분포가 그것을 흉내 낸다.

6. \(n\)의 격자. 분산은 \(n\)에 민감하므로 \(n\in\{10,20,50,100,500\}\) 정도를 훑는다. 앞 문제에서 보았듯 큰 \(n\)에서 오히려 나빠질 수 있어, 큰 \(n\)을 빼면 잘못된 안도감을 준다.

7. 평가 지표. 포함확률, 평균 폭, 좌우 누락의 균형(구간이 참값 아래인지 위인지). 치우친 분포에서는 누락이 한쪽으로 쏠린다.

8. 참 분산을 코드에서 계산하지 말고 손으로 확인한다. 표본에서 추정한 값을 "참값"으로 쓰면 모의실험이 무의미해진다. 흔한 버그다.