# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np


def optimal_path_viterbi(pi,A,B,O):
    """
    :param pi: 初始概率
    :param A: 状态转移概率
    :param B: 观测概率
    :param O: 观测序列
    :return: 最优路径的状态索引、概率最大值矩阵、节点矩阵
    """
    N = A.shape[0] # 状态个数
    T = len(O) #观测序列长度
    # 在时刻t 状态为i 的所有单个路径 (i1,i2,...,it) 中概率最大值 为 delta：δ
    delta = np.zeros([T,N],dtype=np.float64)
    # 在时刻t 状态为i 的所有单个路径 (i1,i2,...,it) 中概率最大 的路径的 前一个节点 t-1 为 psi：ψ
    psi = np.zeros([T,N],dtype=np.int) #
    # 初始化
    for i in range(N):
        delta[0][i] = pi[i]*B[i][0]
    # 递推公式
    for t in range(1,T):
        for i in range(N):
            tmp = [delta[t-1][j]*A[j][i] for j in range(N)]
            psi[t][i] = np.argmax(tmp)
            delta[t][i] = tmp[psi[t][i]] * B[i][O[t]]
    # 最优路径回溯
    i_star=np.zeros([T],dtype=np.int)
    # 初始化最后一个最优节点
    i_star[T-1] = np.argmax(delta[T-1])
    # 从后向前递推
    for t in range(T-2,-1,-1):
        i_star[t] = psi[t+1][i_star[t+1]]
    return i_star, delta, psi


if __name__ == '__main__':

    pi = np.array([0.2,0.4,0.4])

    A = np.array([[0.5, 0.2,0.3],
                 [0.3,0.5,0.2],
                 [0.2,0.3,0.5]])

    B = np.array([[0.5,0.5],
                 [0.4,0.6],
                 [0.7,0.3]])
    O = np.array([0,1,0]) # (红，白，红)

    i_star, delta, psi = optimal_path_viterbi(pi,A,B,O)
    
    print '最优路径索引 i*[t]：\n',i_star
    print '概率最大值矩阵 δ[t][i]：\n',delta
    print '节点矩阵 ψ[t][i]：\n', psi

最终输出：

最优路径索引 i*[t]：
[2 2 2]
概率最大值矩阵 δ[t][i]：
[[0.1     0.16    0.28   ]
 [0.028   0.0504  0.042  ]
 [0.00756 0.01008 0.0147 ]]
节点矩阵 ψ[t][i]：
[[0 0 0]
 [2 2 2]
 [1 1 2]]