왜 HMM이 필요한가

앞서 배운 정렬은 두 서열의 문자를 직접 대응시켰습니다. 그런데 실제 서열에는 직접 관측할 수 없는 상태가 숨어 있습니다. 이 자리는 엑손인가 인트론인가? 이 영역은 GC-rich인가 AT-rich인가? 이 위치는 CpG 섬 안인가 밖인가? 이런 질문에 답하려면 문자 자체가 아니라 문자를 만들어낸 상태의 배열을 뽑아내야 합니다.

숨겨진 마르코프 모델(Hidden Markov Model, HMM)이 정확히 이 문제를 다룹니다. 서열 하나가 관측되면 그 밑에 깔린 상태 배열을 확률 모델로 되짚어 올라갑니다. 이 편에서는 HMM의 세 가지 확률을 정의하고, 그 상태 배열을 찾는 Viterbi 알고리즘을 격자로 유도합니다. Needleman-Wunsch의 격자와 놀랍도록 닮았습니다.

은닉 마르코프 모델의 세 가지 확률

HMM은 다음 다섯 가지 조각으로 정의됩니다.

• 상태 집합 S — 예: GC-rich(H), AT-rich(L)
• 관측 집합 V — 예: DNA 알파벳 A, C, G, T
• 초기 확률 pi — 시퀀스 시작 시 각 상태의 확률
• 전이 확률 a — 상태 i에서 j로 넘어갈 확률
• 방출 확률 b — 상태 s일 때 문자 v가 나올 확률

핵심 가정 두 가지가 있습니다.

1. 마르코프 성질: 다음 상태는 오직 현재 상태에만 의존 (과거 이력 무의미)
2. 관측 독립성: 관측 문자는 오직 현재 상태에만 의존

이 두 가정 덕분에 상태 배열의 결합 확률을 곱셈으로 쪼갤 수 있습니다. 관측 배열과 상태 배열의 결합 확률은 초기 확률 곱하기 첫 방출 곱하기 (전이 곱하기 방출)을 시각 T까지 곱한 것과 같습니다.

간단한 예제 — GC-rich vs AT-rich 검출기

박테리아 게놈에서 CpG 섬을 검출하는 문제를 축약해봅시다. 상태 두 개, 관측 알파벳 네 자리.

• 상태: H (GC-rich, high) · L (AT-rich, low)
• 관측: A, C, G, T

임의로 세팅한 확률표 (실제로는 데이터에서 추정합니다. M20 Baum-Welch 편에서 다룹니다).

초기 확률:

전이 확률 a:

H는 H를 유지하는 관성이 강하고, L도 L을 유지하는 성향이 있습니다. 상태 전환은 그때그때 페널티가 붙는 셈입니다.

방출 확률 b:

H 상태는 C와 G를 더 많이 뱉고, L 상태는 A와 T를 더 많이 뱉습니다. 실제 GC-rich 영역의 통계와 정성적으로 일치합니다.

관측 서열을 GCAT이라고 합시다. 이 서열을 만들 수 있는 상태 배열은 2의 4승 즉 16개. 그중 확률이 가장 큰 배열이 무엇일까요?

Viterbi 알고리즘 — 최대 확률 상태 배열

Viterbi가 던지는 질문은 이렇습니다. 관측 배열이 주어졌을 때, 그것을 만들어낸 상태 배열의 결합 확률이 최대가 되는 배열은?

가능한 상태 배열은 상태 수의 서열 길이 제곱만큼으로 지수적입니다. 정공법은 불가능. 대신 격자 위 최적 경로 문제로 다시 씁니다.

격자 정의: 행은 시간 t, 열은 상태 s. 셀 (t, s)에는 시각 t에 상태 s로 끝나는 최적 부분 경로의 확률을 저장합니다.

점화식 (자연어 서술):

• V(t, s) = 이전 시각의 모든 상태 u에 대해 max로 [ V(t-1, u) 곱하기 a(u, s) ] 그리고 그것에 b(s, o_t) 곱하기
• 경계 조건: V(1, s) = pi(s) 곱하기 b(s, o_1)

이 점화식은 Needleman-Wunsch와 판박이입니다. 차이는 두 가지뿐입니다.

1. 덧셈이 아니라 곱셈
2. 격자 셀 값이 정렬 점수가 아니라 부분 상태 배열의 확률

DP가 성립하는 이유도 완전히 같습니다. 어떤 시각 t의 상태 s로 끝나는 최적 배열은, t-1 시각의 어떤 상태 u로 끝나는 최적 배열 뒤에 (u에서 s로) 전이와 (s에서 문자 o_t로) 방출이 붙은 형태입니다. 최적 부분 구조.

손으로 격자를 채워봅시다

관측 GCAT, 앞의 확률표로 격자를 채웁니다.

t=1, 관측 G:

• V(1, H) = 0.5 곱 0.4 = 0.20
• V(1, L) = 0.5 곱 0.1 = 0.05

t=2, 관측 C:

• V(2, H) = max(0.20 곱 0.7, 0.05 곱 0.4) 곱 0.4 = 0.14 곱 0.4 = 0.056
• V(2, L) = max(0.20 곱 0.3, 0.05 곱 0.6) 곱 0.1 = 0.06 곱 0.1 = 0.006

두 셀 모두 H에서 왔음을 기록합니다 (backpointer).

t=3, 관측 A:

• V(3, H) = max(0.056 곱 0.7, 0.006 곱 0.4) 곱 0.1 = 0.0392 곱 0.1 = 0.00392
• V(3, L) = max(0.056 곱 0.3, 0.006 곱 0.6) 곱 0.4 = 0.0168 곱 0.4 = 0.00672

여기서 처음으로 L 셀이 H 셀을 앞지릅니다. A는 L 상태에서 훨씬 잘 나오기 때문입니다.

t=4, 관측 T:

• V(4, H) = max(0.00392 곱 0.7, 0.00672 곱 0.4) 곱 0.1 = 0.002744 곱 0.1 = 0.0002744
• V(4, L) = max(0.00392 곱 0.3, 0.00672 곱 0.6) 곱 0.4 = 0.004032 곱 0.4 = 0.0016128

표로 정리:

최대 확률은 V(4, L) = 0.0016128. Backpointer를 되짚어가면 최적 상태 배열은 H H L L. 관측 GCAT의 앞 두 글자는 GC-rich 상태, 뒤 두 글자는 AT-rich 상태에서 만들어진 것으로 해석됩니다.

이게 Viterbi의 결과입니다. 관측만 보고 그 뒤에 숨은 상태의 최우도 배열을 격자 최적 경로로 뽑아냈습니다.

로그 스케일링 — 실무의 필수 트릭

위 예제는 서열 길이가 4라서 곱셈이 안전했습니다. 실제 서열은 수천에서 수만입니다. 확률을 계속 곱하면 값이 0에 수렴해서 언더플로우가 발생합니다.

해결은 로그 스케일. 곱셈을 덧셈으로 바꿉니다.

• log V(t, s) = 이전 상태 u에 대해 max로 [ log V(t-1, u) 더하기 log a(u, s) ] 그리고 거기에 log b(s, o_t) 더하기

log(0)은 마이너스 무한대 (파이썬에서 float('-inf'))로 처리합니다. 이 트릭은 HMM뿐 아니라 대부분의 확률 DP에 공통으로 적용됩니다.

파이썬 구현 (30줄)

import math

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

    V = [[float('-inf')] * N for _ in range(T)]
    back = [[0] * N for _ in range(T)]

    for s in range(N):
        V[0][s] = log(start_p[s]) + log(emit_p[s][obs[0]])

    for t in range(1, T):
        for s in range(N):
            best_prev, best_score = 0, float('-inf')
            for u in range(N):
                score = V[t-1][u] + log(trans_p[u][s])
                if score > best_score:
                    best_score, best_prev = score, u
            V[t][s] = best_score + log(emit_p[s][obs[t]])
            back[t][s] = best_prev

    path = [0] * T
    path[T-1] = max(range(N), key=lambda s: V[T-1][s])
    for t in range(T-2, -1, -1):
        path[t] = back[t+1][path[t+1]]

    return [states[s] for s in path], V[T-1][path[T-1]]

states = ['H', 'L']
obs_map = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
obs = [obs_map[c] for c in 'GCAT']
start_p = [0.5, 0.5]
trans_p = [[0.7, 0.3], [0.4, 0.6]]
emit_p  = [[0.1, 0.4, 0.4, 0.1],
           [0.4, 0.1, 0.1, 0.4]]

path, log_prob = viterbi(obs, states, start_p, trans_p, emit_p)
print(path, math.exp(log_prob))

30줄이면 끝입니다. 상태 개수가 10개, 서열 길이가 10,000이 되어도 이 격자의 크기는 100,000셀. 여전히 즉시 계산됩니다.

복잡도

• 시간: T 곱 N의 제곱. T는 서열 길이, N은 상태 수. 각 시각의 각 상태에서 이전 시각의 모든 상태를 훑기 때문
• 공간: T 곱 N. 격자 전체와 backpointer

Needleman-Wunsch와 정확히 같은 클래스입니다. HMM 상태 수가 알파벳 수 대신 들어간 셈.

CS 매핑 — 격자 DP의 확률판

DryBench에서 다룬 DP 격자 최적 경로 문제가 여기 그대로 재현됩니다.

• 격자 좌표 = (시각 t, 상태 s)
• 경로 = 하나의 상태 배열
• 점수 함수 = 로그 확률의 누적 합
• 최적 부분 구조 = t-1까지의 최적 배열 뒤에 마지막 전이와 방출
• Traceback = backpointer로 최적 경로 복원

Needleman-Wunsch가 서열-서열 격자였다면, Viterbi는 시각-상태 격자입니다. 곱셈이 덧셈으로, 정렬 점수가 로그 확률로 바뀌었을 뿐 뼈대는 완전히 같습니다.

이 유사성은 우연이 아닙니다. DP는 확률 곱셈, 최적화 덧셈, 최단 경로 등 여러 옷을 갈아입지만 뼈대는 하나입니다. 이걸 몸으로 익히면 나중에 나올 모든 확률 모델 알고리즘(Forward-Backward, Baum-Welch, Beam Search 등)이 같은 뿌리에서 뻗은 가지로 보입니다.

Rosalind에서 채점받기

Rosalind HMMV 문제를 열어봅시다. HMM 파라미터와 관측이 주어지고 Viterbi 상태 배열을 요구합니다. 위 코드를 그대로 붙이면 통과합니다. 자기 손으로 짠 30줄이 실제 채점을 통과하는 경험이 도구 100번 쓴 것보다 이해를 깊게 합니다.

다음 편으로 이어지는 갈래

• 다음 편 (M19): Forward-Backward 알고리즘 — Viterbi가 "가장 확률 높은 하나의 배열"이라면, Forward-Backward는 "각 시각 각 상태의 확률"을 뽑습니다. 두 알고리즘의 관계가 이 시리즈의 핵심입니다.
• 두 편 뒤 (M20): Baum-Welch EM — 지금은 확률표를 임의로 세팅했지만, 실제로는 관측 데이터만으로 이 표를 학습해야 합니다. EM 알고리즘의 자가 학습 루프.
• 세 편 뒤 (M21): Profile HMM과 gene finding — 단순한 2-상태 모델을 확장해서 유전자 구조(엑손, 인트론, 인터제닉)를 검출하는 실전 편.

더 깊게 파고 싶다면

본문은 BPD가 자체 재구성한 서술입니다. 원리를 이미 이해했다면 아래 정통 강의로 심화해도 좋습니다.

• Rabiner (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE — HMM의 정수를 담은 튜토리얼. 지금도 최고 참고 자료.
• MIT OCW 7.91J — Christopher Burge 교수의 HMM 강의 (자막 완비). 격자 유도를 판서로 봅니다.
• Durbin, Eddy, Krogh, Mitchison Biological Sequence Analysis — HMM을 바이오 서열에 적용하는 표준 교재.
• Rosalind Bioinformatics Textbook Track — HMMV, HMMFB, BW 세 문제를 순서대로 풀면 이 시리즈가 몸에 붙습니다.

다음 편에서 만나기 전에 위 예제를 반드시 손으로 한 번 채워봅시다. 곱셈 여덟 번을 종이에 적고, 결과가 저 위의 표와 맞는지 확인. 이 경험이 HMM의 나머지를 몸에 붙게 합니다.