S²의 표준오차¶

(1) 이론값. 위에서 한 계산을 다시 모아 적는다. \(\text{Uniform}(0,1)\)은 \(\sigma^2 = 1/12\), \(\mu_4 = 1/80\)이므로 불편성에서

이고, 일반 공식에 \(n = 5\)를 넣으면 \((n-3)/(n-1) = 2/4 = 1/2\)이라

를 얻는다. 모집단이 정규라고 가정하면 \(\mu_4 = 3\sigma^4\)이므로

이 되어 \(39\%\) 더 크다. 분산으로는 배율이

다. 균등모집단의 꼬리가 가벼워 \(S^2\)이 정규모집단보다 덜 흔들리는 것이고, 정규를 가정하는 쪽이 보수적(과대평가)이다.

크기 자체도 짚어 둘 만하다. \(\operatorname{SE}/E[S^2] = 0.042492/0.083333 = 0.51\)이다. 표준오차가 추정하려는 값의 절반을 넘는다.

치우침에 대해서는 미리 어림할 것이 하나 있다. \(S^2 \ge 0\)이고 \(\operatorname{SE}\)가 평균의 절반이 넘으므로 \(\mu - 2\operatorname{SE} = 0.0833 - 0.0850 = -0.0017\)이 음수다. 곧 \(\mu - 2\operatorname{SE}\) 아래의 확률은 정확히 \(0\)이고, 그 거리에서는 치우침이 극단적으로 드러나야 한다. \(\mu \pm \operatorname{SE}\)에서는 어떨지 모의실험으로 재어 본다.

(2) 모의실험.

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

np.random.seed(0)

# 크기 5짜리 균등표본을 1만 번 뽑아 그때마다 S^2 을 기록한다.
S_square = []
for _ in range(10_000):
    x = np.random.uniform(size=(5,))
    sigma = x.std(ddof=1)     # ddof=1 이라야 표본표준편차다
    S_square.append(sigma ** 2)

# 1만 개 S^2 값의 **평균**과 **표준편차**를 낸다.
#   평균  -> S^2 의 중심. 참 분산 1/12 ≈ 0.0833 에 가까워야 한다(불편성).
#   표준편차 -> 그것이 곧 S^2 의 **표준오차**다.
# 표준오차란 "통계량의 표집분포의 표준편차"이므로,
# 통계량 값을 잔뜩 모아 그 표준편차를 재면 그것이 표준오차다.
average = np.array(S_square).mean()
standard_error = np.array(S_square).std()

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

# 이론값. Uniform(0,1) 은 sigma^2 = 1/12, mu_4 = 1/80 이다.
n, B = 5, 10_000
sigma2, mu4 = 1 / 12, 1 / 80
var_s2 = (mu4 - (n - 3) / (n - 1) * sigma2 ** 2) / n
se_s2 = np.sqrt(var_s2)
print(f"이론  E[S^2] = {sigma2:.6f},  Var(S^2) = {var_s2:.8f},  SE(S^2) = {se_s2:.6f}")
print(f"몬테카를로 오차  평균 {se_s2 / np.sqrt(B):.6f},  표준오차 {se_s2 / np.sqrt(2 * B):.6f}")
print(f"어긋남을 오차 단위로  평균 {abs(average - sigma2) / (se_s2 / np.sqrt(B)):.2f}배,  "
      f"표준오차 {abs(standard_error - se_s2) / (se_s2 / np.sqrt(2 * B)):.2f}배")
print(f"정규모집단이라면 SE = sigma^2 sqrt(2/(n-1)) = {sigma2 * np.sqrt(2 / (n - 1)):.6f},  "
      f"분산 배율 = {(n - 1) / (2 * n) * (mu4 / sigma2 ** 2 - (n - 3) / (n - 1)):.4f}")
print(f"SE / E[S^2] = {se_s2 / sigma2:.4f}")

# 치우침을 수로 읽는다. 1 표준오차 거리와 2 표준오차 거리를 갈라서 본다.
s = np.array(S_square)
print(f"왜도(S^2) 모의 = {stats.skew(s):.4f},  중앙값 {np.median(s):.5f} < 평균 {sigma2:.5f}")
print(f"mu +- SE  = [{sigma2 - se_s2:+.5f}, {sigma2 + se_s2:+.5f}]   바깥 확률  왼쪽 {(s < sigma2 - se_s2).mean():.4f}  오른쪽 {(s > sigma2 + se_s2).mean():.4f}")
print(f"mu +- 2SE = [{sigma2 - 2 * se_s2:+.5f}, {sigma2 + 2 * se_s2:+.5f}]   바깥 확률  왼쪽 {(s < sigma2 - 2 * se_s2).mean():.4f}  오른쪽 {(s > sigma2 + 2 * se_s2).mean():.4f}")

# 히스토그램에 이론값을 겹쳐 그린다.
fig, ax = plt.subplots(figsize=(12, 3))
ax.set_title("Sampling Distribution of S^2")
ax.hist(S_square, bins=100, density=True, alpha=0.3)
# 평균과 평균 ± 1 표준오차를 세로선으로 표시한다.
# S^2 의 분포는 오른쪽으로 치우쳐 있어 이 구간이 대칭이 아님에 주의하라.
ax.vlines(average, ymin=0, ymax=12, color="k", lw=5, label="Mean")
ax.vlines(average + standard_error, ymin=0, ymax=12,
          color="k", ls="--", label="Mean +/- SE")
ax.vlines(average - standard_error, ymin=0, ymax=12,
          color="k", ls="--")
ax.legend()
plt.show()

출력:

Estimated Mean of S^2:   0.0838
Standard Error of S^2:   0.0425
이론  E[S^2] = 0.083333,  Var(S^2) = 0.00180556,  SE(S^2) = 0.042492
몬테카를로 오차  평균 0.000425,  표준오차 0.000300
어긋남을 오차 단위로  평균 1.18배,  표준오차 0.04배
정규모집단이라면 SE = sigma^2 sqrt(2/(n-1)) = 0.058926,  분산 배율 = 0.5200
SE / E[S^2] = 0.5099
왜도(S^2) 모의 = 0.4122,  중앙값 0.08072 < 평균 0.08333
mu +- SE  = [+0.04084, +0.12583]   바깥 확률  왼쪽 0.1735  오른쪽 0.1718
mu +- 2SE = [-0.00165, +0.16832]   바깥 확률  왼쪽 0.0000  오른쪽 0.0304

두 숫자가 맞는다. 모의 평균 \(0.0838\)이 이론 \(0.083333\)에서 몬테카를로 오차 \(0.000425\)의 \(1.18\)배, 모의 표준오차 \(0.0425\)가 이론 \(0.042492\)에서 오차 \(0.000300\)의 \(0.04\)배 떨어져 있다. 표준오차 쪽은 소수 넷째 자리까지 그대로 맞았다. \(\mu_4\)가 들어간 일반 공식이 옳다는 확인이고, 정규를 가정해 얻은 \(0.058926\)을 썼다면 \(39\%\) 과대평가였을 것이다.

치우침은 \(1\) 표준오차 거리에서 보이지 않는다. 왜도가 \(+0.41\)로 양수이고 중앙값 \(0.0807\)이 평균 \(0.0833\)보다 작으니 분포가 오른쪽으로 끌린 것은 분명하다. 그런데 \(\mu \pm \operatorname{SE}\) 바깥의 확률은 왼쪽 \(0.1735\), 오른쪽 \(0.1718\)로 사실상 같다. 비율의 몬테카를로 오차가 \(\sqrt{0.17 \times 0.83/10^4} = 0.0038\)이므로 두 값의 차 \(0.0017\)은 오차의 절반도 안 된다.

이것이 우연이 아님을 확인하려고 되풀이를 1천만 번으로 늘려 보았다. 왼쪽 \(0.17686\), 오른쪽 \(0.16743\)이 되어 오차 \(0.00012\) 기준으로는 분명히 다르지만, 큰 쪽이 왼쪽이다. 왜도가 \(+0.41\)인 감마분포에 같은 계산을 하면 \(0.1572\)와 \(0.1567\)로 거의 같은 값이 나오니, 한 표준오차 거리에서는 왜도가 좌우 확률을 거의 가르지 않는다는 것이 일반적인 사정이다. 왜도는 3차 적률이어서 훨씬 먼 꼬리에서 결정된다.

\(2\) 표준오차 거리에서는 완전히 드러난다. \(\mu - 2\operatorname{SE} = -0.00165\)가 음수이므로 그 아래의 확률이 정확히 \(0\)이고, 오른쪽은 \(0.0304\)다. 정규분포라면 양쪽이 각각 \(0.0228\)이어야 한다. 한쪽이 \(0\)이고 다른 쪽이 \(3\%\)라는 것은 대칭으로는 설명할 길이 없는 모습이다.

그래서 그림의 점선 두 개를 "전형적인 범위"로 읽을 때 조심해야 한다. 폭만 보면 \(0.0425\)가 평균 \(0.0833\)의 \(51\%\)라 다섯 개짜리 표본의 \(S^2\) 하나로 모분산을 논하기 어렵다는 것이 바로 읽힌다. 그러나 그 구간이 비대칭한 분포를 대칭한 두 선으로 가둔 것이라는 사실은 점선만 보아서는 알 수 없고, 위의 마지막 두 줄처럼 양쪽을 갈라 세어야 드러난다.

\(N(0, 1)\)에서 \(\sigma^2 = 1\)이다.

\(\square\)

\(Q = (n-1)S^2/\sigma^2\)이라 하자. 그러면 \(Q \sim \chi^2(n-1)\)이고 \(\text{Var}(Q) = 2(n-1)\)이다.

\(S^2 = \sigma^2 Q / (n-1)\)이므로:

따라서:

\(\square\)

\(\mu = 1/2\)인 \(X \sim \text{Uniform}(0, 1)\)에 대해:

\(u = x - 1/2\)로 치환하면 \(du = dx\)이고 적분 범위는 \(-1/2\)에서 \(1/2\)이 된다:

\(\square\)

표본분산 \(S^2\)은 아래로 0에서 유계이지만 (원리상) 위로는 유한한 한계가 없다. \(n\)이 작으면 제약 \(S^2 \ge 0\)이 왼쪽 꼬리를 잘라 내는 "바닥"을 만드는 반면, 이따금 나타나는 극단적인 관측값이 \(S^2\)을 큰 값으로 밀어 올려 긴 오른쪽 꼬리를 만든다.

\(n\)이 커지면 두 가지 효과가 분포를 대칭화한다:

1. \(S^2\)에 대한 중심극한정리: \(n\)이 크면 \(S^2\)은 약하게 의존하는 많은 항(제곱편차)의 합에 가까우므로, 중심극한정리 유형의 논증에 의해 그 분포가 정규분포로 수렴한다.
2. 카이제곱의 수렴: 정규모집단에서 \((n-1)S^2/\sigma^2 \sim \chi^2(n-1)\)이다. \(\chi^2(k)\)의 왜도는 \(2\sqrt{2/k}\)이며 \(k = n - 1\)이 커지면 0으로 간다.

두 효과 모두 \(n\)이 클 때 \(S^2\)의 분포가 근사적으로 대칭(정규)이 되고 표준오차가 퍼짐을 잘 요약해 줌을 뜻한다. \(\square\)

import numpy as np
np.random.seed(0)

# 크기 5짜리 지수분포 표본에서 S^2 을 1만 번 계산한다.
# ddof=1 이 표본분산(n-1로 나눔)을 뜻한다.
S_square = [np.random.exponential(size=5).var(ddof=1) for _ in range(10_000)]
empirical_se = np.std(S_square)

# 이론값 Var(S^2) = (1/n)(mu4 - (n-3)/(n-1) * sigma^4)
# Exp(1)에서 sigma^2 = 1, mu4 = 9 이므로 (1/5)(9 - 2/4) = 1.7
theoretical_se = np.sqrt(1.7)
print(f"경험적 SE = {empirical_se:.4f}")
print(f"이론적 SE = {theoretical_se:.4f}")

출력:

경험적 SE = 1.3237
이론적 SE = 1.3038

n = 5로 아주 작은데도 두 값이 잘 맞는다. 다만 4차 적률이 들어가는 공식이라 수렴이 느리므로, 표본이 작으면 이 정도(1.5%) 차이는 정상이다.

\(\sigma^2 = 1\), \(\mu_4 = 9\), \(n = 5\)로 이론값을 계산하면:

비교하자면 정규이론 표준오차는 \(\sigma^2 \sqrt{2/(n-1)} = \sqrt{2/4} = \sqrt{0.5} \approx 0.707\)이다.

지수 모집단의 표준오차가 정규이론 값보다 약 1.84배 큰데, 이는 지수분포의 두꺼운 꼬리(초과첨도 \(= 6\))를 반영한다. 모의실험의 경험적 표준오차는 1.304에 가깝게 나온다. \(\square\)

정규모집단. \(\mu_4 = 3\sigma^4\)을 넣으면

으로 카이제곱 결과와 정확히 일치한다. ✓

초과첨도 형태. \(\gamma_2 = \mu_4/\sigma^4 - 3\)으로 두면 \(\mu_4 = (\gamma_2+3)\sigma^4\)이므로

이다(마지막은 \(n\)이 클 때).

읽는 법. \(S^2\)의 표준오차를 정하는 것은 분산이 아니라 4차 적률, 곧 첨도다.

지수 모집단에서 \(S^2\)의 분산이 정규 이론값의 네 배다. 표준오차로는 두 배다. 이것이 앞 연습문제에서 본 1.84배(\(n=5\)라 \(2n/(n-1)\) 보정이 들어간 값)의 정체다.

실무적 함의가 크다. 평균에 대한 추론은 중심극한정리가 구해 주지만 분산에 대한 추론은 그렇지 않다. 첨도는 자체로 추정하기 어려운 양이라(그 추정값의 정밀도가 8차 적률에 달려 있다), 분산의 신뢰구간을 자료에서 안전하게 만드는 일은 생각보다 까다롭다. 부트스트랩이 유용한 대목이다.

언제나 과소추정하며, \(n\)이 커지면 사라진다.

\(S^2\)은 불편인데 \(S\)는 아닌 이유. 제곱근이 오목함수이기 때문이다. 옌센 부등식에 따라

이다. 등호는 \(S^2\)이 상수일 때만 성립하는데, 그것은 \(n=\infty\)인 경우다.

일반 원리. 불편성은 비선형 변환에서 보존되지 않는다. \(E[\hat\theta]=\theta\)라 해도 \(E[g(\hat\theta)] \ne g(\theta)\)이며, \(g\)가 오목이면 아래로, 볼록이면 위로 치우친다. 앞서 본 \(\hat\lambda = 1/\bar X\)가 볼록 변환이라 과대추정하는 것과 정확히 반대 방향이다.

보정. \(S/c_4\)를 쓰면 불편이 되고, 공정능력지수나 관리도에서 실제로 이 보정을 쓴다. 다만 보정한 추정량이 평균제곱오차 기준으로는 오히려 나쁠 수 있다는 점을 기억할 일이다. 편향을 없애면서 분산이 늘기 때문이다.

실무에서 \(n \ge 25\)면 편향이 1% 미만이라 무시해도 좋다. 문제가 되는 것은 \(n\)이 아주 작을 때, 특히 여러 개의 작은 표본에서 얻은 \(s\)를 평균 내는 경우다. 편향이 상쇄되지 않고 그대로 쌓인다.

\(g(x) = \sqrt x\)에 대해 \(g'(\sigma^2) = 1/(2\sigma)\)이므로

이고

이다. \(\square\)

필요 표본크기. 상대 표준오차가 0.10 이하이려면

이다.

평균과의 대비가 인상적이다. 평균을 상대오차 10%로 추정하려면 \(\sigma/(\sqrt n|\mu|) \le 0.1\), 즉 변동계수가 0.5라면 \(n \ge 25\)면 된다. 그런데 표준편차 자체를 같은 정밀도로 알려면 51개가 필요하고, 분산을 10% 이내로 알려면 \(\sqrt{2/(n-1)} \le 0.1\)에서 \(n \ge 201\)이 필요하다.

산포를 정밀하게 아는 것은 중심을 아는 것보다 훨씬 비싸다. 표본크기 계산이 보통 평균을 기준으로 이루어지므로, 같은 연구에서 분산에 관한 결론은 훨씬 불확실하다는 점을 잊지 말아야 한다.

카이제곱분포의 왜도가 \(\sqrt{8/k}\)이고 선형변환은 왜도를 바꾸지 않으므로, \(k=n-1\)에서

이다.

정규모집단에서 \(\bar X\)는 왜도가 정확히 0인데 \(S^2\)은 \(n=100\)에서도 0.284나 남아 있다.

이 차이가 실무에서 갖는 뜻은 분명하다.

• \(\bar X\)의 신뢰구간은 정규모집단에서 정확하고, 비정규모집단에서도 중심극한정리가 빠르게 구해 준다.
• \(S^2\)의 신뢰구간은 대칭으로 만들면 안 된다. \(\bar s^2 \pm z\cdot\operatorname{SE}\) 꼴은 \(n\)이 꽤 커도 포함확률이 어긋나며, 특히 아래쪽 한계가 음수가 되기도 한다. 앞서 본 카이제곱 기반 구간이 비대칭인 이유가 이것이다.
• 분산에 대한 검정에서 단측·양측의 임계값을 대칭으로 잡으면 틀린다.

한마디로 \(S^2\)은 \(\bar X\)보다 정규성에 훨씬 늦게 다가간다. 분산을 다룰 때는 \(n\)이 크다는 이유만으로 정규근사를 믿으면 안 된다.

방법 1 — 정규이론 대입.

\((n-1)S^2/\sigma^2\sim\chi^2_{n-1}\)을 쓰고 미지의 \(\sigma^2\)을 \(s^2\)으로 대신한 것이다.

• 가정: 모집단이 정규여야 한다. 이것이 강한 가정이다. 연습문제 6에서 보았듯 초과첨도가 6인 모집단에서는 참 표준오차가 이 값의 두 배다.
• 장점: 계산이 한 줄이고 \(n\)이 작아도 쓸 수 있다.

방법 2 — 적률 대입.

연습문제 6의 일반 공식에 표본 4차 적률을 넣은 것이다.

• 가정: 정규성을 요구하지 않는다. 다만 모집단의 8차 적률이 유한해야 \(\hat\mu_4\) 자체가 안정적이다.
• 단점: \(\hat\mu_4\)의 추정오차가 매우 크다. 4제곱이 들어가므로 관측값 하나가 결과를 좌우할 수 있고, \(n\)이 작으면 쓸 수 없다.

방법 3 — 부트스트랩. \(S^2\)을 재표본마다 계산해 그 표준편차를 쓴다.

• 가정: 분포 가정이 거의 없다. 관측이 독립이고 \(n\)이 너무 작지 않으면 된다.
• 장점: 자동으로 첨도의 영향을 반영한다. 세 방법 가운데 대체로 가장 믿을 만하지만, 꼬리가 두꺼운 모집단에서는 \(n = 30\)쯤에서도 포함확률이 명목값에 한참 못 미친다(지수모집단 쪽 연습문제 10이 이를 재어 본다).

권고. 정규성을 확신할 수 있으면 방법 1이 가장 효율적이다. 그렇지 않으면 부트스트랩이 안전하다. 방법 2는 이론적으로 옳지만 실무에서 \(\hat\mu_4\)의 불안정성 때문에 잘 쓰이지 않는다.