1 개요
점추정 개요에서 추정량 탐색과 평가의 전체 구조를 조망했다. 이 포스트에서는 추정량 탐색 방법(Methods of Finding Estimators) 을 세분화하여, 네 가지 체계적 방법의 원리, 예시, 장단점을 깊이 있게 다룬다 (Casella & Berger, 2002, Ch.7.2).
직관만으로 좋은 추정량을 찾는 것은 단순한 모형에서만 가능하다. 복잡한 모형에서는 체계적인 방법이 필요하며, 각 방법은 서로 다른 철학과 장단점을 가진다.
| 방법 | 철학 | 핵심 아이디어 |
|---|---|---|
| 적률법 | 표본 \(\approx\) 모집단 | 표본 적률 = 모집단 적률로 놓고 풀기 |
| MLE | 관측 데이터가 가장 그럴듯한 모수 | 우도함수 최대화 |
| 베이즈 | 사전 믿음 + 데이터 → 사후 믿음 | 사후분포에서 추출 |
| EM | 잠재 변수가 있으면 반복으로 풀기 | E-step/M-step 반복 |
이 네 방법은 경쟁 관계가 아니라 상보적이다 — 적률법으로 초기값을 구하고, MLE를 수치적으로 최적화하고, 베이즈로 불확실성을 정량화하고, EM으로 잠재 변수 문제를 해결한다.
2 적률법 (Method of Moments)
2.1 원리
Karl Pearson(1890년대)이 제안한 가장 오래된 추정 방법이다. 모집단 적률을 표본 적률과 같다고 놓고, 연립방정식을 풀어 모수를 추정한다.
\(k\) 개의 미지 모수 \((\theta_1, \ldots, \theta_k)\) 가 있을 때, \(j\) 번째 표본 적률과 모집단 적률을 정의한다:
\[ m_j = \frac{1}{n}\sum_{i=1}^n X_i^j \quad (\text{표본}), \qquad \mu'_j(\boldsymbol{\theta}) = E[X^j] \quad (\text{모집단}) \]
적률법 추정량 \((\tilde{\theta}_1, \ldots, \tilde{\theta}_k)\) 는 연립방정식
\[ m_j = \mu'_j(\theta_1, \ldots, \theta_k), \quad j = 1, \ldots, k \]
의 해이다.
직관: “대수의 법칙에 의해 \(m_j \xrightarrow{P} \mu'_j\) 이므로, \(m_j = \mu'_j\) 로 놓는 것은 대표본에서 합리적이다.”
2.2 예시 1: 정규분포
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\mu, \sigma^2)\) — 모수 2개이므로 처음 2개 적률을 사용한다.
\[ m_1 = \bar{X} = \mu, \quad m_2 = \frac{1}{n}\sum X_i^2 = \mu^2 + \sigma^2 \]
풀면
\[ \tilde{\mu} = \bar{X}, \qquad \tilde{\sigma}^2 = m_2 - m_1^2 = \frac{1}{n}\sum(X_i - \bar{X})^2 \]
직관과 일치한다 — 표본평균으로 모평균을, 표본분산(비편향 아닌 버전)으로 모분산을 추정한다.
2.3 예시 2: 이항분포 (미지의 \(k\) 와 \(p\))
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} \text{Binomial}(k, p)\) 에서 \(k\) 와 \(p\) 가 모두 미지일 때 — 범죄 통계에서 미신고 범죄의 총 발생 건수(\(k\))와 신고율(\(p\))을 추정하는 상황이다.
\(E[X] = kp\), \(\text{Var}(X) = kp(1-p)\) 이므로
\[ \bar{X} = kp, \quad \frac{1}{n}\sum(X_i - \bar{X})^2 = kp(1-p) \]
풀면
\[ \tilde{p} = 1 - \frac{\frac{1}{n}\sum(X_i - \bar{X})^2}{\bar{X}}, \qquad \tilde{k} = \frac{\bar{X}}{\tilde{p}} \]
문제점: \(\bar{X}\) 가 표본분산보다 작으면 \(\tilde{p}\) 가 음수가 된다. 모수 범위(\(0 < p < 1\), \(k \geq 1\))를 벗어나는 추정값이 산출될 수 있다는 것이 적률법의 대표적 한계이다.
2.4 예시 3: Satterthwaite 근사
적률법은 추정량 탐색 외에 분포 근사에도 사용된다. 독립인 카이제곱 확률변수의 가중합 \(\sum a_i Y_i\) (\(Y_i \sim \chi^2_{r_i}\))의 분포를 \(\chi^2_\nu / \nu\) 로 근사할 때, 첫 두 적률을 맞추면
\[ \hat{\nu} = \frac{(\sum a_i Y_i)^2}{\sum \frac{a_i^2}{r_i} Y_i^2} \]
이것이 Satterthwaite(1946) 근사이다. Welch의 \(t\)-검정 등에서 현재까지 널리 사용된다.
구체적 수치 예시 (Welch’s t-검정): 두 그룹의 표본이 \(n_1 = 10\), \(s_1^2 = 4.0\), \(n_2 = 15\), \(s_2^2 = 1.5\) 일 때, 각 그룹의 분산 추정량 \(S_i^2/(n_i - 1) \sim \sigma_i^2/(n_i - 1) \cdot \chi^2_{n_i-1}/(n_i-1)\) 을 이용한다. Welch의 \(t\)-검정에서는 \(W = S_1^2/n_1 + S_2^2/n_2\) 를 분모로 사용하며,
\[ a_1 = \frac{1}{n_1} = 0.1, \quad a_2 = \frac{1}{n_2} = 0.0\overline{6}, \quad Y_1 = S_1^2, \; r_1 = n_1 - 1 = 9, \quad Y_2 = S_2^2, \; r_2 = 14 \]
를 대입하면
\[ \hat{\nu} = \frac{(0.1 \times 4.0 + 0.0\overline{6} \times 1.5)^2}{\frac{(0.1 \times 4.0)^2}{9} + \frac{(0.0\overline{6} \times 1.5)^2}{14}} = \frac{(0.5)^2}{0.16/9 + 0.01/14} \approx \frac{0.25}{0.0178 + 0.000714} \approx 13.5 \]
실제 자유도 9나 14 대신 13.5를 사용하여 \(t\)-분포의 임계값을 구한다. 두 그룹의 분산이 다를수록 \(\hat{\nu}\) 는 작아지며, 이는 더 보수적인 기각역으로 이어진다.
2.5 장단점 요약
| 장점 | 단점 |
|---|---|
| 계산이 단순하다 — 연립방정식만 풀면 된다 | 최적이 아닌 경우가 많다 (MLE보다 분산이 크다) |
| 거의 항상 해가 존재한다 | 모수 범위를 벗어나는 추정값이 가능하다 |
| MLE/EM의 초기값으로 유용하다 | 고차 적률의 분산이 커서 불안정할 수 있다 |
| 분포 근사에도 활용된다 | 충분통계량을 활용하지 않는다 |
3 최대우도추정법 (Maximum Likelihood Estimation)
3.1 정의
가장 널리 사용되는 추정 방법이다. 우도원리와 직결된다.
각 표본점 \(\mathbf{x}\) 에 대해, 우도함수 \(L(\theta|\mathbf{x}) = \prod_{i=1}^n f(x_i|\theta)\) 를 \(\theta\) 의 함수로서 최대화하는 값 \(\hat{\theta}(\mathbf{x})\) 를 MLE 라 한다.
직관: “이 데이터가 발생할 확률(밀도)을 가장 높게 만드는 모수값을 선택한다.”
MLE의 치역은 모수 공간 \(\Theta\) 와 일치한다 — 정의상 \(\Theta\) 위에서 최대화하기 때문이다. 이것은 적률법과의 중요한 차이이다.
3.2 MLE 찾기: 미분법
우도함수가 미분 가능하면, MLE 후보는 스코어 방정식의 해이다:
\[ \frac{\partial}{\partial \theta_i} L(\theta|\mathbf{x}) = 0, \quad i = 1, \ldots, k \]
실무에서는 곱의 미분 대신 로그우도의 미분이 훨씬 쉽다:
\[ \frac{\partial}{\partial \theta_i} \ell(\theta|\mathbf{x}) = 0, \quad \ell(\theta|\mathbf{x}) = \log L(\theta|\mathbf{x}) = \sum_{j=1}^n \log f(x_j|\theta) \]
\(\log\) 가 순증가함수이므로 \(L\) 과 \(\ell\) 의 극값은 일치한다.
1도 조건은 극값의 필요조건이지 충분조건이 아니다. 극소점, 변곡점, 국소 극대점이 모두 스코어 방정식을 만족한다. 경계(boundary)에서의 극값은 1도 미분이 0이 아닐 수 있다. 반드시 전역 최대임을 확인해야 한다.
3.3 예시 1: 정규분포 (\(\theta\) 만 미지)
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\theta, 1)\) 에서
\[ L(\theta|\mathbf{x}) = (2\pi)^{-n/2} \exp\!\left(-\frac{1}{2}\sum(x_i - \theta)^2\right) \]
스코어 방정식: \(\sum(x_i - \theta) = 0\) → \(\hat{\theta} = \bar{x}\)
전역 최대 확인: \(\sum(x_i - a)^2 \geq \sum(x_i - \bar{x})^2\) (모든 \(a\))이므로 \(\exp(-\frac{1}{2}\sum(x_i - \theta)^2)\) 는 \(\theta = \bar{x}\) 에서 최대이다. 경계(\(\pm \infty\))에서 우도는 0이다.
3.4 예시 2: 베르누이 분포
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} \text{Bernoulli}(p)\) 에서 \(y = \sum x_i\) 로 놓으면
\[ \ell(p|\mathbf{x}) = y \log p + (n-y) \log(1-p) \]
\(0 < y < n\) 이면: \(d\ell/dp = y/p - (n-y)/(1-p) = 0\) → \(\hat{p} = y/n\)
\(y = 0\) 이면: \(\ell = n\log(1-p)\) 은 \(p\) 에 대해 순감소 → \(\hat{p} = 0\)
\(y = n\) 이면: \(\ell = n\log p\) 은 \(p\) 에 대해 순증가 → \(\hat{p} = 1\)
모든 경우 \(\hat{p} = y/n = \sum X_i / n\) 이다.
3.5 예시 3: 정규분포 (\(\mu, \sigma^2\) 모두 미지)
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\mu, \sigma^2)\) 에서 로그우도를 \(\mu\) 와 \(\sigma^2\) 에 대해 편미분하면
\[ \hat{\mu} = \bar{X}, \qquad \hat{\sigma}^2 = \frac{1}{n}\sum(X_i - \bar{X})^2 \]
MLE \(\hat{\sigma}^2\) 는 \(n\) 으로 나누므로 비편향이 아니다 (\(E[\hat{\sigma}^2] = \frac{n-1}{n}\sigma^2\)). 비편향 추정량은 \(S^2 = \frac{1}{n-1}\sum(X_i - \bar{X})^2\) 이다. 이것은 MLE가 반드시 비편향이 아니라는 사실을 보여준다.
프로파일 우도(profile likelihood) 기법: \(\mu\) 를 먼저 최대화하면 (\(\hat{\mu} = \bar{X}\)), 남은 문제가 \(\sigma^2\) 에 대한 1차원 최대화로 축소된다. 이 방법은 다변량 MLE 문제에서 2도 편미분 조건을 직접 확인하는 것보다 훨씬 효율적이다.
3.6 모수 범위 제한이 있는 경우
\(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\theta, 1)\) 에서 \(\theta \geq 0\) 이 알려져 있다면
\[ \hat{\theta} = \begin{cases} \bar{X} & \text{if } \bar{X} \geq 0 \\ 0 & \text{if } \bar{X} < 0 \end{cases} \]
\(\bar{X} < 0\) 이면 우도함수가 \(\theta \geq 0\) 에서 순감소하므로 경계값 \(\theta = 0\) 이 최대점이다. MLE는 항상 모수 범위 안에 있다.
3.7 MLE의 불변성 원리
\(\hat{\theta}\) 가 \(\theta\) 의 MLE이면, 임의의 함수 \(\tau(\theta)\) 에 대해 \(\tau(\hat{\theta})\) 가 \(\tau(\theta)\) 의 MLE이다.
증명 핵심: \(\tau\) 가 일대일이 아닐 수 있으므로, 유도 우도(induced likelihood) 를 정의한다:
\[ L^*(\eta|\mathbf{x}) = \sup_{\{\theta : \tau(\theta) = \eta\}} L(\theta|\mathbf{x}) \]
\(L^*\) 의 최대점 \(\hat{\eta}\) 에 대해
\[ L^*(\hat{\eta}|\mathbf{x}) = \sup_\eta \sup_{\{\theta:\tau(\theta)=\eta\}} L(\theta|\mathbf{x}) = \sup_\theta L(\theta|\mathbf{x}) = L(\hat{\theta}|\mathbf{x}) = L^*(\tau(\hat{\theta})|\mathbf{x}) \]
따라서 \(\hat{\eta} = \tau(\hat{\theta})\) 이다.
응용 예시: 정규분포에서 \(\hat{\mu} = \bar{X}\) 이면, \(\mu^2\) 의 MLE는 \(\bar{X}^2\), \(\sin(\mu)\) 의 MLE는 \(\sin(\bar{X})\) 이다.
3.8 MLE의 수치적 불안정성
MLE는 우도함수의 최대화로 구하므로, 수치적 민감도(numerical sensitivity) 문제가 있을 수 있다.
\(X_1, \ldots, X_5 \overset{\text{iid}}{\sim} \text{Binomial}(k, p)\), \(k\) 와 \(p\) 모두 미지일 때 (Olkin, Petkau, & Zidek, 1981):
- 데이터 \((16, 18, 22, 25, 27)\) → \(\hat{k} = 99\)
- 데이터 \((16, 18, 22, 25, 28)\) → \(\hat{k} = 190\)
27을 28로 바꾸었을 뿐인데 \(\hat{k}\) 가 거의 두 배로 변한다. 이런 현상은 우도함수가 최대점 근방에서 매우 평탄할 때 발생한다.
3.9 MLE 장단점 요약
| 장점 | 단점 |
|---|---|
| 점근적으로 최적 (일치, 정규, 효율) | 소표본에서 편향이 있을 수 있다 |
| 불변성 원리가 성립한다 | 전역 최대를 보장하기 어려울 수 있다 |
| 모수 범위를 항상 만족한다 | 수치적으로 불안정할 수 있다 |
| 우도원리와 자연스럽게 연결된다 | 정칙 조건이 필요하다 |
4 베이즈 추정 (Bayesian Estimation)
4.1 빈도주의와의 철학적 차이
빈도주의에서 \(\theta\) 는 고정된 미지의 상수이다. 베이지안에서 \(\theta\) 는 변동하는 확률변수이며, 그 분포로 불확실성을 표현한다.
4.2 사전-사후 업데이트
\[ \pi(\theta|\mathbf{x}) = \frac{f(\mathbf{x}|\theta) \, \pi(\theta)}{m(\mathbf{x})}, \quad m(\mathbf{x}) = \int f(\mathbf{x}|\theta)\pi(\theta) \, d\theta \]
- \(\pi(\theta)\): 사전분포 — 데이터 이전의 \(\theta\) 에 대한 믿음
- \(f(\mathbf{x}|\theta) = L(\theta|\mathbf{x})\): 우도 — 데이터가 제공하는 정보
- \(\pi(\theta|\mathbf{x})\): 사후분포 — 데이터 이후 업데이트된 믿음
- \(m(\mathbf{x})\): 주변 우도 — 정규화 상수
사후분포에서 점추정량을 추출한다:
| 추정량 | 최적화 기준 | 수식 |
|---|---|---|
| 사후 평균 | 제곱 손실 최소화 | \(E[\theta|\mathbf{x}]\) |
| 사후 중앙값 | 절대 손실 최소화 | \(\text{med}(\theta|\mathbf{x})\) |
| 사후 최빈값 (MAP) | 0-1 손실 최소화 | \(\arg\max_\theta \pi(\theta|\mathbf{x})\) |
4.3 예시 1: 이항-베타 (켤레 모형)
\(Y = \sum X_i \sim \text{Binomial}(n, p)\), 사전분포 \(p \sim \text{Beta}(\alpha, \beta)\) 일 때
사후분포는
\[ \pi(p|y) = \frac{\Gamma(n+\alpha+\beta)}{\Gamma(y+\alpha)\Gamma(n-y+\beta)} p^{y+\alpha-1}(1-p)^{n-y+\beta-1} \]
이것은 \(\text{Beta}(y+\alpha, n-y+\beta)\) 이다. 베이즈 추정량(사후 평균)은
\[ \hat{p}_B = \frac{y + \alpha}{\alpha + \beta + n} \]
이것을 가중평균으로 재배열하면 구조가 명확해진다:
\[ \hat{p}_B = \underbrace{\frac{n}{\alpha+\beta+n}}_{\text{데이터 가중치}} \cdot \underbrace{\frac{y}{n}}_{\text{표본비율}} + \underbrace{\frac{\alpha+\beta}{\alpha+\beta+n}}_{\text{사전 가중치}} \cdot \underbrace{\frac{\alpha}{\alpha+\beta}}_{\text{사전 평균}} \]
베이즈 추정량은 사전 평균과 표본비율의 가중평균이다. \(n\) 이 커지면 데이터 가중치가 1에 가까워져 MLE에 수렴한다.
4.4 예시 2: 정규-정규 (켤레 모형)
\(X \sim N(\theta, \sigma^2)\), 사전분포 \(\theta \sim N(\mu, \tau^2)\) (\(\sigma^2, \mu, \tau^2\) 기지)일 때
\[ E[\theta|x] = \frac{\tau^2}{\tau^2 + \sigma^2} x + \frac{\sigma^2}{\sigma^2 + \tau^2} \mu, \qquad \text{Var}(\theta|x) = \frac{\sigma^2 \tau^2}{\sigma^2 + \tau^2} \]
사후분포도 정규분포이다 — 정규족은 자신의 켤레 가족이다.
사전 분산 \(\tau^2\) 의 역할:
- \(\tau^2 \to \infty\) (사전 정보 모호): \(E[\theta|x] \to x\) (MLE로 수렴)
- \(\tau^2 \to 0\) (사전 정보 강함): \(E[\theta|x] \to \mu\) (사전 평균에 고정)
- \(\sigma^2 > \tau^2\) (사전이 데이터보다 정밀): 사전 평균에 더 많은 가중
4.5 켤레 가족 (Conjugate Family)
사전분포의 클래스 \(\Pi\) 가 표본 분포 가족 \(\mathcal{F}\) 의 켤레 가족(conjugate family) 이란, 모든 사전분포 \(\pi \in \Pi\) 와 모든 표본 \(\mathbf{x}\) 에 대해 사후분포도 \(\Pi\) 에 속하는 것이다.
| 표본 분포 | 켤레 사전분포 | 사후분포 |
|---|---|---|
| \(\text{Binomial}(n, p)\) | \(\text{Beta}(\alpha, \beta)\) | \(\text{Beta}(y+\alpha, n-y+\beta)\) |
| \(\text{Poisson}(\lambda)\) | \(\text{Gamma}(\alpha, \beta)\) | \(\text{Gamma}(\sum x_i + \alpha, n + \beta)\) |
| \(N(\mu, \sigma^2_0)\) | \(N(\mu_0, \tau^2)\) | \(N(\mu_{\text{post}}, \sigma^2_{\text{post}})\) |
| \(\text{Exp}(\theta)\) | \(\text{Gamma}(\alpha, \beta)\) | \(\text{Gamma}(n + \alpha, \sum x_i + \beta)\) |
켤레 가족의 장점은 사후분포가 닫힌 형태(closed-form) 로 구해진다는 것이다. 사전 → 사후 업데이트가 모수(hyperparameter)의 업데이트로 환원된다.
5 EM 알고리즘 (Expectation-Maximization Algorithm)
5.1 동기: 불완전 자료에서의 MLE
관측 데이터 \(\mathbf{Y}\) 와 관측되지 않은 잠재 데이터 \(\mathbf{X}\) 가 있을 때, 완전 데이터 \((\mathbf{Y}, \mathbf{X})\) 의 우도는 다루기 쉽지만 불완전 데이터 \(\mathbf{Y}\) 만의 우도는 다루기 어려운 상황이 빈번하다.
\[ g(\mathbf{y}|\theta) = \int f(\mathbf{y}, \mathbf{x}|\theta) \, d\mathbf{x} \]
\(L(\theta|\mathbf{y}) = g(\mathbf{y}|\theta)\) 를 직접 최대화하기 어렵지만, \(L(\theta|\mathbf{y}, \mathbf{x}) = f(\mathbf{y}, \mathbf{x}|\theta)\) 는 최대화하기 쉬울 수 있다. EM 알고리즘은 이 아이디어를 활용한다.
5.2 E-step과 M-step
초기값 \(\theta^{(0)}\) 에서 시작하여, 다음을 반복한다:
E-step (Expectation): 현재 추정값 \(\theta^{(t)}\) 하에서, 완전 데이터 로그우도의 조건부 기댓값을 계산한다:
\[ Q(\theta|\theta^{(t)}) = E_{\mathbf{X}|\mathbf{y}, \theta^{(t)}}[\log L(\theta|\mathbf{y}, \mathbf{X})] \]
M-step (Maximization): \(Q\) 함수를 \(\theta\) 에 대해 최대화한다:
\[ \theta^{(t+1)} = \arg\max_\theta Q(\theta|\theta^{(t)}) \]
5.3 왜 작동하는가: 핵심 항등식
로그우도를 분해하면
\[ \log L(\theta|\mathbf{y}) = \log L(\theta|\mathbf{y}, \mathbf{x}) - \log k(\mathbf{x}|\theta, \mathbf{y}) \]
여기서 \(k(\mathbf{x}|\theta, \mathbf{y}) = f(\mathbf{y}, \mathbf{x}|\theta) / g(\mathbf{y}|\theta)\) 는 \(\mathbf{X}\) 의 조건부 분포이다.
양변에 \(E_{\mathbf{X}|\mathbf{y}, \theta^{(t)}}\) 를 취하면, 좌변의 \(\log L(\theta|\mathbf{y})\) 는 \(\mathbf{X}\) 에 무관하므로 그대로이고
\[ \log L(\theta|\mathbf{y}) = Q(\theta|\theta^{(t)}) - H(\theta|\theta^{(t)}) \]
\(H(\theta|\theta^{(t)}) = E[\log k(\mathbf{X}|\theta, \mathbf{y})|\theta^{(t)}, \mathbf{y}]\) 에 대해 Jensen 부등식에 의해 \(H(\theta|\theta^{(t)}) \leq H(\theta^{(t)}|\theta^{(t)})\) 이므로, \(Q\) 를 증가시키면 \(L\) 도 증가한다.
5.4 단조 수렴 정리
EM 반복 수열 \(\{\theta^{(t)}\}\) 는 다음을 만족한다:
\[ L(\theta^{(t+1)}|\mathbf{y}) \geq L(\theta^{(t)}|\mathbf{y}) \]
등호는 연속된 반복이 동일한 \(Q\) 최대값을 산출할 때만 성립한다.
불완전 데이터 우도가 매 반복마다 증가하거나 유지된다. 따라서 수열은 수렴하지만, 전역 최대가 아닌 국소 최대에 수렴할 수 있다. 다른 초기값에서 여러 번 실행하는 것이 권장된다.
5.5 예시: 가우시안 혼합 모형 (Gaussian Mixture Model)
\(K\) 개의 정규분포 혼합에서의 EM이 가장 대표적인 응용이다.
\[ f(x|\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\sigma}^2) = \sum_{k=1}^K \pi_k \, \phi(x|\mu_k, \sigma_k^2) \]
- 잠재 변수: \(Z_i \in \{1, \ldots, K\}\) — 각 관측값의 클러스터 소속
- E-step: \(\gamma_{ik} = P(Z_i = k | x_i, \theta^{(t)})\) — 사후 소속 확률 계산
- M-step: 가중 표본평균/분산으로 \(\mu_k, \sigma_k^2, \pi_k\) 업데이트
\[ \mu_k^{(t+1)} = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}, \quad (\sigma_k^2)^{(t+1)} = \frac{\sum_i \gamma_{ik}(x_i - \mu_k^{(t+1)})^2}{\sum_i \gamma_{ik}} \]
5.6 EM의 주요 응용
| 응용 | 관측 데이터 | 잠재 변수 |
|---|---|---|
| 혼합 모형 | 관측값 \(x_i\) | 클러스터 소속 \(z_i\) |
| 결측치 처리 | 관측된 값 | 결측된 값 |
| Hidden Markov Model | 관측 시퀀스 | 은닉 상태 시퀀스 |
| 요인 분석 | 관측 변수 | 잠재 요인 |
| 생존 분석 (절단 데이터) | 관측 시간 | 절단된 사건 시간 |
6 네 방법의 비교
| 기준 | 적률법 | MLE | 베이즈 | EM |
|---|---|---|---|---|
| 철학 | 빈도주의 | 빈도주의 | 베이지안 | 빈도주의 (MLE 보조) |
| 계산 난이도 | 낮음 | 중간~높음 | 중간~높음 | 높음 (반복) |
| 점근 최적성 | 아니오 | 예 (효율적) | — | MLE에 수렴 |
| 모수 범위 보장 | 아니오 | 예 | 예 | 예 |
| 사전 정보 활용 | 아니오 | 아니오 | 예 | 아니오 |
| 잠재 변수 처리 | 아니오 | 어려움 | 가능 | 전문 |
| 소표본 성능 | 보통 | 보통 | 좋음 (좋은 사전 시) | MLE 수준 |
7 코드 예시
7.1 Step 1: 순수 Python 구현 (네 방법 비교)
감마분포에서 네 방법을 모두 구현하여 비교한다.
import math
import random
random.seed(42)
# 감마분포 Gamma(alpha=3, beta=2) 에서 표본 생성
# 역변환법 대신 미리 준비된 데이터 사용
data = [4.2, 6.1, 3.8, 7.3, 5.5, 2.9, 8.1, 4.7, 6.8, 3.2,
5.1, 7.9, 4.4, 6.3, 3.6, 5.8, 7.1, 4.0, 6.5, 3.4]
n = len(data)
x_bar = sum(data) / n
x_var = sum((x - x_bar)**2 for x in data) / n
# === 1. 적률법 ===
# Gamma(alpha, beta): E[X] = alpha*beta, Var(X) = alpha*beta^2
beta_mom = x_var / x_bar
alpha_mom = x_bar / beta_mom
print("=== 적률법 (Method of Moments) ===")
print(f" alpha = {alpha_mom:.4f}, beta = {beta_mom:.4f}")
# === 2. MLE (고정점 반복) ===
# MLE for beta = x_bar / alpha
# alpha의 프로파일: psi(alpha) - log(alpha) + log(x_bar) - mean(log(x)) = 0
mean_log_x = sum(math.log(x) for x in data) / n
alpha_mle = alpha_mom # 적률법으로 초기값
for _ in range(100):
# digamma 근사: psi(a) ≈ log(a) - 1/(2a) - 1/(12a^2)
psi_a = math.log(alpha_mle) - 1/(2*alpha_mle) - 1/(12*alpha_mle**2)
score = psi_a - math.log(x_bar) + math.log(alpha_mle) - mean_log_x
# trigamma: psi'(a) ≈ 1/a + 1/(2a^2) + 1/(6a^3)
deriv = 1/alpha_mle + 1/(2*alpha_mle**2) + 1/(6*alpha_mle**3) + 1/alpha_mle
alpha_new = alpha_mle - score / deriv
if abs(alpha_new - alpha_mle) < 1e-12:
break
alpha_mle = max(alpha_new, 0.01)
beta_mle = x_bar / alpha_mle
print(f"\n=== MLE (뉴턴-랩슨) ===")
print(f" alpha = {alpha_mle:.4f}, beta = {beta_mle:.4f}")
# === 3. 베이즈 (약정보 사전분포) ===
# Gamma 우도의 켤레 사전분포는 복잡하므로, 단순화하여
# alpha 기지(=MLE), beta에 대해 Gamma(a0, b0) 사전 → Gamma 사후
# beta의 사후: Gamma(n*alpha + a0, sum(x) + b0)
a0, b0 = 0.01, 0.01 # 약정보 사전
alpha_bayes = alpha_mle # alpha는 MLE 사용
a_post = n * alpha_bayes + a0
b_post = sum(data) + b0
beta_bayes = a_post / b_post # 사후 평균 (역모수화 주의)
# 실제로 Gamma(shape, rate) 켤레에서 beta = shape/rate
beta_bayes = (n * alpha_bayes + a0) / (sum(data) + b0)
print(f"\n=== 베이즈 (약정보 사전) ===")
print(f" alpha = {alpha_bayes:.4f} (MLE 사용)")
print(f" beta_posterior_mean = {beta_bayes:.4f}")
# === 비교 ===
print(f"\n=== 비교 (참값: alpha=3.0, beta=2.0) ===")
print(f" {'방법':12s} | {'alpha':>8s} | {'beta':>8s}")
print(f" {'-'*12}-+-{'-'*8}-+-{'-'*8}")
print(f" {'적률법':12s} | {alpha_mom:8.4f} | {beta_mom:8.4f}")
print(f" {'MLE':12s} | {alpha_mle:8.4f} | {beta_mle:8.4f}")
print(f" {'참값':12s} | {'3.0000':>8s} | {'2.0000':>8s}")7.2 Step 2: scipy/numpy 구현 (EM 알고리즘 — 가우시안 혼합)
import numpy as np
from scipy.stats import norm
np.random.seed(42)
# 2-component Gaussian Mixture 데이터 생성
n1, n2 = 200, 300
mu_true = [2.0, 7.0]
sigma_true = [1.0, 1.5]
pi_true = [0.4, 0.6]
x1 = np.random.normal(mu_true[0], sigma_true[0], n1)
x2 = np.random.normal(mu_true[1], sigma_true[1], n2)
data = np.concatenate([x1, x2])
np.random.shuffle(data)
n = len(data)
# EM 알고리즘
K = 2
# 초기값 (적률법 기반: 데이터를 반으로 나눠 추정)
mu = np.array([np.mean(data[:n//2]), np.mean(data[n//2:])])
sigma2 = np.array([np.var(data[:n//2]), np.var(data[n//2:])])
pi_k = np.array([0.5, 0.5])
print("EM 알고리즘: 2-component Gaussian Mixture")
print(f"참값: mu={mu_true}, sigma={sigma_true}, pi={pi_true}\n")
for iteration in range(50):
# === E-step: 사후 소속 확률 계산 ===
gamma = np.zeros((n, K))
for k in range(K):
gamma[:, k] = pi_k[k] * norm.pdf(data, mu[k], np.sqrt(sigma2[k]))
gamma /= gamma.sum(axis=1, keepdims=True)
# === M-step: 모수 업데이트 ===
N_k = gamma.sum(axis=0)
mu_new = (gamma.T @ data) / N_k
sigma2_new = np.array([
np.sum(gamma[:, k] * (data - mu_new[k])**2) / N_k[k]
for k in range(K)
])
pi_new = N_k / n
# 수렴 확인
if np.max(np.abs(mu_new - mu)) < 1e-8:
break
mu, sigma2, pi_k = mu_new, sigma2_new, pi_new
print(f"수렴: {iteration+1}회 반복")
print(f"\n {'':5s} | {'mu':>8s} | {'sigma':>8s} | {'pi':>8s}")
print(f" {'-'*5}-+-{'-'*8}-+-{'-'*8}-+-{'-'*8}")
for k in range(K):
print(f" k={k+1:2d} | {mu[k]:8.4f} | {np.sqrt(sigma2[k]):8.4f} | {pi_k[k]:8.4f}")
print(f"\n 참값:")
for k in range(K):
print(f" k={k+1:2d} | {mu_true[k]:8.4f} | {sigma_true[k]:8.4f} | {pi_true[k]:8.4f}")
# 로그우도 계산
log_lik = np.sum(np.log(
sum(pi_k[k] * norm.pdf(data, mu[k], np.sqrt(sigma2[k])) for k in range(K))
))
print(f"\n 최종 로그우도: {log_lik:.2f}")8 관련 주제
선행 지식
상위 주제
기존 세부 포스트
후속 주제
- 추정량 평가 방법 — MSE, Cramer-Rao, Rao-Blackwell, UMVUE
- 점근적 평가 (Ch.10) — MLE의 점근적 성질
9 참고 문헌
- Casella, G. & Berger, R. L. (2002). Statistical Inference (2nd ed.). Duxbury. Chapter 7, Section 7.2.
- Dempster, A. P., Laird, N. M. & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS-B, 39, 1-38.
- Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2, 110-114.
- Olkin, I., Petkau, A. J. & Zidek, J. V. (1981). A comparison of \(n\) estimators for the binomial distribution. JASA, 76, 637-642.
- Zehna, P. W. (1966). Invariance of maximum likelihood estimators. Annals of Mathematical Statistics, 37, 744.