확률표는 어디에서 오는가

M18과 M19에서는 확률표(π, a, b)를 임의로 주고 시작했습니다. 실제로는 이 표가 처음부터 있지는 않습니다. 관측 서열만 있고 상태 라벨은 없는 상황이 대부분입니다. 데이터에서 이 표를 어떻게 배울 수 있을까요?

정답이 있는 지도학습이라면 카운팅 한 번으로 끝납니다. 각 상태에서 각 문자가 몇 번 나왔는지 세면 됩니다. 그런데 상태를 볼 수 없다면? Baum-Welch 알고리즘이 이 문제를 다룹니다. Expectation-Maximization(EM)의 특수 사례로, 상태 배열을 확률적으로 추정하면서 파라미터를 반복 재추정합니다. 이 편은 HMM 시리즈의 자가 학습 심장입니다.

EM의 두 단계 — 감으로 잡기

EM은 이렇게 굴러갑니다.

1. 초기화: 파라미터 표에 아무 값이나 넣는다 (완전 랜덤 or 균등)
2. E-step (Expectation): 현재 파라미터로 각 시각·상태의 확률(γ)과 전이 확률(ξ)을 계산
3. M-step (Maximization): 그 확률을 카운트인 것처럼 취급해서 파라미터를 재추정
4. E와 M을 log likelihood가 수렴할 때까지 반복

핵심 통찰은 이겁니다. 상태를 확실히 모르니까 확률로 처리하고, 그 확률을 마치 진짜 카운트인 것처럼 써서 파라미터를 다시 계산한다. 자기 자신을 조금 더 맞게 만드는 부트스트래핑입니다.

E-step — Forward-Backward를 재료로 두 값을 뽑기

M19에서 계산한 값 두 개가 여기서 재료가 됩니다.

γ_t(s) = 시각 t에서 상태 s일 posterior 확률

$\gamma_t(s) = \frac{\alpha_t(s) \beta_t(s)}{P(O)}$

ξ_t(i, j) = 시각 t에서 상태 i, t+1에서 상태 j로 전이할 posterior 확률 (신규)

$\xi_t(i, j) = \frac{\alpha_t(i) \cdot a(i \to j) \cdot b(j, o_{t+1}) \cdot \beta_{t+1}(j)}{P(O)}$

ξ는 γ의 "전이" 버전입니다. 한 상태에서 다음 상태로 넘어가는 확률을 두 시각의 α와 β로 조립합니다. 이 두 값이 있으면 M-step은 카운팅과 같은 방식으로 됩니다.

M-step — 카운트 재추정

초기 확률 재추정:

$\pi'(s) = \gamma_1(s)$

첫 시각의 posterior가 곧 새 초기 확률입니다.

전이 확률 재추정:

$a'(i \to j) = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)}$

분자는 "i에서 j로 넘어간 기대 횟수", 분모는 "i에서 어디로든 넘어간 기대 횟수". 두 카운트의 비가 조건부 확률의 정의 그대로입니다.

방출 확률 재추정:

$b'(s, v_k) = \frac{\sum_{t: o_t = v_k} \gamma_t(s)}{\sum_{t=1}^{T} \gamma_t(s)}$

분자는 "상태 s에서 문자 v_k가 관측된 기대 횟수", 분모는 "상태 s가 나온 총 기대 횟수". 카운팅과 완벽히 같은 구조입니다.

이 세 식이 M-step의 전부입니다. E-step에서 뽑은 γ와 ξ를 카운트처럼 다뤄서 파라미터의 최우도 추정을 재계산할 뿐입니다.

왜 이게 수렴하나

EM에는 강력한 보장이 있습니다. 매 반복마다 데이터의 log likelihood가 단조 증가한다. 절대 나빠지지 않습니다. 유한한 반복 후 수렴합니다.

증명은 Jensen 부등식과 Q-함수 유도가 필요해서 이 편의 범위를 넘습니다. 하지만 감으로 잡을 수는 있습니다. E-step은 현재 파라미터가 만든 posterior를 사용하고, M-step은 그 posterior 하에서 최적 파라미터를 뽑습니다. 매 단계가 "고정된 다른 쪽 하에서 자기 쪽만 최적화"하는 좌표하강법의 확률 버전입니다.

다만 조심할 점이 있습니다. EM은 지역 최적에 수렴할 뿐 전역 최적을 보장하지 않습니다. 초기값이 나쁘면 나쁜 지역 최적에 갇힙니다. 실무에서는 다른 초기값으로 여러 번 돌려서 가장 좋은 결과를 고르거나, K-means 같은 사전 클러스터링으로 초기값을 잡아줍니다.

손으로 한 반복 돌려봅시다

M19에서 이미 계산한 α, β, γ를 그대로 씁니다. 관측 GCAT, 초기 확률표 M18과 같음.

M19의 posterior (γ):

(t=1, 2의 γ는 α×β/P(O) 공식으로 손 계산 가능; 검증 연습으로 남깁니다.)

ξ 계산 예 (t=1, H→H): $\xi_1(H, H) = \frac{\alpha_1(H) \cdot a(H\to H) \cdot b(H, C) \cdot \beta_2(H)}{P(O)}$ $= \frac{0.20 \times 0.7 \times 0.4 \times 0.0469}{0.003676} = \frac{0.002627}{0.003676} = 0.715$

이런 식으로 (i, j) 조합 4개 × 시각 3개 = 12개의 ξ가 나옵니다. 전이 재추정에서

$a'(H \to H) = \frac{\xi_1(H,H) + \xi_2(H,H) + \xi_3(H,H)}{\gamma_1(H) + \gamma_2(H) + \gamma_3(H)}$

이 계산을 마치면 새 전이 확률표가 나옵니다. 기존 0.7보다 아마 살짝 다른 값. 이 새 표로 다음 E-step을 돌리고, 또 M-step. Log likelihood가 반올림 오차 이내로 안정될 때까지 반복합니다.

파이썬 구현 (수렴 루프 포함)

import math

def logsumexp(vals):
    m = max(vals)
    return m + math.log(sum(math.exp(v - m) for v in vals)) if m != float('-inf') else float('-inf')

def forward_backward(obs, N, start_p, trans_p, emit_p):
    T = len(obs)
    log = math.log

    alpha = [[log(start_p[s]) + log(emit_p[s][obs[0]]) if t == 0 else 0.0
              for s in range(N)] for t in range(T)]
    for t in range(1, T):
        for s in range(N):
            terms = [alpha[t-1][sp] + log(trans_p[sp][s]) for sp in range(N)]
            alpha[t][s] = logsumexp(terms) + log(emit_p[s][obs[t]])
    log_P_O = logsumexp(alpha[T-1])

    beta = [[0.0] * N for _ in range(T)]
    for t in range(T-2, -1, -1):
        for s in range(N):
            terms = [log(trans_p[s][sp]) + log(emit_p[sp][obs[t+1]]) + beta[t+1][sp]
                     for sp in range(N)]
            beta[t][s] = logsumexp(terms)

    return alpha, beta, log_P_O

def baum_welch(obs, N, M, start_p, trans_p, emit_p, max_iter=100, tol=1e-5):
    T = len(obs)
    prev_ll = float('-inf')

    for iteration in range(max_iter):
        alpha, beta, log_P_O = forward_backward(obs, N, start_p, trans_p, emit_p)

        # E-step
        gamma = [[math.exp(alpha[t][s] + beta[t][s] - log_P_O) for s in range(N)]
                 for t in range(T)]
        xi = [[[0.0] * N for _ in range(N)] for _ in range(T-1)]
        for t in range(T-1):
            for i in range(N):
                for j in range(N):
                    num = alpha[t][i] + math.log(trans_p[i][j]) + \
                          math.log(emit_p[j][obs[t+1]]) + beta[t+1][j]
                    xi[t][i][j] = math.exp(num - log_P_O)

        # M-step
        start_p = [gamma[0][s] for s in range(N)]
        for i in range(N):
            denom = sum(gamma[t][i] for t in range(T-1))
            for j in range(N):
                trans_p[i][j] = sum(xi[t][i][j] for t in range(T-1)) / denom
        for s in range(N):
            denom = sum(gamma[t][s] for t in range(T))
            for k in range(M):
                num = sum(gamma[t][s] for t in range(T) if obs[t] == k)
                emit_p[s][k] = num / denom

        # 수렴 판정
        if abs(log_P_O - prev_ll) < tol:
            break
        prev_ll = log_P_O

    return start_p, trans_p, emit_p, log_P_O, iteration + 1

# 랜덤 초기값으로 시작
import random
random.seed(42)
N, M = 2, 4
start_p = [0.6, 0.4]
trans_p = [[0.5, 0.5], [0.5, 0.5]]
emit_p = [[random.random() for _ in range(M)] for _ in range(N)]
for s in range(N):
    total = sum(emit_p[s])
    emit_p[s] = [v/total for v in emit_p[s]]

# 관측 (여러 개 이어붙임)
obs = [{'A':0,'C':1,'G':2,'T':3}[c] for c in 'GCATGCATGGCCAATT']
sp, tp, ep, ll, iters = baum_welch(obs, N, M, start_p, trans_p, emit_p)

print(f'수렴 반복 수: {iters}, 최종 log P(O): {ll:.4f}')
print('학습된 전이:', tp)
print('학습된 방출:', ep)

실행하면 몇 십 번 반복 만에 수렴하고, GC-rich vs AT-rich에 가까운 방출 확률표가 등장합니다. 정답 라벨을 준 적이 없는데 알고리즘이 문자 통계만으로 두 개의 잠재 상태를 발견해냅니다. 이게 unsupervised learning의 원형입니다.

지역 최적의 함정 — 반드시 알아둘 것

Baum-Welch로 여러 번 실행해보면 결과가 매번 다릅니다.

• 초기값을 이상하게 잡으면 두 상태가 거의 동일한 방출 분포로 수렴합니다 (상태가 하나로 붕괴)
• 초기값이 우연히 데이터 구조와 정반대면 나쁜 지역 최적에 갇힙니다
• Log likelihood가 낮아도 알고리즘은 만족하며 멈춥니다

실무에서는 20～100개의 랜덤 시작점으로 반복 실행하고 가장 높은 log likelihood를 뽑는 것이 표준입니다. 또는 K-means로 관측을 사전 클러스터링해서 초기 방출 확률표를 그럴싸하게 잡습니다. 이런 트릭들이 없으면 Baum-Welch는 실전에서 안정적이지 못합니다.

복잡도

한 번의 EM 반복 = O(T × N²) (Forward-Backward와 같음). 수렴까지 보통 수십～수백 반복. 총 O(T × N² × I). T가 수만이어도 실용적이지만, 서열이 여러 개일 때는 배치 병렬화가 자연스럽습니다.

CS 매핑 — EM의 정수

DryBench에서 언급된 좌표하강법과 부트스트래핑의 확률 버전입니다.

• 좌표하강: 파라미터와 잠재변수를 번갈아가며 최적화
• 부트스트래핑: 자기가 만든 posterior를 카운트로 재사용
• 단조 수렴: log likelihood는 절대 감소하지 않음
• 지역 최적: 초기값에 종속. 여러 재시작 필요

이 알고리즘의 아이디어는 다른 곳에 반복 등장합니다. GMM(Gaussian Mixture Model)의 EM, LDA(토픽 모델)의 variational EM, 반지도학습의 pseudo-labeling. 관측된 부분에서 잠재된 부분을 추정하고, 그 추정을 카운트처럼 사용해서 파라미터를 다시 배운다 — 이 패턴 하나가 확률 모델 학습의 척추입니다.

Rosalind에서 채점받기

Rosalind BW 문제. 관측 서열과 초기 파라미터가 주어지고, N번 Baum-Welch 반복 후 최종 파라미터를 요구합니다. 위 코드의 for 루프 반복 수만 문제의 지정 값으로 바꾸면 됩니다.

다음 편으로 이어지는 갈래

• 다음 편 (M21): Profile HMM + gene finding — HMM을 실전 유전자 예측에 붙입니다. 상태가 두 개가 아니라 엑손·인트론·스플라이스 사이트 등 수십 개로 커집니다. HMMER 도구.
• 관련 확장: 온라인 EM — 전체 데이터를 매번 훑지 않고 스트리밍으로 파라미터 업데이트. 대용량 서열 학습에 필수.
• CRF 학습 — 로그 우도의 기울기를 직접 계산해서 파라미터 학습. HMM보다 유연하지만 계산이 무겁습니다.

더 깊게 파고 싶다면

• Baum, Petrie, Soules, Weiss (1970), A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains, Annals of Mathematical Statistics — 원 논문. 반세기가 지났지만 정수는 살아있습니다.
• Dempster, Laird, Rubin (1977), Maximum Likelihood from Incomplete Data via the EM Algorithm, JRSS-B — EM을 일반 프레임워크로 정립한 논문.
• UC Berkeley CS176 — Yun S. Song 교수의 HMM 학습 강의. 판서로 Q-함수 유도를 봅니다.
• Bishop Pattern Recognition and Machine Learning Chapter 9 — GMM과 HMM EM을 하나의 프레임으로 통합해서 유도.
• Durbin et al. Biological Sequence Analysis Chapter 3.3 — 바이오 문맥의 Baum-Welch 튜토리얼.

다음 편으로 넘어가기 전에, 위 파이썬 코드를 실행해서 초기값을 여러 번 바꿔봅시다. 어떤 초기값이 좋은 지역 최적으로 수렴하고, 어떤 초기값이 나쁜 곳에 갇히는지 감을 잡는 것이 이 편의 진짜 실습입니다.