viterbi算法实现

文章由 http://huangjian.info/blog/232/viterbi/ 整理获得。

Viterbi以它发明者的名字命名。

Viterbi算法：给出一个观测序列o1,o2,o3 …，我们希望找到观测序列背后的隐藏状态序列s1, s2, s3, …；这样一种由动态规划的方法来寻找出现概率最大的隐藏状态序列（被称为Viterbi路径）的方法。

隐藏状态序列被成为：隐马可夫序列（HMM，Hidden Markov Model）。

Markov随机过程具有如下的性质：在任意时刻，从当前状态转移到下一个状态的概率与当前状态之前的那些状态没有关系。Viterbi要用到这个性质。

HMM的一个任务是通过观测序列来找到背后的隐藏序列。另外两个任务是：a) 给定一个HMM，计算一个观测序列出现的可能性；b)已知一个观测序列，HMM参数不定，如何优化这些参数使得观测序列的出现概率最大。解决前一个问题可以用与Viberbi结构非常类似的Forward算法来解决（实际上在下面合二为一），而后者可以用Baum-Welch/EM算法来迭代逼近。

Wiki上一个例子。

假设你有一个朋友在外地，每天你可以通过电话来了解他每天的活动。他每天只会做三种活动之一——Walk, Shop, Clean。你的朋友从事哪一种活动的概率与当地的气候有关，这里，我们只考虑两种天气——Rainy, Sunny。我们知道，天气与运动的关系如下：

	Rainy	Sunny
Walk	0.1	0.6
Shop	0.4	0.3
Clean	0.5	0.1

例如，在下雨天出去散步的可能性是0.1。

而天气之前互相转换的关系如下，（从行到列）

	Rainy	Sunny
Rainy	0.7	0.3
Sunny	0.4	0.6

例如，从今天是晴天而明天就开始下雨的可能性是0.4 。

同时为了求解问题我们假设初始情况：通话开始的第一天的天气有0.6的概率是Rainy，有0.4概率是Sunny。OK，现在的问题是，如果连续三天，你发现你的朋友的活动是：Walk->Shop->Clean；那么，如何判断你朋友那里这几天的天气是怎样的？

python代码：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
 X = ['Rainy', 'Sunny'] = [0, 1]
 y = ['Walk', 'Shop', 'Clean'] = [0, 1, 2]

  天气状态转移:
  Rainy -> Rainy = 0.7, Rainy -> Sunny = 0.3
  Sunny -> Rainy = 0.4, Sunny -> Sunny = 0.6
  tp = [[0.7, 0.3], [0.4, 0.6]]

  某天气下某行动的概率：
  Walk|Rainy = 0.1, Shop|Rainy = 0.4, Clean|Rainy = 0.5
  Walk|Sunny = 0.6, Shop|Sunny = 0.3, Clean|Sunny = 0.1
  ep = [[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]]

  第一天天气的概率：
  Rainy = 0.6, Sunny = 0.4
  sp = [0.6, 0.4]
"""

def forward_viterbi(y, X, sp, tp, ep):
    T = {}
    for state in X:
        ##  (prob,  V.path, V.prob )
        T[state] = (sp[state], [state], sp[state])
    for output in y:
        U = {}
        for next_state in X:
            total = 0
            argmax = None
            valmax = 0
            for source_state in X:
                (prob, v_path, v_prob) = T[source_state]
                p = ep[source_state][output] * tp[source_state][next_state]
                prob *= p
                v_prob *= p
                total += prob
                if v_prob > valmax:
                    argmax = v_path + [next_state]
                    valmax = v_prob
            U[next_state] = (total, argmax, valmax)
        T = U

    ## apply sum/max to the final states:
    total = 0
    argmax = None
    valmax = 0
    for state in X:
        (prob, v_path, v_prob) = T[state]
        total += prob
        if v_prob > valmax:
            argmax = v_path
            valmax = v_prob
    return (total, argmax, valmax)

if __name__ == '__main__':
    X = [0, 1]
    y = [0, 1, 2]
    sp = [0.6, 0.4]
    tp = [[0.7, 0.3], [0.4, 0.6]]
    ep = [[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]]
    ret = forward_viterbi(y, X, sp, tp, ep)
    print(ret)

output：(0.033611999999999996, [1, 0, 0, 0], 0.009407999999999998)

所以朋友那边这几天最可能的天气情况是Sunny->Rainy->Rainy->Rainy，它有0.009408的概率出现。而我们算法的另一个附带的结论是，我们所观察到的朋友这几天的活动序列：Walk->Shop->Clean在我们的隐马可夫模型之下出现的总概率是0.033612（也是Forward算法的输出）。

几点说明：

1. 算法对于每一个状态要记录一个三元组：(prob, v_path, v_prob)，其中，prob是从开始状态到当前状态所有路径（不仅仅 是最有可能的viterbi路径）的概率加在一起的结果（作为算法附产品，它可以输出一个观察序列在给定HMM下总的出现概率，即forward算法的输出），v_path是从开始状态一直到当前状态的viterbi路径，v_prob则是该路径的概率。
2. 算法开始，初始化T （T是一个Map，将每一种可能状态映射到上面所说的三元组上）
3. 三重循环，对每个一活动y，考虑下一步每一个可能的状态next_state，并重新计算若从T中的当前状态state跃迁到next_state概率会有怎样的变化。跃迁主要考虑天气转移（tp[source_state][next_state]）与该天气下从事某种活动（ep[source_state][output]）的联合概率。所有下一步状态考虑完后，要从T中找出最优的选择viterbi路径——即概率最大的viterbi路径，即上面更新Map U的代码U[next_state] = (total, argmax, valmax)。
4. 算法最后还要对T中的各种情况总结，对total求和，选择其中一条作为最优的viterbi路径。
5. 算法输出四个天气状态，这是因为，计算第三天的概率时，要考虑天气转变到下一天的情况。

还有一个没有定义变量的C++代码，和python很像，仅仅作为参考：

#include <iostream>
using namespace std;

void forward_viterbi(const vector<string> & ob, viterbi_triple_t & vtriple)
{
    //alias

    map<string, double>& sp = start_prob;
    map<string, map<string, double> > & tp = transition_prob;
    map<string, map<string, double> > & ep = emission_prob;

    // initialization
    InitParameters();

    map<string, viterbi_triple_t> T;

    for (vector<string>::iterator it = states.begin(); it != states.end(); ++it)
    {
        viterbi_triple_t foo;
        foo.prob = sp[*it];
        foo.vpath.push_back(*it);
        foo.vprob = sp[*it];

        T[*it] = foo;
    }

    map<string, viterbi_triple_t> U;
    double total = 0;
    vector<string> argmax;
    double valmax = 0;
    double p = 0;

    for (vector<string>::const_iterator itob = ob.begin(); itob != ob.end(); ++itob)
    {
        cout << “observation=” << *itob << endl;
        U.clear();
        for (vector<string>::iterator itNextState = states.begin(); itNextState != states.end(); ++itNextState)
        {
            cout << “\tnext_state=” << *itNextState << endl;
            total = 0;
            argmax.clear();
            valmax = 0;
            for (vector<string>::iterator itSrcState = states.begin(); itSrcState != states.end(); ++itSrcState)
            {
                cout << “\t\tstate=” << *itSrcState << endl;
                viterbi_triple_t foo = T[*itSrcState];
                p = ep[*itSrcState][*itob] * tp[*itSrcState][*itNextState];
                cout << “\t\t\tp=” << p << endl;
                foo.prob *= p;
                foo.vprob *= p;
                cout << “\t\t\ttriple=” << foo << endl;
                total += foo.prob;
                if (foo.vprob > valmax)
                {
                    foo.vpath.push_back(*itNextState);
                    argmax = foo.vpath;
                    valmax = foo.vprob;
                }
            }

            U[*itNextState] = viterbi_triple_t(total, argmax, valmax);
            cout << “\tUpdate U[" << *itNextState << "]=” << U[*itNextState] << “” << endl;
        }
        T.swap(U);
    }

    total = 0;
    argmax.clear();
    valmax = 0;
    for (vector<string>::iterator itState = states.begin(); itState != states.end(); ++itState)
    {
        viterbi_triple_t foo = T[*itState];
        total += foo.prob;
        if (foo.vprob > valmax)
        {
            argmax.swap(foo.vpath);
            valmax = foo.vprob;
        }
    }

    vtriple.prob = total;
    vtriple.vpath = argmax;
    vtriple.vprob = valmax;
    cout << “final triple=” << vtriple << endl;
}