분산추정량의 평균제곱오차¶

(a) \(W = \sum(X_i - \bar X)^2 \sim \sigma^2 \chi^2_{n-1}\)일 때: \(\mathbb{E}[\hat\sigma^2_c] = (n-1)\sigma^2/c\), \(\mathrm{Var}(\hat\sigma^2_c) = 2(n-1)\sigma^4/c^2\).

\(\mathrm{MSE}(c) = (\sigma^4/c^2)[(n-1-c)^2 + 2(n-1)]\).

(b) 미분하여 0으로 놓으면 \(c^* = n + 1\).

(c) \(n = 10\), \(\sigma = 1\): \(\mathrm{MSE}(9) \approx 0.222\), \(\mathrm{MSE}(10) = 0.190\), \(\mathrm{MSE}(11) \approx 0.182\). \(c^* = 11\)이 이긴다.

최소 평균제곱오차: 분모 \(n+1\)이 MLE와 불편추정량을 모두 이긴다. 맞바꿈은 이렇다: 작은 편향을 받아들여 더 큰 분산 감소를 얻는다.

평균제곱오차에서 이기는데도 \(n+1\) 분모가 거의 쓰이지 않는 이유는 \(S^2\)의 해석이 더 깔끔하고(불편) 적당한 \(n\)에서는 평균제곱오차 차이가 작기 때문이다.

MSE는 제곱오차 손실 \(L(\hat\theta, \theta) = (\hat\theta - \theta)^2\)이다. 널리 쓰이는 이유:

• 수학적으로 다루기 쉽다(기댓값의 선형성, 편향 + 분산으로의 분해).
• 어디서나 미분 가능하다 — 해석적 최적화가 가능하다.
• 큰 오차에 큰 벌점을 준다(제곱).

대안:

• 절대오차: \(|\hat\theta - \theta|\). 사후 중앙값에서 최소가 되며 이상점에 더 로버스트하다.
• 분위 손실: \(\rho_\tau(\hat\theta - \theta) = (\hat\theta - \theta)(\tau - \mathbf 1\{\hat\theta < \theta\})\). 분위회귀에서 쓰인다.
• 0-1 손실: \(\mathbf 1\{\hat\theta \ne \theta\}\). 분류에서 쓰인다.
• 로그 손실 / KL 발산: 밀도추정, 확률적 예측에서 쓰인다.

손실의 선택은 어떤 종류의 오차에 벌점을 줄지를 반영한다. MSE는 기본값이지만 언제나 옳은 선택은 아니다.

(관측값 하나로) \(\mathbf X \sim N(\boldsymbol\mu, I_p)\)에서 \(\boldsymbol\mu \in \mathbb{R}^p\)를 추정한다. MLE: \(\hat{\boldsymbol\mu} = \mathbf X\). 평균제곱오차: \(\mathbb{E}[\|\mathbf X - \boldsymbol\mu\|^2] = p\).

James-Stein 추정량: \(\hat{\boldsymbol\mu}_{\text{JS}} = (1 - (p-2)/\|\mathbf X\|^2) \mathbf X\) (0 쪽으로 축소).

Stein(1956)이 보인 것: \(p \ge 3\)일 때 모든 \(\boldsymbol\mu\)에 대해 \(\mathrm{MSE}(\hat{\boldsymbol\mu}_{\text{JS}}) < p\)이다. MLE가 지배당한다 — 언제나 더 나은 추정량이 있다.

직관: 임의의 점 쪽으로 축소하더라도 JS 추정량이 평균제곱오차에서 MLE보다 낫다. 고차원에서는 도입된 편향보다 분산 감소가 더 크다.

현대 통계학: 고차원 문제에서 축소추정량(능형, 라소, James-Stein)이 평균제곱오차 기준으로 MLE를 일상적으로 이긴다. 경험적 Bayes와 계층 모형이 이를 형식화한다.

\(\mathrm{MSE}(S^2) = 2\sigma^4/(n-1)\). 제곱근을 취하면 \(\sqrt{\mathrm{MSE}} = \sigma^2 \sqrt{2/(n-1)}\).

\(n' = 4n\)이면 \(\sqrt{\mathrm{MSE}'} = \sigma^2 \sqrt{2/(4n - 1)} \approx \sigma^2 \sqrt{2/(n-1)}/2 = \sqrt{\mathrm{MSE}}/2\).

\(n\)이 네 배가 되면 절반이 된다. \(\bar X\)에서와 같은 \(\sqrt n\) 축척이며, 여기서는 분산추정의 RMSE에 적용된 것이다.

세 가지 모두:

• \(\mathrm{MSE}(\hat\sigma^2_{\text{MLE}})/\sigma^4 = (2n-1)/n^2 \to 2/n\).
• \(\mathrm{MSE}(S^2)/\sigma^4 = 2/(n-1) \to 2/n\).
• \(\mathrm{MSE}(\hat\sigma^2_{c^*})/\sigma^4 = 2/(n+1) \to 2/n\).

점근적으로 동등하다 — \(1/n\)의 1차 항까지 보면 셋 모두 MSE \(\approx 2\sigma^4/n\)이다.

함의: 분모의 선택은 유한표본에서만 중요하다. 큰 \(n\)에서는 편한 것을 쓰면 된다. 차이가 무시할 수 없는 작은 \(n\)에서는 손실함수에 따라 고르라: 불편성(\(n-1\)), MLE의 편의(\(n\)), 최소 평균제곱오차(\(n+1\)).

평균제곱오차.

여기서 \(V = \operatorname{Var}(S^2) = 2\sigma^4/(n-1)\)이고 \(\Delta = \sigma^2-\sigma_0^2\)이다.

최적화. \(w\)로 미분해 0으로 두면

읽는 법.

• \(\Delta=0\)(목표가 정확히 맞음): \(w^*=0\)으로 자료를 전혀 쓰지 않는다. 목표값이 참값이면 그것을 쓰는 것이 당연하다.
• \(\Delta\)가 크면 \(w^*\to1\)로 표본분산을 그대로 쓴다.
• \(V\)가 크면(\(n\)이 작으면) \(w^*\)가 작아져 더 많이 축소한다.

최소 MSE.

언제나 \(S^2\)보다 낫다(\(V\)가 \(S^2\)의 MSE이므로). 다만 \(w^*\)가 미지의 \(\Delta\)에 의존해 직접 쓸 수 없다.

실무의 해법.

• \(\Delta\)를 추정한다. \(\hat\Delta^2 = \max(0,\ (S^2-\sigma_0^2)^2-V)\) 같은 추정값을 넣는다. 제임스-스타인의 발상과 같다.
• 경험 베이즈. 여러 분산을 동시에 추정할 때 그들의 흩어짐에서 \(\Delta\)를 추정한다. 유전자 발현 분석의 limma가 정확히 이 방식으로 유전자별 분산을 전체 평균 쪽으로 축소한다.
• 르두아-울프. 공분산행렬에서 표본공분산을 구조화된 목표(대각행렬이나 등상관 행렬)로 축소하며, 최적 \(w\)를 자료에서 해석적으로 추정한다. 자산 수가 관측 수에 가까울 때 필수적이다.

문제의 성격. 평균제곱오차는 척도에 의존한다. \(\sigma^2\)의 MSE를 최소로 하는 추정량이 \(\ln\sigma^2\)의 MSE도 최소로 하리라는 보장이 없다.

\(\ln\sigma^2\)의 경우. \(\hat\sigma^2_c = \text{SS}/c\)에 대해 \(\ln\hat\sigma^2_c = \ln\text{SS}-\ln c\)이므로, \(c\)를 바꾸는 것은 상수를 더하는 것이다. 따라서

• 분산은 \(c\)와 무관하다.
• 편향만 \(c\)가 결정한다. \(E[\ln\text{SS}] = \ln(2\sigma^2)+\psi\{(n-1)/2\}\)이므로, 불편이 되려면

$$ \ln c = \ln 2+\psi!\left(\frac{n-1}{2}\right) $$

즉 \(c^* = 2\exp\{\psi((n-1)/2)\}\)다.

\(n=10\)이면 \(\psi(4.5)=1.3888\)이므로 \(c^* = 2e^{1.3888} = 8.01\)로, \(n-1=9\)보다도 작다.

\(1/\sigma^2\)의 경우. \(E[1/\text{SS}]\)를 역카이제곱 적률로 계산하면

이므로 불편이 되려면 \(c = n-3\)이다. MSE를 최소로 하려면 더 작은 값이 나온다.

정리.

교훈. "최적 분모"는 무엇을 추정하느냐에 따라 달라진다. 이것이 불편성과 MSE 최적성이 모두 모수화에 의존한다는 사실의 구체적인 예다. \(\ln\sigma^2\) 척도에서 분산이 \(c\)에 무관하다는 점이 흥미로운데, 이는 로그 척도가 분산안정화 변환이기 때문이다.

상황. \(k\)개 집단에서 각각 \(S_i^2\)을 얻었다. 참 분산 \(\sigma_i^2\)들이 서로 다를 수 있다.

세 가지 선택.

축소가 이득인 조건. 앞 연습문제의 결과에서 \(w_i^*\)가 \(\Delta_i^2/(V_i+\Delta_i^2)\)이므로

• \(n_i\)가 작으면(\(V_i\)가 크면) 더 많이 축소한다.
• 집단 간 분산의 차이가 작으면 더 많이 축소한다.

제임스-스타인의 적용. 로그 척도에서 하면 깔끔하다. \(\ln S_i^2\)이 근사적으로 정규이고 분산이 \(2/(n_i-1)\)로 \(\sigma_i\)에 무관하므로, 정규 평균의 축소 문제로 환원된다. 따라서 \(k\ge3\)이면 제임스-스타인류의 축소가 모든 \(\{\sigma_i^2\}\)에서 총 손실을 줄인다.

실제로 이것이 분산의 계층모형이며

를 두고 경험 베이즈로 \((\mu_0,\tau^2)\)을 추정한 뒤 각 \(\sigma_i^2\)을 축소한다.

대표적 응용. 유전자 발현 분석이다. 유전자가 2만 개인데 표본이 6개뿐이면 유전자별 \(S_i^2\)이 극도로 불안정하다. 우연히 작게 나온 분산이 \(t\) 통계량을 폭발시켜 거짓 발견을 양산한다. limma의 조절 \(t\) 통계량은 분산을 전체 쪽으로 축소해 이 문제를 해결하며, 실질적으로 자유도를 빌려 온다.

주의. 축소는 개별 \(\sigma_i^2\)의 추정을 나쁘게 할 수 있다. 총 손실이 줄어드는 것이지 모든 집단이 개선되지는 않는다. 특정 집단의 분산 자체가 관심사라면 신중해야 한다.

공정한 비교의 전제. 네 모집단의 참 분산을 같게 맞춰야 MSE를 직접 비교할 수 있다. 모두 \(\sigma^2=1\)로 표준화한다.

비교할 추정량. \(\text{SS}/c\) 계열(\(c=n-1,n,n+1\))과 강건 대안(MAD 기반, IQR 기반, 절사분산).

설계.

rng = np.random.default_rng(0)
B = 50_000
for n in [5, 10, 30, 100]:
    for name, draw in populations.items():
        ss = np.empty(B)
        for b in range(B):
            x = draw(n)
            ss[b] = ((x - x.mean()) ** 2).sum()
        for c in [n - 1, n, n + 1]:
            est = ss / c
            print(n, name, c, est.mean() - 1, est.var(), np.mean((est - 1) ** 2))

볼 것.

1. 편향과 MSE를 따로. \(n-1\)이 모든 모집단에서 불편인지 확인한다(그래야 한다. 불편성은 분포에 무관하다).
2. 최적 \(c\)가 어떻게 변하는지. \(c\)를 연속적으로 바꿔 MSE 곡선을 그려 최소점을 찾는다. 이론값 \(c^*\approx n+1+n\gamma_2/2\)와 맞는지 확인한다.
3. 강건 추정량과의 비교. \(t_5\)에서 MAD 기반 추정량이 \(S^2\)보다 나은지 본다.
4. \(n\)에 따른 수렴. \(n\)이 커지면 세 분모의 차이가 사라지는지 확인한다.

주의할 점.

• \(B\)를 충분히. MSE 추정값의 상대 표준오차가 \(\sqrt{2/B}\) 수준이므로, 5% 차이를 구별하려면 \(B\ge10^4\)이 필요하다.
• 공통난수. 같은 표본에서 여러 \(c\)를 계산하면 차이의 분산이 크게 줄어 훨씬 적은 \(B\)로 구별할 수 있다. 위 코드가 그렇게 되어 있다.
• \(t_5\)의 4차 적률. \(t_5\)는 4차 적률이 유한하지만(\(\nu>4\)) 아슬아슬하다. MSE 추정값이 잘 수렴하는지 확인해야 하며, \(t_3\)을 넣으면 4차 적률이 없어 MSE 자체가 무한임을 기억한다.