Day 046-Q5 — トレリスDP + Viterbi 最尤経路（HMM・log-space）

2026-05-30 赤色 Master / Phase 8+ ★★★★★★★★★ Viterbi / HMM / Forward / log-sum-exp

問題

隠れマルコフモデル（HMM）が与えられる：

  • 状態数 $K$（状態 $1, \ldots, K$）、観測記号数 $V$
  • 初期確率 $\pi_k$（状態 $k$ から始まる確率）
  • 遷移確率 $A_{ij}$（状態 $i$ から $j$ へ遷移する確率）
  • 発行確率 $B_{kv}$（状態 $k$ が記号 $v$ を発行する確率）
  • 観測列 $O_1, O_2, \ldots, O_T$

以下の2つを求めよ：

  1. 最尤状態系列（Viterbi）: $\arg\max_{q_1,\ldots,q_T} P(O, q | \lambda)$
  2. 観測列の対数確率（Forward）: $\log P(O | \lambda)$（小数点以下3桁）

制約

$1 \le K \le 50$
$1 \le V \le 50$
$1 \le T \le 10^5$
全確率は正（$> 0$）
時間制限: 2秒

入出力例

入力例 1

2 2 4
0.6 0.4
0.7 0.3 0.4 0.6
0.5 0.5 0.4 0.6
1 2 1 2

出力例 1

1 1 1 1
-5.411

概念図: HMM トレリス（K=2, T=4）

t=1 t=2 t=3 t=4 q₁=1 q₂=1 q₃=1 q₄=1 q₁=2 q₂=2 q₃=2 q₄=2 O₁=1 O₂=2 O₃=1 O₄=2 Viterbi

ヒント（段階的開示）

ヒント1: 方向性
確率の積が非常に小さくなるため、log-space で計算します。積 → 和、max → max（そのまま）、sum → log_sum_exp に変換します。
ヒント2: Viterbi アルゴリズム
  • $\delta_t(k) = \max_{q_1,\ldots,q_{t-1}} \log P(O_{1:t}, q_t=k | \lambda)$
  • $\delta_0(k) = \log \pi_k + \log B_{k,O_0}$
  • $\delta_t(k) = \max_j [\delta_{t-1}(j) + \log A_{jk}] + \log B_{k,O_t}$
  • $\psi_t(k) = \arg\max_j$ を保存してバックトレース
ヒント3: log_sum_exp の実装 
import math
def log_sum_exp(values):
    m = max(values)
    if m == -math.inf:
        return -math.inf
    return m + math.log(sum(math.exp(v - m) for v in values
                            if v > -math.inf))

# Forward アルゴリズム: Viterbi の max → log_sum_exp に置換

模範解答 (Python)

import sys
import math
input = sys.stdin.readline
NEG_INF = -math.inf

def log_sum_exp(vals):
    m = max(vals)
    if m == NEG_INF:
        return NEG_INF
    return m + math.log(sum(math.exp(v - m) for v in vals if v > NEG_INF))

def main():
    K, V, T = map(int, input().split())
    pi = list(map(float, input().split()))
    A_flat = list(map(float, input().split()))
    A = [[A_flat[i*K+j] for j in range(K)] for i in range(K)]
    B_flat = list(map(float, input().split()))
    B = [[B_flat[i*V+j] for j in range(V)] for i in range(K)]
    obs = [int(x)-1 for x in input().split()]

    log_pi = [math.log(x) if x > 0 else NEG_INF for x in pi]
    log_A  = [[math.log(A[i][j]) if A[i][j] > 0 else NEG_INF for j in range(K)] for i in range(K)]
    log_B  = [[math.log(B[k][v]) if B[k][v] > 0 else NEG_INF for v in range(V)] for k in range(K)]

    # === Viterbi ===
    delta = [[NEG_INF]*K for _ in range(T)]
    psi   = [[-1]*K for _ in range(T)]
    for k in range(K):
        delta[0][k] = log_pi[k] + log_B[k][obs[0]]

    for t in range(1, T):
        for k in range(K):
            best_val = NEG_INF
            best_j = -1
            lbk = log_B[k][obs[t]]
            for j in range(K):
                v = delta[t-1][j] + log_A[j][k]
                if v > best_val:
                    best_val = v
                    best_j = j
            delta[t][k] = best_val + lbk
            psi[t][k] = best_j

    q_star = [0] * T
    q_star[T-1] = max(range(K), key=lambda k: delta[T-1][k])
    for t in range(T-2, -1, -1):
        q_star[t] = psi[t+1][q_star[t+1]]

    viterbi_seq = ' '.join(str(q+1) for q in q_star)

    # === Forward (log-sum-exp) ===
    alpha = [NEG_INF] * K
    for k in range(K):
        alpha[k] = log_pi[k] + log_B[k][obs[0]]

    for t in range(1, T):
        new_alpha = [NEG_INF] * K
        for k in range(K):
            vals = [alpha[j] + log_A[j][k] for j in range(K)]
            new_alpha[k] = log_sum_exp(vals) + log_B[k][obs[t]]
        alpha = new_alpha

    log_prob = log_sum_exp(alpha)

    print(viterbi_seq)
    print(f"{log_prob:.3f}")

main()

Step-by-Step 解説

1HMM の基本構造
隠れ状態 $q_t$ が Markov 連鎖を形成し、各時刻に観測 $O_t$ を発行。確率の積が $T \to 10^5$ でアンダーフローするため log-space が必須。
2Viterbi アルゴリズム $O(TK^2)$
$\delta_t(k)$ を dp 配列として更新。Viterbi は「max + 加算」の DP。バックトレース配列 $\psi$ を保存して最後から逆追いする。
3Forward アルゴリズム $O(TK^2)$
Viterbi の max を log_sum_exp に置き換えたもの。$\log P(O|\lambda) = \text{lse}(\alpha_{T-1})$。
4log_sum_exp の数値安定性
最大値 $m$ を引いてから exp → sum → log → $m$ を戻す。オーバーフロー・アンダーフロー防止の標準テクニック。
5計算量の確認
内側ループ $K^2$、外側 $T$ 回。$T=10^5, K=50$ で $2.5 \times 10^8$ 演算。PyPy 推奨、Python では numpy 活用が現実的。

計算量

Viterbi: $O(TK^2)$
Forward: $O(TK^2)$
バックトレース: $O(TK)$
空間: $O(TK)$（delta/psi 配列）

よくあるミス

ミス原因正しい書き方
確率のまま積を取る$T$ が大きいとアンダーフローlog-space で和を取る
log_sum_exp に $-\infty$ が混入math.log(0) エラーif v > NEG_INF でフィルタ
Viterbi と Forward を混同max vs sumViterbi=max、Forward=lse
バックトレース方向を逆にする$T-1$ から $0$ へ逆追いq_star[t] = psi[t+1][q_star[t+1]]

次のステップ

  • 発展問題: Baum-Welch アルゴリズム（EM法によるHMMパラメータ推定）
  • 関連: Day039 Q5 吸収マルコフ連鎖（確率行列のガウス消去）
  • 応用: 音声認識・自然言語処理（品詞タグ付け）

自己評価