EM 알고리즘 (The EM Algorithm)

1 개요

추정량 탐색 방법에서 EM 알고리즘을 네 가지 방법 중 하나로 요약했다. 이 포스트에서는 EM(Expectation-Maximization) 알고리즘을 심층적으로 다룬다 (Casella & Berger, 2002, Ch.7.2.4; Dempster, Laird, & Rubin, 1977).

EM 알고리즘은 다른 추정 방법과 근본적으로 다르다. 적률법이나 MLE가 “추정량을 구하는 공식”을 제공하는 반면, EM은 MLE에 수렴하는 반복 알고리즘을 제공한다. 핵심 아이디어는 단순하다:

하나의 어려운 우도 최대화를 일련의 쉬운 최대화로 대체한다.

이것은 관측되지 않은 잠재 데이터(latent data) 가 있을 때 특히 유용하다. 잠재 데이터가 관측되었다면 완전 데이터 우도는 쉽게 최대화할 수 있다. 잠재 데이터가 없으므로, 현재 추정값 하에서 잠재 데이터의 기댓값을 대입하고(E-step), 이를 기반으로 최대화한다(M-step).

---

2 완전 데이터와 불완전 데이터

2.1 두 가지 우도 문제

EM 알고리즘에서는 두 가지 우도 문제를 구분한다:

	완전 데이터	불완전 데이터
관측	\((\mathbf{Y}, \mathbf{X})\) — 관측 + 잠재	\(\mathbf{Y}\) — 관측만
우도	\(L(\theta|\mathbf{y}, \mathbf{x}) = f(\mathbf{y}, \mathbf{x}|\theta)\)	\(L(\theta|\mathbf{y}) = g(\mathbf{y}|\theta)\)
최대화	쉬움 (닫힌 형태 해가 자주 존재)	어려움 (적분/합산 필요)

두 우도의 관계:

\[ g(\mathbf{y}|\theta) = \int f(\mathbf{y}, \mathbf{x}|\theta) \, d\mathbf{x} \]

(이산 잠재 변수의 경우 적분을 합산으로 대체한다.)

우리의 목표는 \(L(\theta|\mathbf{y})\) 를 최대화하는 것이지만, \(L(\theta|\mathbf{y}, \mathbf{x})\) 를 통해 우회한다.

2.2 “결측 데이터”의 넓은 해석

“잠재 데이터”는 반드시 물리적으로 결측된 값일 필요가 없다. 계산을 편리하게 하기 위해 인위적으로 도입한 변수일 수도 있다:

상황	관측 데이터 \(\mathbf{Y}\)	잠재 데이터 \(\mathbf{X}\)
물리적 결측	관측된 값	결측된 관측값
혼합 모형	관측값 \(y_i\)	클러스터 소속 \(z_i\)
절단 데이터	관측된 생존 시간	절단된 사건 시간
Hidden Markov	관측 시퀀스	은닉 상태 시퀀스

---

3 EM 알고리즘의 형식적 정의

3.1 핵심 항등식

잠재 데이터의 조건부 분포를 \(k(\mathbf{x}|\theta, \mathbf{y}) = f(\mathbf{y}, \mathbf{x}|\theta) / g(\mathbf{y}|\theta)\) 로 정의하면

\[ \log L(\theta|\mathbf{y}) = \log L(\theta|\mathbf{y}, \mathbf{x}) - \log k(\mathbf{x}|\theta, \mathbf{y}) \]

양변에 \(E_{\mathbf{X}|\mathbf{y}, \theta'}[\cdot]\) (현재 추정값 \(\theta'\) 하에서의 조건부 기댓값)을 취하면

\[ \log L(\theta|\mathbf{y}) = \underbrace{E[\log L(\theta|\mathbf{y}, \mathbf{X})|\theta', \mathbf{y}]}_{Q(\theta|\theta')} - \underbrace{E[\log k(\mathbf{X}|\theta, \mathbf{y})|\theta', \mathbf{y}]}_{H(\theta|\theta')} \]

좌변의 \(\log L(\theta|\mathbf{y})\) 는 \(\mathbf{X}\) 에 무관하므로 기댓값을 취해도 변하지 않는다.

3.2 E-step과 M-step

초기값 \(\theta^{(0)}\) 에서 출발하여, \(r = 0, 1, 2, \ldots\) 에 대해 반복한다:

E-step (Expectation Step)

현재 추정값 \(\theta^{(r)}\) 하에서, 완전 데이터 로그우도의 조건부 기댓값을 계산한다:

\[ Q(\theta|\theta^{(r)}) = E_{\mathbf{X}|\mathbf{y}, \theta^{(r)}}[\log L(\theta|\mathbf{y}, \mathbf{X})] \]

M-step (Maximization Step)

\(Q\) 함수를 \(\theta\) 에 대해 최대화한다:

\[ \theta^{(r+1)} = \arg\max_\theta Q(\theta|\theta^{(r)}) \]

직관: E-step에서 “잠재 데이터의 자리에 기댓값을 대입”하고, M-step에서 “마치 완전 데이터가 관측된 것처럼 MLE를 구한다.”

---

4 왜 작동하는가: 단조 수렴 정리

4.1 정리

정리 7.2.20: EM 단조 수렴

EM 반복 수열 \(\{\theta^{(r)}\}\) 는 다음을 만족한다:

\[ L(\theta^{(r+1)}|\mathbf{y}) \geq L(\theta^{(r)}|\mathbf{y}) \]

등호는 연속된 반복이 동일한 \(Q\) 최대값을 산출할 때만 성립한다.

4.2 증명 스케치

핵심 항등식 \(\log L(\theta|\mathbf{y}) = Q(\theta|\theta^{(r)}) - H(\theta|\theta^{(r)})\) 에서

\[ \log L(\theta^{(r+1)}|\mathbf{y}) - \log L(\theta^{(r)}|\mathbf{y}) = [Q(\theta^{(r+1)}|\theta^{(r)}) - Q(\theta^{(r)}|\theta^{(r)})] - [H(\theta^{(r+1)}|\theta^{(r)}) - H(\theta^{(r)}|\theta^{(r)})] \]

첫 번째 차이 \(\geq 0\): M-step의 정의에 의해 \(Q(\theta^{(r+1)}|\theta^{(r)}) \geq Q(\theta^{(r)}|\theta^{(r)})\)

두 번째 차이 \(\leq 0\): Jensen 부등식(또는 Gibbs 부등식/KL 발산의 비음수성)에 의해

\[ H(\theta|\theta') = E_{\theta'}[\log k(\mathbf{X}|\theta, \mathbf{y})] \leq H(\theta'|\theta') \]

이것은 \(\text{KL}(k(\cdot|\theta', \mathbf{y}) \| k(\cdot|\theta, \mathbf{y})) \geq 0\) 의 결과이다.

따라서 \(\log L(\theta^{(r+1)}|\mathbf{y}) - \log L(\theta^{(r)}|\mathbf{y}) \geq 0 - 0 = 0\). \(\square\)

해석: 불완전 데이터 우도가 매 반복마다 증가하거나 유지된다. 우도가 위로 유계이면 수열은 수렴한다.

---

5 예시 1: 다중 포아송 비율 (결측 데이터)

5.1 설정

\(X_i \sim \text{Poisson}(\tau_i)\) 와 \(Y_i \sim \text{Poisson}(\beta\tau_i)\) 를 독립으로 관측한다 (\(i = 1, \ldots, n\)). \(Y_i\) 는 지역 \(i\) 의 질병 발생 건수이고, \(\tau_i\) 는 기저 비율, \(\beta\) 는 전체 효과이다.

완전 데이터 우도: \(((\mathbf{x}, \mathbf{y}))\) 를 모두 관측했으면

\[ \hat{\beta} = \frac{\sum y_i}{\sum x_i}, \qquad \hat{\tau}_j = \frac{x_j + y_j}{\hat{\beta} + 1} \]

불완전 데이터: \(x_1\) 이 결측되었다면, 불완전 데이터 우도 \(L(\beta, \tau_1, \ldots, \tau_n | y_1, (x_2, y_2), \ldots, (x_n, y_n))\) 을 직접 최대화해야 한다.

5.2 EM 적용

E-step: \(x_1\) 의 조건부 기대를 \(\theta^{(r)}\) 하에서 계산한다. \(X_1 \sim \text{Poisson}(\tau_1)\) 이므로 \(E[X_1|\theta^{(r)}] = \tau_1^{(r)}\)

완전 데이터 로그우도에서 \(x_1\) 을 \(\tau_1^{(r)}\) 로 대체한다.

M-step: 대체된 완전 데이터 MLE를 계산한다.

\[ \hat{\beta}^{(r+1)} = \frac{\sum_{i=1}^n y_i}{\tau_1^{(r)} + \sum_{i=2}^n x_i}, \quad \hat{\tau}_1^{(r+1)} = \frac{\tau_1^{(r)} + y_1}{\hat{\beta}^{(r+1)} + 1} \]

\[ \hat{\tau}_j^{(r+1)} = \frac{x_j + y_j}{\hat{\beta}^{(r+1)} + 1}, \quad j = 2, \ldots, n \]

이 수열은 불완전 데이터 MLE에 수렴한다.

---

6 예시 2: 가우시안 혼합 모형 (Gaussian Mixture Model)

EM의 가장 대표적인 응용이다.

6.1 설정

\(K\) 개의 정규분포 혼합에서 \(n\) 개의 관측값 \(y_1, \ldots, y_n\) 이 iid로 추출되었다:

\[ f(y|\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\sigma}^2) = \sum_{k=1}^K \pi_k \, \phi(y|\mu_k, \sigma_k^2) \]

• \(\pi_k\): 혼합 비율 (\(\sum \pi_k = 1\), \(\pi_k > 0\))
• \(\mu_k, \sigma_k^2\): \(k\) 번째 성분의 평균, 분산
• \(\phi(\cdot|\mu, \sigma^2)\): 정규 pdf

6.2 잠재 변수

\(Z_i \in \{1, \ldots, K\}\): \(y_i\) 가 어느 성분에서 나왔는지를 나타내는 잠재 변수. \(P(Z_i = k) = \pi_k\)

완전 데이터: \((y_i, z_i)\) 를 관측했다면, 성분별로 데이터를 분리하여 각 성분의 평균/분산을 독립적으로 추정할 수 있다.

6.3 E-step 유도

현재 추정값 \(\theta^{(t)} = (\boldsymbol{\pi}^{(t)}, \boldsymbol{\mu}^{(t)}, (\boldsymbol{\sigma}^2)^{(t)})\) 하에서, \(Z_i\) 의 사후 확률(responsibility)을 계산한다:

\[ \gamma_{ik}^{(t)} = P(Z_i = k | y_i, \theta^{(t)}) = \frac{\pi_k^{(t)} \, \phi(y_i|\mu_k^{(t)}, (\sigma_k^2)^{(t)})}{\sum_{j=1}^K \pi_j^{(t)} \, \phi(y_i|\mu_j^{(t)}, (\sigma_j^2)^{(t)})} \]

\(\gamma_{ik}^{(t)}\) 는 “\(y_i\) 가 성분 \(k\) 에서 나왔을 확률”에 대한 현재 최선의 추정이다.

6.4 M-step 유도

\(\gamma_{ik}^{(t)}\) 를 가중치로 사용하여 가중 MLE를 계산한다:

\[ N_k = \sum_{i=1}^n \gamma_{ik}^{(t)} \quad \text{(성분 $k$ 의 유효 표본 크기)} \]

\[ \mu_k^{(t+1)} = \frac{1}{N_k} \sum_{i=1}^n \gamma_{ik}^{(t)} y_i \]

\[ (\sigma_k^2)^{(t+1)} = \frac{1}{N_k} \sum_{i=1}^n \gamma_{ik}^{(t)} (y_i - \mu_k^{(t+1)})^2 \]

\[ \pi_k^{(t+1)} = \frac{N_k}{n} \]

직관: 각 관측값을 “부분적으로” 각 성분에 할당한다. \(\gamma_{ik}\) 가 할당 비중이고, 이 비중으로 가중 평균/분산을 계산한다.

6.5 수렴 판단

실무적 수렴 기준:

\[ |\ell(\theta^{(t+1)}|\mathbf{y}) - \ell(\theta^{(t)}|\mathbf{y})| < \epsilon \quad \text{또는} \quad \|\theta^{(t+1)} - \theta^{(t)}\| < \epsilon \]

\(\epsilon\) 은 보통 \(10^{-6}\) ~ \(10^{-8}\) 수준이다.

---

7 EM의 실무적 고려사항

7.1 초기값 민감성

EM은 국소 최대에 수렴할 수 있다. 초기값에 따라 다른 국소 최대에 도달할 수 있으므로:

전략	설명
다중 시작(multiple starts)	여러 랜덤 초기값에서 실행, 최대 우도 선택
K-means 초기화	K-means 결과를 EM의 초기값으로 사용
적률법 초기화	데이터의 적률로 초기 모수 추정
계층적 초기화	작은 \(K\) 에서 시작하여 점진적으로 증가

7.2 수렴 속도

EM의 수렴 속도는 일반적으로 선형(linear) 이다. 뉴턴-랩슨의 이차(quadratic) 수렴보다 느리다. 수렴 속도는 불완전 정보의 비율에 의존한다:

\[ \text{수렴 속도} \approx 1 - \frac{I_{\text{obs}}(\theta)}{I_{\text{comp}}(\theta)} \]

관측 정보가 완전 정보에 비해 적을수록 (잠재 데이터의 정보 비중이 클수록) 수렴이 느리다.

7.3 변형과 확장

변형	핵심 아이디어
ECM (Expectation Conditional Maximization)	M-step을 조건부 최대화로 분할
ECME	ECM에서 일부 M-step을 실제 우도로 최대화
PX-EM (Parameter-Expanded EM)	확장된 모수 공간에서 더 빠른 수렴
Monte Carlo EM	E-step의 기댓값을 몬테카를로로 근사
Variational EM	E-step을 변분 근사로 대체
온라인 EM	미니배치로 순차 업데이트

---

8 EM과 다른 방법의 비교

기준	EM	뉴턴-랩슨	경사하강법
수렴 보장	단조 증가 보장	보장 없음	보장 없음 (볼록이면 보장)
수렴 속도	선형	이차	선형 (스텝 크기 의존)
구현 난이도	중간 (해석적 E/M)	높음 (헤시안 필요)	낮음
잠재 변수	자연스럽게 처리	직접 처리 어려움	주변화 필요
국소 최대 위험	있음	있음	있음
스텝 크기 조정	불필요	불필요 (근방)	필요

---

9 코드 예시

9.1 Step 1: 순수 Python 구현 (1차원 GMM EM)

2-성분 가우시안 혼합 모형의 EM을 순수 Python으로 구현한다.

import math
import random

random.seed(42)

# 데이터 생성: 2-성분 혼합 N(2,1) + N(7,1.5)
data = []
for _ in range(200):
    if random.random() < 0.4:
        data.append(random.gauss(2.0, 1.0))
    else:
        data.append(random.gauss(7.0, 1.5))

n = len(data)
K = 2

# 초기값
mu = [1.0, 8.0]
sigma2 = [2.0, 2.0]
pi_k = [0.5, 0.5]

def normal_pdf(x, mu, sigma2):
    return (1 / math.sqrt(2 * math.pi * sigma2)) * \
           math.exp(-0.5 * (x - mu)**2 / sigma2)

def log_likelihood(data, mu, sigma2, pi_k, K):
    ll = 0
    for x in data:
        p = sum(pi_k[k] * normal_pdf(x, mu[k], sigma2[k]) for k in range(K))
        ll += math.log(max(p, 1e-300))
    return ll

print("=== 2-성분 GMM: EM 알고리즘 ===")
print(f"참값: mu=[2.0, 7.0], sigma=[1.0, 1.5], pi=[0.4, 0.6]\n")

for iteration in range(100):
    ll_old = log_likelihood(data, mu, sigma2, pi_k, K)

    # === E-step ===
    gamma = []
    for i in range(n):
        row = [pi_k[k] * normal_pdf(data[i], mu[k], sigma2[k]) for k in range(K)]
        total = sum(row)
        gamma.append([r / total for r in row])

    # === M-step ===
    for k in range(K):
        N_k = sum(gamma[i][k] for i in range(n))
        mu[k] = sum(gamma[i][k] * data[i] for i in range(n)) / N_k
        sigma2[k] = sum(gamma[i][k] * (data[i] - mu[k])**2
                        for i in range(n)) / N_k
        pi_k[k] = N_k / n

    ll_new = log_likelihood(data, mu, sigma2, pi_k, K)

    if iteration < 5 or iteration % 10 == 9:
        print(f"  반복 {iteration+1:3d}: log L = {ll_new:.4f}, "
              f"mu = [{mu[0]:.3f}, {mu[1]:.3f}], "
              f"pi = [{pi_k[0]:.3f}, {pi_k[1]:.3f}]")

    if abs(ll_new - ll_old) < 1e-10:
        print(f"\n  수렴: {iteration+1}회 반복")
        break

print(f"\n  최종 결과:")
for k in range(K):
    print(f"    성분 {k+1}: mu={mu[k]:.4f}, "
          f"sigma={math.sqrt(sigma2[k]):.4f}, pi={pi_k[k]:.4f}")

9.2 Step 2: scipy/sklearn 구현 (GMM + 초기값 비교)

다중 초기값에서의 EM 결과를 비교하여 국소 최대 문제를 시연한다.

import numpy as np
from scipy.stats import norm

np.random.seed(42)

# 3-성분 혼합 모형 데이터 생성
n_samples = 500
mu_true = [-3.0, 2.0, 7.0]
sigma_true = [1.0, 0.8, 1.2]
pi_true = [0.3, 0.4, 0.3]

# 데이터 생성
z = np.random.choice(3, size=n_samples, p=pi_true)
data = np.array([np.random.normal(mu_true[zi], sigma_true[zi]) for zi in z])

K = 3

def em_gmm(data, K, mu_init, max_iter=200, tol=1e-8):
    """EM for Gaussian Mixture Model"""
    n = len(data)
    mu = np.array(mu_init, dtype=float)
    sigma2 = np.ones(K) * np.var(data) / K
    pi_k = np.ones(K) / K

    log_liks = []

    for iteration in range(max_iter):
        # 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 = (gamma.T @ data) / N_k
        for k in range(K):
            sigma2[k] = np.sum(gamma[:, k] * (data - mu[k])**2) / N_k[k]
            sigma2[k] = max(sigma2[k], 1e-6)  # 분산 하한
        pi_k = N_k / n

        # 로그우도
        ll = np.sum(np.log(sum(pi_k[k] * norm.pdf(data, mu[k], np.sqrt(sigma2[k]))
                               for k in range(K))))
        log_liks.append(ll)

        if len(log_liks) > 1 and abs(log_liks[-1] - log_liks[-2]) < tol:
            break

    return mu, np.sqrt(sigma2), pi_k, log_liks

# 다양한 초기값으로 실행
print("=== 3-성분 GMM: 초기값 민감도 분석 ===")
print(f"참값: mu={mu_true}, sigma={sigma_true}, pi={pi_true}\n")

inits = [
    ("좋은 초기값", [-2.0, 1.0, 6.0]),
    ("랜덤 초기값 1", [0.0, 0.0, 0.0]),
    ("랜덤 초기값 2", [-5.0, -5.0, 10.0]),
    ("K-means 근사", [data[data < 0].mean(),
                     data[(data > 0) & (data < 5)].mean(),
                     data[data > 5].mean()]),
]

best_ll = -np.inf
best_result = None

for name, mu_init in inits:
    mu_est, sigma_est, pi_est, lls = em_gmm(data, K, mu_init)
    final_ll = lls[-1]
    converged_in = len(lls)

    if final_ll > best_ll:
        best_ll = final_ll
        best_result = (mu_est, sigma_est, pi_est)

    # 정렬 (mu 기준)
    order = np.argsort(mu_est)
    print(f"  {name:15s}: {converged_in:3d}회, log L = {final_ll:.2f}")
    for k in order:
        print(f"    mu={mu_est[k]:7.3f}, sigma={sigma_est[k]:.3f}, pi={pi_est[k]:.3f}")
    print()

print(f"  최적 결과 (최대 로그우도: {best_ll:.2f}):")
mu_b, sigma_b, pi_b = best_result
order = np.argsort(mu_b)
for k in order:
    print(f"    mu={mu_b[k]:7.3f}, sigma={sigma_b[k]:.3f}, pi={pi_b[k]:.3f}")

# 단조 수렴 확인
_, _, _, lls_check = em_gmm(data, K, [-2.0, 1.0, 6.0])
diffs = np.diff(lls_check)
print(f"\n  단조 수렴 확인: 모든 ll 차이 >= 0? {np.all(diffs >= -1e-10)}")
print(f"  최소 차이: {diffs.min():.2e}, 최대 차이: {diffs.max():.2e}")

---

10 응용 분야

분야	EM 활용	잠재 변수
클러스터링	가우시안 혼합 모형	클러스터 소속
결측치 대체	다중 대체(MI)의 모수 추정	결측된 값
자연어처리	HMM 기반 품사 태깅, 토픽 모델	은닉 상태, 토픽 할당
컴퓨터 비전	이미지 분할, 배경/전경 분리	픽셀의 클래스
유전체학	인구 구조 분석 (STRUCTURE)	조상 소속
신호 처리	소스 분리 (ICA 변형)	소스 신호
생존 분석	절단/중도절단 데이터의 MLE	절단된 사건 시간

---

11 관련 주제

선행 지식

• 최대우도추정량 — EM이 수렴하는 목표
• 우도원리 — 우도함수의 본질
• 충분성 원리 — 완전 데이터의 충분통계량 활용

상위 주제

• 점추정 개요
• 추정량 탐색 방법

후속 주제

• 베이즈 추정량 — 베이지안 EM, 변분 EM
• Monte Carlo 시뮬레이션 — Monte Carlo EM

---

12 참고 문헌

• Casella, G. & Berger, R. L. (2002). Statistical Inference (2nd ed.). Duxbury. Chapter 7, Section 7.2.4.
• Dempster, A. P., Laird, N. M. & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS-B, 39, 1-38.
• McLachlan, G. J. & Krishnan, T. (2008). The EM Algorithm and Extensions (2nd ed.). Wiley.
• Wu, C. F. J. (1983). On the convergence properties of the EM algorithm. Annals of Statistics, 11, 95-103.
• Meng, X. L. & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm. Biometrika, 80, 267-278.