隐马尔科夫模型 Viterbi算法 Python实现

已知初始状态概率向量 pi，状态转移概率矩阵 P，发射概率矩阵 B，求观察序列为 O 的条件下状态序列为 Q 的概率。代码如下：

def HMM_ProCond(pi, P, B, Q, O):
    Q, O = [q - 1 for q in Q], [o - 1 for o in O]
    N, proNume, proDeno, proCond = len(pi[0]), - numpy.log(pi[0][Q[0]]) - numpy.log(B[Q[0]][O[0]]), 0, []
    for i in range(N):
        proCond.append(- numpy.log(pi[0][i]) - numpy.log(B[i][O[0]]))
    for i in range(1, len(Q)):
        proNume -= numpy.log(P[Q[i - 1]][Q[i]]) + numpy.log(B[Q[i]][O[i]])
        proCondi = []
        for j in range(N):
            proCondi.append(0)
            for k in range(N):
                proCondi[-1] += numpy.exp(numpy.log(P[k][j]) - proCond[k])
            proCondi[-1] = - numpy.log(proCondi[-1]) - numpy.log(B[j][O[i]])
        proCond = proCondi
    for i in range(N):
        proDeno += numpy.exp(- proCond[i])
    return numpy.exp(- proNume - numpy.log(proDeno))

已知 pi, P, B 如上所述，求观察序列为 O 的条件下概率最大的状态序列 Q。该问题用 Viterbi 算法解决，代码如下：

def HMM_Viterbi(pi, P, B, O):
    N, Q, V = len(pi[0]), [], []
    for i in range(N):
        Q.append([i + 1])
        V.append(- numpy.log(pi[0][i]) - numpy.log(B[i][O[0] - 1]))
    for Oi in O[1 :]:
        Qi, Vi = [], []
        for j in range(N):
            Qi.append([])
            Vi.append(-1)
            for k in range(N):
                Vk = V[k] - numpy.log(P[k][j]) - numpy.log(B[j][Oi - 1])
                if Vi[-1] == -1 or Vi[-1] > Vk:
                    Qi[-1], Vi[-1] = Q[k] + [j + 1], Vk
        Q, V = Qi, Vi
    m = 0
    for i in range(N):
        if V[m] > V[i]:
            m = i
    return Q[m]

Viterbi 算法本质上是动态规划，设 V[i, j] = Max{P(Q[1 … i] | O[1 … i] and Q[i] = j)}，则有状态转移方程 V[i, j] = B[j, O[i]] * Max{V[i - 1, k] * P[k, j]}，边界 V[1, j] = B[j, O[1]] * pi[j]，于是 Max{V[len(O), j]} 即为所求 Max{P(Q, O)}，由条件概率公式得 P(Q|O) 最大时 P(Q, O) 亦为最大，每次迭代时记录状态即可得到 Q。

给个例子，pi = $\left[\begin{array}{cc}0.6& 0.4\end{array}\right]$, P = $\left[\begin{array}{cc}0.6& 0.4\\ 0.3& 0.7\end{array}\right]$, B = $\left[\begin{array}{ccc}0.6& 0.3& 0.1\\ 0.1& 0.4& 0.5\end{array}\right]$，当 O = [1, 2, 3], Q = [1, 2, 2] 时，求 P(Q|O)；当 O = [2, 1, 3, 2, 3, 2, 2, 3, 3, 1, 2, 1, 1, 1, 2, 3, 3, 3, 3, 2] 时，求 Q 使 P(Q|O) 最大。代码如下：

import numpy

pi, P, B = [[0.6, 0.4]], [[0.6, 0.4], [0.3, 0.7]], [[0.6, 0.3, 0.1], [0.1, 0.4, 0.5]]
N = len(pi[0])

print(HMM_ProCond(pi, P, B, [1, 2, 2], [1, 2, 3]))
print(HMM_Viterbi(pi, P, B, [2, 1, 3, 2, 3, 2, 2, 3, 3, 1, 2, 1, 1, 1, 2, 3, 3, 3, 3, 2]))