흔한 분포에 대한 적률법¶

(a) 표본 적률: \(\bar x = 2.33\), \(m_2 \approx 2.77\).

\(\mathbb{E}[X] = \alpha\beta\), \(\mathrm{Var}(X) = \alpha\beta^2\)인 Gamma\((\alpha, \beta)\)에 대해:

\(\hat\beta = m_2/\bar x \approx 1.19\), \(\hat\alpha = \bar x^2/m_2 \approx 1.96\).

(b) scipy.stats.gamma.fit(x, floc=0)으로 구한 MLE: 전형적인 출력은 \(\hat\alpha \approx 1.69\), \(\hat\beta \approx 1.38\).

(c) 히스토그램에 두 밀도를 겹쳐 그린다. \(n = 10\)에서는 두 추정량 모두 분산이 크고 곡선이 비슷하지만 동일하지는 않다. 구별하려면 자료가 더 필요하다.

모집단 적률: \(\mathbb{E}[X] = (a+b)/2\), \(\mathrm{Var}(X) = (b-a)^2/12\). 따라서 \(\mathbb{E}[X^2] = (a+b)^2/4 + (b-a)^2/12\)이다.

표본을 모집단과 같다고 두면 \(\bar X = (a+b)/2\)이고 \(m_2 = (b-a)^2/12\)이며, \(m_2\)는 (\(n\)을 쓴) 표본분산이다.

두 번째 식에서 \(b - a = \sqrt{12 m_2}\)이다. \(a + b = 2\bar X\)와 결합하면:

\(\hat a_{\text{MoM}} = \bar X - \sqrt{3 m_2}\), \(\hat b_{\text{MoM}} = \bar X + \sqrt{3 m_2}\).

문제점: 어떤 관측값이 \([\hat a, \hat b]\) 밖에 있으면 이 추정값은 불가능하다. MLE는 순서통계량으로 이를 처리한다: \(\hat a_{\text{MLE}} = \min X_i\), \(\hat b_{\text{MLE}} = \max X_i\)로 언제나 실현 가능하다.

MLE가 지키는 경계 제약을 적률법이 어길 수 있음을 보여 준다.

Beta\((\alpha, \beta)\): \(\mathbb{E}[X] = \alpha/(\alpha+\beta)\), \(\mathrm{Var}(X) = \alpha\beta/[(\alpha+\beta)^2(\alpha+\beta+1)]\).

\(m = \bar X\), \(v = m_2\)(표본분산)라 하고 풀면:

\(\alpha + \beta = m(1 - m)/v - 1\). 따라서 \(\hat\alpha = m \cdot (m(1-m)/v - 1)\), \(\hat\beta = (1-m)(m(1-m)/v - 1)\).

타당성: \(v < m(1-m)\)이 필요하다(베타분포의 분산은 같은 평균을 갖는 베르누이분포의 분산으로 위에서 유계이다). 표본에서 \(v > m(1-m)\)이면 적률법이 실패한다(맞는 베타분포가 없다). 과대산포가 심한 자료에는 적률법을 쓸 수 없다는 뜻이다.

\(\mathbb{E}[X] = e^{\mu + \sigma^2/2}\), \(\mathrm{Var}(X) = (e^{\sigma^2} - 1) e^{2\mu + \sigma^2}\).

변동계수의 제곱을 계산하면 \(\mathrm{CV}^2 = \mathrm{Var}(X)/(\mathbb{E}[X])^2 = e^{\sigma^2} - 1\)이다.

따라서 \(\hat\sigma^2_{\text{MoM}} = \ln(1 + \mathrm{CV}^2_{\text{sample}}) = \ln(1 + s^2/\bar X^2)\)이다.

\(\hat\mu_{\text{MoM}} = \ln \bar X - \hat\sigma^2_{\text{MoM}}/2\).

더 쉬운 대안 (MLE 방식): 먼저 변환한다. \(Y_i = \ln X_i\)가 i.i.d. \(N(\mu, \sigma^2)\)이므로 표준적으로 \(\hat\mu = \bar Y\), \(\hat\sigma^2 = s_Y^2\)이다. 적률법보다 깔끔하고 효율적이다. 원래 척도에서 적률법을 쓰면 편리한 로그변환을 포기하는 셈이다.

적률법: \(\bar X = \theta/2 \Rightarrow \hat\theta_{\text{MoM}} = 2\bar X\). \(\mathrm{Var}(\hat\theta_{\text{MoM}}) = 4 \mathrm{Var}(\bar X) = 4 \theta^2/(12n) = \theta^2/(3n)\).

MLE: \(\hat\theta_{\text{MLE}} = X_{(n)}\). 편향 보정한 \(((n+1)/n) X_{(n)}\)은 \(\mathrm{Var} = \theta^2/[n(n+2)]\)이다.

비: \(\mathrm{Var}(\hat\theta_{\text{MoM}})/\mathrm{Var}(\hat\theta_{\text{MLE}}) = (n+2)/3 \to \infty\).

이 문제에서 \(n \to \infty\)일 때 적률법은 MLE보다 무한히 나쁘다. 표본평균은 최댓값에 관한 정보를 버리지만 MLE는 그것을 활용한다.

일반적인 교훈: 적률법은 소박하며 저차 적률만 사용한다. MLE는 완전한 가능도를 사용하므로, 분포가 적률 너머의 "구조"를 가질 때, 특히 지지집합이 유계이거나(균등분포의 끝점) 꼬리가 두꺼울 때 훨씬 효율적일 수 있다.

적률법은 MLE를 다루기 어렵거나 반복적 MLE 최적화의 출발점이 필요할 때 여전히 유용하다.

\(\theta \in \mathbb{R}^k\)이고 \(j = 1, \ldots, m\)에 대해 \(m \ge k\)개의 적률 조건 \(\mathbb{E}[g_j(X; \theta)] = 0\)이 있다고 하자.

GMM 추정량:

여기서 \(W\)는 양의 정부호 \(m \times m\) 가중행렬이다.

최적 \(W\): 적률 조건의 공분산행렬의 역행렬 \(W^* = \Omega^{-1}\)이다. 이것이 점근분산을 최소화한다.

2단계 GMM:

1. \(W = I\)(단위행렬)로 초기 \(\hat\theta^{(1)}\)을 구한다.
2. \(\hat\theta^{(1)}\)에서 \(\Omega\)를 추정하고 \(W = \hat\Omega^{-1}\)로 둔다.
3. 다시 최적화한다.

GMM은 실증경제학의 토대이며(Hansen의 1982년 논문은 노벨상을 받았다) 내생변수보다 도구변수가 많을 때의 도구변수 추정을 떠받친다.

적률. 앞서 본 대로

이다.

\(\lambda\)를 소거한다. 변동계수를 만들면 \(\lambda\)가 약분된다.

닫힌 해가 없는 이유. 우변이 감마함수의 비라 \(k\)에 대해 대수적으로 풀리지 않는다. \(h\)는 \(k\)의 순감소함수(\(k\)가 크면 분포가 덜 퍼진다)이므로 해는 유일하다.

푸는 방법.

1. 표본에서 \(\hat c = s/\bar x\)(변동계수)를 계산한다.
2. \(h(k) = \hat c^2\)을 수치적으로 푼다. 이분법이나 뉴턴법이면 충분하며, \(h\)가 단조라 수렴이 안정적이다.
3. 그다음 \(\hat\lambda = \bar x/\Gamma(1+1/\hat k)\)로 얻는다.

from scipy.optimize import brentq
from scipy.special import gamma
h = lambda k: gamma(1 + 2/k) / gamma(1 + 1/k)**2 - 1
k_hat = brentq(lambda k: h(k) - (s/xbar)**2, 0.1, 50)
lam_hat = xbar / gamma(1 + 1/k_hat)

어림 공식. \(\hat k \approx \hat c^{-1.086}\)이라는 근사가 널리 쓰인다(\(0.1 < \hat c < 1\) 범위에서 꽤 정확하다). 수치해의 좋은 초기값이 된다.

덧붙임. 와이불의 MLE도 닫힌 해가 없어 수치적으로 풀어야 하므로, 적률법의 "계산이 쉽다"는 장점이 여기서는 크지 않다. 그래도 MLE 반복의 초기값으로 쓰기에 좋고, 중도절단이 없는 완전 자료에서는 결과가 크게 다르지 않다.

1단계 — 차이가 표집변동 수준인가. 부트스트랩으로 두 추정값을 함께 재표집해 차이 \(\hat\alpha_{\text{MoM}}-\hat\alpha_{\text{MLE}}\)의 분포를 본다. 그 분포가 0을 넉넉히 담으면 우연으로 설명된다. \(n=10\)이면 두 추정값의 표준오차가 모두 커서 대개 이 경우에 해당한다.

2단계 — 모형 적합도를 본다. 두 추정값으로 각각 감마 밀도를 그려 히스토그램이나 경험분포와 겹쳐 본다. Q-Q 그림이 더 예민하다. 한쪽이 뚜렷이 나쁘면 그쪽을 버린다.

3단계 — 이상치를 확인한다. 적률법은 \(s^2\)을 쓰므로 큰 값 하나에 크게 흔들린다. 최댓값을 빼고 다시 계산해 보아 적률법 추정값만 크게 움직이면 그것이 원인이다.

4단계 — 원칙적으로는 MLE. 감마 모형이 맞다면 MLE가 점근효율적이고, \(n=10\)처럼 작은 표본에서도 대개 적률법보다 MSE가 작다. 다만 MLE도 소표본에서 \(\alpha\)를 과대추정하는 편향이 알려져 있어, 편향 보정(\(\hat\alpha_{\text{corr}} = \hat\alpha(n-3)/n\) 꼴)을 고려할 만하다.

5단계 — 그러나 모형을 의심한다. 두 방법의 차이가 크고 적합도도 나쁘면 감마분포 자체가 틀렸을 가능성이 높다. 로그정규, 와이불, 역가우스를 함께 적합해 AIC로 비교한다.

보고. \(n=10\)에서는 어느 방법을 쓰든 추정값의 신뢰구간이 매우 넓다. 점추정값의 차이를 논하기 전에 구간을 먼저 보여 주는 것이 정직하다. 두 구간이 크게 겹친다면 "방법에 따라 다르다"는 논의 자체가 과잉 해석이다.

모형. \(f(x) = \pi\phi(x;\mu_1,\sigma_1^2)+(1-\pi)\phi(x;\mu_2,\sigma_2^2)\)로 모수가 5개(\(\pi,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2\))다.

적률법. 5개 모수이므로 1차부터 5차까지의 적률을 맞춰야 한다.

피어슨의 방법(1894). 그는 이 연립방정식을 대수적으로 정리해 9차 다항방정식 하나로 줄였다. 게의 이마 너비 자료가 두 종의 혼합인지를 판정하려는 것이 동기였고, 손으로 그 방정식을 풀었다. 통계사에서 혼합모형의 출발점으로 꼽히는 작업이다.

문제점.

• 5차 적률을 써야 한다. 표집변동이 엄청나서 추정값이 극도로 불안정하다. \(n\)이 수백이어도 믿기 어렵다.
• 9차 방정식의 실근이 여럿일 수 있고, 그중 유효한 것을 고르는 기준이 명확하지 않다.
• 모수공간을 벗어나기 쉽다. 음의 분산이나 \(\pi \notin [0,1]\)이 흔히 나온다.

현대적 방법 — EM 알고리즘.

• E 단계: 각 관측값이 성분 1에 속할 사후확률 \(\gamma_i = \dfrac{\pi\phi(x_i;\mu_1,\sigma_1^2)}{f(x_i)}\)를 계산한다.
• M 단계: \(\gamma_i\)를 가중치로 삼아 가중평균·가중분산으로 모수를 갱신하고 \(\hat\pi = \bar\gamma\)로 둔다.

비교.

그래도 적률법의 쓸모. EM의 초기값으로 쓸 수 있고, 성분이 잘 분리된 경우에는 빠르게 합리적인 답을 준다. 또 고차 적률을 쓰지 않는 변형(분위수 기반, 특성함수 기반)이 계산 통계에서 다시 연구되고 있다.

절차.

rng = np.random.default_rng(0)
alpha_true, beta_true = 2.0, 3.0
ns = [10, 20, 50, 100, 200, 500, 1000]
B = 5000
for n in ns:
    est = np.empty((B, 2))
    for b in range(B):
        x = rng.gamma(alpha_true, beta_true, n)
        m, v = x.mean(), x.var(ddof=1)
        est[b] = [m**2 / v, v / m]          # (alpha_hat, beta_hat)
    bias = est.mean(axis=0) - [alpha_true, beta_true]
    sd = est.std(axis=0, ddof=1)
    print(n, bias, sd)

그려야 할 것.

1. \(\log n\) 대 \(\log|\text{편향}|\). 기울기가 편향의 수렴 속도를 준다. 적률법의 편향은 보통 \(O(1/n)\)이므로 기울기 \(-1\)인 직선이 나와야 한다. 그렇지 않으면 일치성이 의심스럽다.
2. \(\log n\) 대 \(\log(\text{표준편차})\). 기울기가 \(-1/2\)이어야 한다. 정칙 문제의 표준 속도다.
3. \(n\cdot\operatorname{Var}\) 대 \(n\). 이 값이 상수로 안정되면 그 상수가 점근분산이고, 이론값과 맞는지 확인할 수 있다. 가장 예민한 진단이다.
4. 각 \(n\)에서의 추정값 분포. 상자그림이나 히스토그램. 치우침이 언제 사라지는지, 정규근사가 언제부터 통하는지를 눈으로 본다.
5. MLE와의 MSE 비. \(n\)에 대해 그리면 점근상대효율로 수렴하는지 볼 수 있다.

주의할 점.

• 모수공간을 벗어난 추정값을 어떻게 처리했는지 반드시 기록한다. 버리면 편향이 왜곡되고, 자르면 그것도 하나의 추정량이다.
• \(B\)를 충분히 크게. 편향은 분산보다 작은 양이므로 정밀하게 재려면 \(B\)가 커야 한다. 편향의 표준오차가 \(\operatorname{SD}/\sqrt B\)임을 기억한다.
• 여러 참값에서 반복한다. \(\alpha=2\)에서의 결론이 \(\alpha=0.5\)나 \(\alpha=20\)에서도 성립하는지 확인한다. 감마분포는 \(\alpha\)에 따라 모양이 크게 달라진다.