理解隐马尔科夫（HMM）模型

前言

在李航的《统计学方法》第十章有对隐马尔科夫模型（Hidden Markov Model，HMM）比较详细的介绍和推导公式，我参考公式结合中文分词应用实现了隐马模型观测序列的生成、前向算法、维特比算法。

本文在此针对HMM模型在中文分词中的应用，讲讲实现原理。我尽可能的撇开公式，撇开推导。结合实际开源代码作为例子，争取做到雅俗共赏，童叟无欺。

没有公式，就没有伤害。

理解一个算法，我认为需要做到：会其意，知其形。本文回答的，其实主要是第一点。但是这一点呢，恰恰是最重要，而且很多书上不会讲的。

正如你在追一个姑娘，姑娘对你说“你什么都没做错！”你要是只看姑娘的表达形式呢，认为自己什么都没做错，显然就理解错了。你要理会姑娘的意思，“你赶紧给我道歉！”这样当你看到对应的表达形式呢，赶紧认错，跪地求饶就对了。数学也是一样，你要是不理解意思，光看公式，往往一头雾水。不过呢，数学的表达顶多也就是晦涩了点，姑娘的表达呢，有的时候就完全和本意相反。所以俺一直认为理解姑娘比理解数学难多了。

介绍

定义

wiki上有定义：

隐马尔可夫模型（Hidden Markov Model，HMM）是统计模型，它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。然后利用这些参数来作进一步的分析，例如模式识别。

马可夫模型的概率

这里用x表示状态， y表示观察值

假设观察到的结果为 Y:

Y=y(0),y(1),…,y(L-1)

隐藏条件为X:

X=x(0),x(1),…,x(L-1)

长度为 L，则马可夫模型的概率可以表达为：

详解

我们抛开教材上拗口的红球白球与盒子模型吧，来看一个简单的掷骰子的例子。

示例说明

假设我手里有三个不同的骰子。第一个骰子是我们平常见的骰子（称这个骰子为D6），6个面，每个面（1，2，3，4，5，6）出现的概率是1/6。第二个骰子是个四面体（称这个骰子为D4），每个面（1，2，3，4）出现的概率是1/4。第三个骰子有八个面（称这个骰子为D8），每个面（1，2，3，4，5，6，7，8）出现的概率是1/8。

假设我们开始掷骰子，我们先从三个骰子里挑一个，挑到每一个骰子的概率都是1/3。然后我们掷骰子，得到一个数字，1，2，3，4，5，6，7，8中的一个。不停的重复上述过程，我们会得到一串数字，每个数字都是1，2，3，4，5，6，7，8中的一个。

例如我们可能得到这么一串数字（掷骰子10次）：1 6 3 5 2 7 3 5 2 4

这串数字叫做可见状态链（对应上面公式的观察到的结果为 Y）。但是在隐马尔可夫模型中，我们不仅仅有这么一串可见状态链，还有一串隐含状态链（对应上面公式的隐藏的为 X）。

在这个例子里，这串隐含状态链就是你用的骰子的序列。比如，隐含状态链有可能是：D6 D8 D8 D6 D4 D8 D6 D6 D4 D8

一般来说，

HMM中说到的马尔可夫链其实是指隐含状态链，隐含状态（骰子）之间存在转换概率（transition probability）。

在我们这个例子里，D6的下一个状态是D4，D6，D8的概率都是1/3。D4的下一个状态是D4，D6，D8的转换概率也都一样是1/3。这样设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，我们可以这样定义，D6后面不能接D4，D6后面是D6的概率是0.9，是D8的概率是0.1，这样就是一个新的HMM，一般情况权重设定也确实是不一样的。

同样的，

尽管可见状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率（emission probability）。

就我们的例子来说，六面骰（D6）产生1的输出概率是1/6。产生2，3，4，5，6的概率也都是1/6。我们同样可以对输出概率进行其他定义。比如，我有一个被赌场动过手脚的六面骰子，掷出来是1的概率更大，是1/2，掷出来是2，3，4，5，6的概率是1/10。

而三个骰子之间也是可以相互转换的，其转换关系示意图如下所示。

其实对于HMM来说，如果提前知道所有隐含状态之间的转换概率和所有隐含状态到所有可见状态之间的输出概率，做模拟是相当容易的。

但是应用HMM模型时候呢，往往是缺失了一部分信息的，有时候你知道骰子有几种，每种骰子是什么，但是不知道掷出来的骰子序列；有时候你只是看到了很多次掷骰子的结果，剩下的什么都不知道。

如果应用算法去估计这些缺失的信息，就成了一个很重要的问题。这应该如何做呢？下面就来说明。

描述问题

这里就要顺带着说明与HMM模型相关的算法了，算法分为三类，分别对应着解决三种问题：

评估问题。

知道骰子有几种（隐含状态数量），每种骰子是什么（转换概率），根据掷骰子掷出的结果（可见状态链），我想知道掷出这个结果的概率。

这个问题看似意义不大，因为你掷出来的结果很多时候都对应了一个比较大的概率，问这个问题的目的呢，其实是检测观察到的结果和已知的模型是否吻合。如果很多次结果都对应了比较小的概率，那么就说明我们已知的模型很有可能是错的，有人偷偷把我们的骰子给换了。通常利用前向算法，分别计算每个产生给定观测序列的概率，然后从中选出最优的HMM模型。
解码问题。
知道骰子有几种（隐含状态数量），每种骰子是什么（转换概率），根据掷骰子掷出的结果（可见状态链），我想知道每次掷出来的都是哪种骰子（隐含状态链）。
学习问题。
知道骰子有几种（隐含状态数量），不知道每种骰子是什么（转换概率），观测到很多次掷骰子的结果（可见状态链），我想反推出每种骰子是什么（转换概率）。

这个问题很重要，因为这是最常见的情况。很多时候我们只有可见结果，不知道HMM模型里的参数，我们需要从可见结果估计出这些参数，这是建模的一个必要步骤，通常使用Baum-Welch算法解决。

解决方案

1.计算结果概率

知道骰子有几种，每种骰子是什么，每次掷的都是什么骰子，根据掷骰子掷出的结果，求产生这个结果的概率。

解法无非就是概率相乘：

2.计算隐含状态概率

计算不可见的隐含状态概率，破解骰子序列，这里有两种解法。

第一种解法，解最大似然路径问题

举个栗子：

我知道我有三个骰子，六面骰，四面骰，八面骰。我也知道我掷了十次的结果（1 6 3 5 2 7 3 5 2 4），我不知道每次用了那种骰子，我想知道最有可能的骰子序列。

其实最简单而暴力的方法就是穷举所有可能的骰子序列，然后依照第零个问题的解法把每个序列对应的概率算出来。然后我们从里面把对应最大概率的序列挑出来就行了。如果马尔可夫链不长，当然可行。如果长的话，穷举的数量太大，就很难完成了。

第二种解法，维特比算法（Viterbi algorithm）

维特比（Viterbi）算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划(dynamic programming)求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。

请不太理解动态规划算法的同学查看我之前的动态规划算法博文，现在我们来看看如何利用Vertibi算法计算骰子出现的概率。

还是举个栗子：

首先，如果我们只掷一次骰子：

看到结果为1。对应的最大概率骰子序列就是D4，因为D4产生1的概率是1/4，高于1/6和1/8。

把这个情况拓展，我们掷两次骰子：

结果为1，6。这时问题变得复杂起来，我们要计算三个值，分别是第二个骰子是D6，D4，D8的最大概率。显然，要取到最大概率，第一个骰子必须为D4。这时，第二个骰子取到D6的最大概率是：

同样的，我们可以计算第二个骰子是D4或D8时的最大概率。我们发现，第二个骰子取到D6的概率最大。而使这个概率最大时，第一个骰子为D4。所以最大概率骰子序列就是D4 D6。

继续拓展，我们掷三次骰子：

同样，我们计算第三个骰子分别是D6，D4，D8的最大概率。我们再次发现，要取到最大概率，第二个骰子必须为D6。这时，第三个骰子取到D4的最大概率是

同上，我们可以计算第三个骰子是D6或D8时的最大概率。我们发现，第三个骰子取到D4的概率最大。而使这个概率最大时，第二个骰子为D6，第一个骰子为D4。所以最大概率骰子序列就是D4 D6 D4。

小结

写到这里，大家应该看出点规律了。既然掷骰子一二三次可以算，掷多少次都可以以此类推。我们发现，我们要求最大概率骰子序列时要做这么几件事情。

首先，不管序列多长，要从序列长度为1算起，算序列长度为1时取到每个骰子的最大概率。

然后，逐渐增加长度，每增加一次长度，重新算一遍在这个长度下最后一个位置取到每个骰子的最大概率。因为上一个长度下的取到每个骰子的最大概率都算过了，重新计算的话其实不难。当我们算到最后一位时，就知道最后一位是哪个骰子的概率最大了。

最后，我们要把对应这个最大概率的序列从后往前推出来。

其他示例

还有一个来自wiki的经典的老王治病的HMM例子，请同学们自行跳转链接。

应用

中文分词中的应用

下面结合中文分词来说明HMM，HMM的典型介绍就是这个模型是一个五元组：

StatusSet: 状态值集合（隐状态）
ObservedSet: 观察值集合（输出文字集合）
TransProbMatrix: 转移概率矩阵（隐状态）
EmitProbMatrix: 发射概率矩阵（隐状态表现为显状态的概率）
InitStatus: 初始状态概率（隐状态）

结合参数说明HMM解决的三种问题：

参数(StatusSet, TransProbMatrix, EmitRobMatrix, InitStatus)已知的情况下，求解观察值序列。(Forward-backward算法)
参数(ObservedSet, TransProbMatrix, EmitRobMatrix, InitStatus)已知的情况下，求解状态值序列。(Viterbi算法)
参数(ObservedSet)已知的情况下，求解(TransProbMatrix, EmitRobMatrix, InitStatus)。(Baum-Welch算法)

其中，第二种问题是Viterbi算法求解状态值序列最常用，语音识别、中文分词、新词发现、词性标注都有它的一席之地。

五元组参数在中文分词中的具体含义，直接给五元组参数赋予具体含义：

StatusSet & ObservedSet

状态值集合为(B, M, E, S): {B:begin, M:middle, E:end, S:single}

分别代表每个状态代表的是该字在词语中的位置，B代表该字是词语中的起始字，M代表是词语中的中间字，E代表是词语中的结束字，S则代表是单字成词。

观察值集合为就是所有汉字字符串(“小明硕士毕业于中国科学院计算所”)，甚至包括标点符号所组成的集合。

状态值也就是我们要求的值，在HMM模型中文分词中，我们的输入是一个句子(也就是观察值序列)，输出是这个句子中每个字的状态值。

举个栗子：

小明硕士毕业于中国科学院计算所

输出的状态序列为：
BEBEBMEBEBMEBES

根据这个状态序列我们可以进行切词:
BE/BE/BME/BE/BME/BE/S

所以切词结果如下:
小明/硕士/毕业于/中国/科学院/计算/所

同时我们可以注意到：

B后面只可能接(M or E)，不可能接(B or S)。而M后面也只可能接(M or E)，不可能接(B, S)。

没错，就是这么简单，现在输入输出都明确了，下文讲讲输入和输出之间的具体过程，里面究竟发生了什么不可告人的秘密？

上面只介绍了五元组中的两元【StatusSet, ObservedSet】，下面介绍剩下的三元【InitStatus, TransProbMatrix, EmitProbMatrix】。

这五元的关系是通过一个叫Viterbi的算法串接起来， ObservedSet序列值是Viterbi的输入，而StatusSet序列值是Viterbi的输出，输入和输出之间Viterbi算法还需要借助三个模型参数，分别是InitStatus, TransProbMatrix, EmitProbMatrix，接下来一一讲解：

InitStatus：初始状态概率分布是最好理解的，可以示例如下：

#B
-0.26268660809250016
#E
-3.14e+100
#M
-3.14e+100
#S
-1.4652633398537678

PS：示例数值是对概率值取对数之后的结果(可以让概率相乘的计算变成对数相加)，其中-3.14e+100作为负无穷，也就是对应的概率值是0。下同。

也就是句子的第一个字属于{B,E,M,S}这四种状态的概率，如上可以看出，E和M的概率都是0，这和实际相符合，开头的第一个字只可能是词语的首字(B)，或者是单字成词(S)。

TransProbMatrix：转移概率是马尔科夫链很重要的一个知识点，大学里面学过概率论的人都知道，马尔科夫链最大的特点就是当前T=i时刻的状态Status(i)，只和T=i时刻之前的n个状态有关。也就是:

{Status(i-1), Status(i-2), Status(i-3), ... Status(i - n)}

更进一步的说，HMM模型有三个基本假设作为模型的前提，其中有个【有限历史性假设】，也就是马尔科夫链的n=1。即Status(i)只和Status(i-1)相关，这个假设能大大简化问题。

回过头看TransProbMatrix，其实就是一个4x4(4就是状态值集合的大小)的二维矩阵，示例如下：

矩阵的横坐标和纵坐标顺序是BEMS x BEMS。(数值是概率求对数后的值，别忘了。)

-3.14e+100 -0.510825623765990 -0.916290731874155 -3.14e+100
-0.5897149736854513 -3.14e+100 -3.14e+100 -0.8085250474669937
-3.14e+100 -0.33344856811948514 -1.2603623820268226 -3.14e+100
-0.7211965654669841 -3.14e+100 -3.14e+100 -0.6658631448798212

比如TransProbMatrix[0][0]代表的含义就是从状态B转移到状态B的概率，由TransProbMatrix[0][0] = -3.14e+100可知，这个转移概率是0，这符合常理。由状态各自的含义可知，状态B的下一个状态只可能是ME，不可能是BS，所以不可能的转移对应的概率都是0，也就是对数值负无穷，在此记为-3.14e+100。

由上TransProbMatrix矩阵可知，对于各个状态可能转移的下一状态，且转移概率对应如下：

#B
#E:-0.510825623765990,M:-0.916290731874155
#E
#B:-0.5897149736854513,S:-0.8085250474669937
#M
#E:-0.33344856811948514,M:-1.2603623820268226
#S
#B:-0.7211965654669841,S:-0.6658631448798212

EmitProbMatrix：这里的发射概率(EmitProb)其实也是一个条件概率而已，根据HMM模型三个基本假设里的“观察值独立性假设”，观察值只取决于当前状态值，也就是:

P(Observed[i], Status[j]) = P(Status[j]) * P(Observed[i]|Status[j])

其中P(Observed[i]|Status[j])这个值就是从EmitProbMatrix中获取。

EmitProbMatrix示例如下：

#B
耀:-10.460283,涉:-8.766406,谈:-8.039065,伊:-7.682602,洞:-8.668696,...
#E
耀:-9.266706,涉:-9.096474,谈:-8.435707,伊:-10.223786,洞:-8.366213,...
#M
耀:-8.47651,涉:-10.560093,谈:-8.345223,伊:-8.021847,洞:-9.547990,....
#S
蘄:-10.005820,涉:-10.523076,唎:-15.269250,禑:-17.215160,洞:-8.369527...

虽然EmitProbMatrix也称为矩阵，这个矩阵太稀疏了，实际工程中一般是将上面四行发射转移概率存储为4个Map，详见我的github代码。

到此，已经介绍完HMM模型的五元参数，假设现在手头上已经有这些参数的具体概率值，并且已经加载进来，(也就是该模型的字典)，那么我们只剩下Viterbi这个算法函数，这个模型就算可以开始使用了。

贴一个jieba分词（java版）的Viterbi算法代码：

public void viterbi(String sentence, List<String> tokens) {
        Vector<Map<Character, Double>> v = new Vector<Map<Character, Double>>();
        Map<Character, Node> path = new HashMap<Character, Node>();

        v.add(new HashMap<Character, Double>());
        for (char state : states) {
            Double emP = emit.get(state).get(sentence.charAt(0));
            if (null == emP)
                emP = MIN_FLOAT;
            v.get(0).put(state, start.get(state) + emP);
            path.put(state, new Node(state, null));
        }

        for (int i = 1; i < sentence.length(); ++i) {
            Map<Character, Double> vv = new HashMap<Character, Double>();
            v.add(vv);
            Map<Character, Node> newPath = new HashMap<Character, Node>();
            for (char y : states) {
                Double emp = emit.get(y).get(sentence.charAt(i));
                if (emp == null)
                    emp = MIN_FLOAT;
                Pair<Character> candidate = null;
                for (char y0 : prevStatus.get(y)) {
                    Double tranp = trans.get(y0).get(y);
                    if (null == tranp)
                        tranp = MIN_FLOAT;
                    tranp += (emp + v.get(i - 1).get(y0));
                    if (null == candidate)
                        candidate = new Pair<Character>(y0, tranp);
                    else if (candidate.freq <= tranp) {
                        candidate.freq = tranp;
                        candidate.key = y0;
                    }
                }
                vv.put(y, candidate.freq);
                newPath.put(y, new Node(y, path.get(candidate.key)));
            }
            path = newPath;
        }
        double probE = v.get(sentence.length() - 1).get(‘E‘);
        double probS = v.get(sentence.length() - 1).get(‘S‘);
        Vector<Character> posList = new Vector<Character>(sentence.length());
        Node win;
        if (probE < probS)
            win = path.get(‘S‘);
        else
            win = path.get(‘E‘);

        while (win != null) {
            posList.add(win.value);
            win = win.parent;
        }
        Collections.reverse(posList);

        int begin = 0, next = 0;
        for (int i = 0; i < sentence.length(); ++i) {
            char pos = posList.get(i);
            if (pos == ‘B‘)
                begin = i;
            else if (pos == ‘E‘) {
                tokens.add(sentence.substring(begin, i + 1));
                next = i + 1;
            }
            else if (pos == ‘S‘) {
                tokens.add(sentence.substring(i, i + 1));
                next = i + 1;
            }
        }
        if (next < sentence.length())
            tokens.add(sentence.substring(next));
    }

对算法有疑问的可以参考这段动画，将代码单步一遍，什么都明白了：