Matlab统计学工具箱之（隐）马尔可夫模型：Markov Models

原文：https://chunqiu.blog.ustc.edu.cn/?p=454

目录

马尔科夫链

隐马尔科夫模型（HMM）

HMM简介

分析隐马尔科夫模型

产生测试序列

估计状态序列

估计状态转移矩阵和输出矩阵

估计后验状态概率

改变初始状态分布

 

马尔科夫链

         如下图的马尔科夫模型状态转移矩阵

 

         图中的方框代表你要建模过程的可能状态，箭头表示状态间的转移情况。箭头上的标签表示转移概率。在过程中的每一步，模型会产生一个输出（决定于所处的状态）然后转移到另一个状态。马尔科夫模型的一个重要特性是下一个状态只与当前的状态有关，与过去的转移无关。

         例如，对于一个顺序抛掷硬币的过程，存在两个状态：正和反。最近一次的投掷决定了模型当前的状态，下一次投掷决定了如何转移到下一状态。如果硬币是公平的，则转移概率都是二分之一，输出就是当前的状态。在更复杂的模型中，每一个状态可以以一个随机过程产生输出，如可以通过投掷骰子来决定输出。

         马尔科夫链是离散状态马尔科夫模型的数学描述。马尔科夫链如下定义：

l  一系列状态{1,2,…,M}；

l  一个MXM的状态转移矩阵T，其中Tij表示从状态i到j的状态转移概率。矩阵每行的和为1，因为它表示从一个给定状态到其他各状态的所有转移概率的和。

l  一系列可能的输出{s1,s2,…,sN}，默认时可以当做{1,2,…,N}，N为可能输出的个数，当然你可以选择不同的数据集或符号来表示；

l  一个MXN的输出矩阵E，其中Eik表示模型在状态i时的输出sk的概率。

马尔科夫链的过程：从第0步的初始状态i0开始；然后以概率T1i1转移到状态i1，并以概率Ei1k1输出sk1；最后，以概率T1i1 Ei1k1 Ti1i2 Ei2k2 … Tir-1ir Eirk观察到状态i1i2…ir，同时输出sk1sk2…skr。

 

隐马尔科夫模型（HMM）

HMM简介

         隐马尔科夫模型是你只观察到输出，但不知道模型产生输出所经历的状态。隐马尔科夫模型分析就是从观察的数据恢复相关的状态序列。

         例如，考虑一个有2个状态6个可能输出的马尔科夫模型，模型中用到

n  一个红骰子，6个面，分别标号1~6；

n  一个绿骰子，12个面，5个面标号2~6，其他面标记1；

n  一个非均衡的红硬币，以0.9的概率掷出正，0.1的概率掷出反；

n  一个非均衡的绿硬币，以0.95的概率掷出正，0.05的概率掷出反。

模型以下面的规则产生一系列来自{1,2,3,4,5,6}的数字序列：

n  先掷红骰子，记下骰子的输出；

n  然后掷红硬币：

l  如果是正，则掷红骰子记下结果；

l  否则掷绿骰子记结果；

n  接下来的下一步，掷与上一步所用骰子同颜色的硬币，如果是正，则掷同颜色的骰子；否则掷另一个骰子。

模型的状态转移图如下

 

         这个模型中，通过掷与状态同颜色的骰子来决定输出，通过掷与状态同颜色的硬币来决定转移到哪一状态。

         状态转移矩阵为

 

         输出矩阵为

 

         这个模型是非隐的，因为你从硬币、骰子的颜色可以知道状态序列。如果假设有人产生了一系列输出，但不给你展示硬币和骰子，所有你能看到的只是输出序列。如果你开始看到1出现的次数比其他都多，你可能会假设模型处于绿状态，但你不能确定，因为你不能看到所掷骰子的颜色。

         隐马尔科夫模型提出了如下问题：

1)         给定输出序列，求最可能的状态序列？

2)         给定输出序列，你如何估计模型的转移概率和输出概率？

3)         模型得出给定序列的先验概率？

4)         模型在序列某一位置是特定状态的后验概率？

分析隐马尔科夫模型

         与隐马尔科夫模型有关的统计工具箱函数如下

l  Hmmgenerate——从一个马尔科夫模型产生状态、输出序列；

l  Hmmestimate——从输出序列和已知的状态序列计算转移、输出序列概率的极大似然估计；

l  Hmmtrain——从输出序列计算转移、输出序列概率的极大似然估计；

l  Hmmviterbi——计算隐马尔科夫模型的最大可能状态序列；

l  Hmmdecode——计算输出序列的状态后验概率。

下面来讲如何使用这些函数分析隐马尔科夫模型。

产生测试序列

         下面的命令创建状态转移矩阵和输出矩阵

TRANS = [.9 .1; .05 .95;];

 

EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...

7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

        要产生一组状态和输出的随机序列，使用函数hmmgenerate

[seq,states] = hmmgenerate(1000,TRANS,EMIS);

输出参数seq为输出序列，states为状态序列。

         Hmmgenerate函数第0步开始于状态1，然后第一步转移到状态i1，而i1作为states的第一个元素。如果要改变初始状态，看后面的改变初始状态分布一节。

估计状态序列

         给定状态转移矩阵TRANS和输出矩阵EMIS，函数hmmviterbi使用Viterbi算法计算产生给定输出序列seq的极大可能状态序列

likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates与seq长度相同。

         为了测试hmmviterbi的精确度，计算实际状态序列states与估计状态序列likelystates吻合的百分比

sum(states==likelystates)/1000

ans =

   0.8200

         在这个例子中，估计的状态和实际吻合达82%

估计状态转移矩阵和输出矩阵

函数hmmestimate和hmtrain可以用来根据给定的输出序列估计状态转移矩阵和输出矩阵。

使用hmmestimate

         此函数需要你知道产生输出序列的状态序列。基于输出序列和状态序列，如下估计转移矩阵和输出矩阵

[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)

 

TRANS_EST =

0.8989    0.1011

0.0585    0.9415

 

EMIS_EST =

0.1721    0.1721    0.1749    0.1612    0.1803    0.1393

0.5836    0.0741    0.0804    0.0789    0.0726    0.1104

         与原始的序列对比

TRANS

TRANS =

0.9000    0.1000

0.0500    0.9500

 

EMIS

EMIS =

0.1667    0.1667    0.1667    0.1667    0.1667    0.1667

0.5833    0.0833    0.0833    0.0833    0.0833    0.0833

使用hmmtrain

         如果你不知道状态序列，但知道TRANS和EMIS矩阵的初始猜测，可以使用hmmtrain函数估计。

         假设已知TRANS和EMIS矩阵的初始猜测如下

TRANS_GUESS = [.85 .15; .1 .9];

EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

下面做矩阵估计

[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)

 

TRANS_EST2 =

0.2286    0.7714

0.0032    0.9968

 

EMIS_EST2 =

0.1436    0.2348    0.1837    0.1963    0.2350    0.0066

0.4355    0.1089    0.1144    0.1082    0.1109    0.1220

         Hmmtrain函数使用迭代算法调整TRANS_GUESS和EMIS_GUESS，使得基于调整的矩阵产生更加类似于观察序列seq的输出。当两次迭代中矩阵的变化不大时迭代结束。

         如果如果算法在默认的最大迭代次数100次以内仍没结束，将返回迭代的最后结果并产生一个警告。这时可以通过增加最大迭代次数来增大迭代的次数

hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)

maxiter是新设的最大迭代次数。

还可以改变默认的最小容忍值

hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

tol为新设的最小容忍值，增大它可以让算法更快的停止，但结果将不精确。

         两个因素将降低hmmtrain函数估计结果的可靠性：

l  算法收敛到一个局部极大值，但此极大值并不是真正的状态转移矩阵与输出矩阵的位置。对此，可以使用多个不同的初始猜测值来估计，使用最好的结果；

l  Seq序列太短，以至于不足以训练矩阵，对此，可以使用更长的seq序列。

估计后验状态概率

输出序列seq的后验状态概率是模型给出seq序列时所处状态的后验概率。可以使用hmmdecode函数计算后验状态概率：

PSTATES = hmmdecode(seq,TRANS,EMIS)

         随着序列长度的增加，序列的概率会趋向于0，所以一个足够长的序列的概率将小于计算机的最小精度，为了避免此问题，hmmdecode还返回了概率的对数值。

改变初始状态分布

         默认地，统计工具箱的隐马尔科夫模型函数是从状态1开始的。换句话说，初始状态的分布集中在状态1处，其他位置为0。如果想给初始状态赋予一个概率分布p=[p1,p2,…,pM]，可以如下实现

1. 创建一个M+1XM+1的辅助状态转移矩阵如下

 

其中的T是真正的状态转移矩阵，第一列包含M+1个0，向量p的和为1。

1. 创建M+1XN的辅助输出矩阵如下

 

         如果状态转移矩阵与输出矩阵式TRANS和EMIS，可以使用如下命令创建辅助矩阵

TRANS_HAT = [0 p; zeros(size(TRANS,1),1) TRANS];

 EMIS_HAT = [zeros(1,size(EMIS,2)); EMIS];

 

原文：http://www.aiseminar.cn/bbs/forum.php?mod=viewthread&tid=596

此文讲述的内容在Matlab 7.0、7.5（R2007b）中均有――马尔可夫工具箱，主要内容如下。

简介：马尔可夫处理是随机处理的一个典型例子――此种处理根据特定的概率产生随机输出或状态序列。马尔可夫处理的特别之处在于它的无记忆性――他的下一个状态仅依赖他的当前状态，不考虑导致他们的历史。马尔可夫处理的模型在实际应用中使用非常广泛，从每日股票价格到染色体中的基因位置都有应用。

马尔可夫链

马尔可夫模型用状态图可视化描述如下。
<ignore_js_op> 

在图中，矩形代表你要描述的模型在处理中可能出现的状态，箭头描述了状态之间的转换。每个箭头上的标签表明了该转换出现的概率。在处理的每一步，模型都可能根据当前状态产生一种output或emission，然后做一个到另一状态的转换。马尔可夫模型的一个重要特点是：他的下个状态仅仅依赖当前状态，而与导致其成为当前状态的历史变换无关。

马尔可夫链是马尔可夫模型的一组离散状态集合的数学描述形式。马尔可夫链特征归纳如下：

1. 一个状态的集合{1, 2, ..., M}
2. 一个M * M的转移矩阵T，(i, j)位置的数据是从状态i转到状态j的概率。T的每一行值的和必然是1，因为这是从一个给定状态转移到其他所有可能状态的概率之和。
3. 一个可能的输出（output）或发布（emissions）的集合{S1, S2, ..., SN}。默认情况下，发布的集合是{1, 2, ..., N}，这里N是可能的发布的个数，当然，你也可以选择一个不同的数字或符号的集合。
4. 一个M * N的发布矩阵E，(i, k)入口给出了从状态i得到发布的标志Sk的概率。

马尔可夫链在第0步，从一个初始状态i0开始。接着，此链按照T(1, i1)概率转移到状态i1，且按概率E(i1, k1)概率发布一个输出S(k1)。因此，在前r步，状态序列(i1, i2, i3, ..., ir)和发布序列(Sk1, Sk2, ..., Skr)的可能的观测结果是
T(1, i1)E(i, k1), T(i1, i2)E(i2, k2), ..., T(ir-1, ir)E(ir, k)


隐马尔可夫模型

简介：隐马尔可夫模型相对马尔可夫模型的不同之处在于，你观测到一组发布序列，但是却不知道模型通过什么样的状态序列得到这样的发布。隐马尔可夫模型分析就是要从观测数据恢复出这一状态序列。

一个例子：考虑一个拥有2个状态和6个位置发布的马尔可夫模型。模型描述如下：

1. 一个红色骰子，有6个面，标记为1到6。
2. 一个绿色骰子，有12个面，有5个标记为2到6，余下的全标记为1。
3. 一个加权的红硬币，正面的概率是0.9，背面的概率是0.1。
4. 一个加权的绿硬币，正面的概率为0.95，背面的概率为0.05。

这个模型按照下面的规则从集合{1, 2, 3, 4, 5, 6}产生一个数字序列：

1. 从滚动红骰子开始，记下出现的数字，作为发布结果。
2. 投红硬币：
如果结果是正面，滚动红骰子，记下结果；
如果结果是背面，滚动绿骰子，记下结果。
3. 在下面的每一步，你都抛和前一步所投骰子相同颜色的硬币。如果硬币出现正面，滚和硬币相同颜色的骰子。如果硬币出现反面，改为投另种颜色的骰子。

这个模型的状态图如下，有两个状态，红和绿：
<ignore_js_op> 
通过投掷和状态一样颜色的骰子来解决输出（发布），通过抛同样颜色的硬币来决定状态的转移。

状态转移概率矩阵为：
>> T = [0.9 0.1; 0.05 0.95]
T =
    0.9000    0.1000
    0.0500    0.9500

输出或发布概率矩阵为：
>> E = [1/6 1/6 1/6 1/6 1/6 1/6; 7/12 1/12 1/12 1/12 1/12 1/12]
E =
    0.1667    0.1667    0.1667    0.1667    0.1667    0.1667
    0.5833    0.0833    0.0833    0.0833    0.0833    0.0833

这个模型并不是隐藏的，因为我们从硬币和骰子的颜色已经知道状态序列。假设，有其他人产生了一个发布结果，而没有向你展示硬币和骰子，你能看到的只有结果。当你看到1比其他数字多时，你也许猜测这个模型是在绿色状态，但是因为你不能看到被投骰子的颜色，所以你并不能确定。

隐马尔可夫模型提出了如下问题：

1. 给定一个输出序列，最有可能的状态过程是什么？
2. 给定一个输出序列，如何估计模型的转移和发布概率？
3. 如何求这个模型产生一个给定序列的先验概率？
4. 如何求这个模型产生一个给定序列的后验概率？


隐马尔可夫模型分析

统计工具箱包含5个与隐马尔可夫模型相关的函数：
hmmgenerate 从一个马尔可夫模型产生一个状态序列和输出序列；
hmmestimate 计算转移和输出的极大似然估计；
hmmtrain 从一个输出序列计算转移和输出概率的极大似然估计；
hmmviterbi 计算一个隐马尔可夫模型最可能的状态变化过程；
hmmdecode 计算一个给定输出序列的后验状态概率。

下面部分介绍如何使用这些函数来分析隐马尔可夫模型。

1. 产生一个测试序列

下面代码产生上面简介中模型的转移和输出矩阵：
TRANS = [.9 .1; .05 .95;];
EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...
7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

要从模型产生一个随机的状态序列和输出序列，使用hmmgenerate：
[seq,states] = hmmgenerate(1000,TRANS,EMIS);

输出中，seq是输出序列，states是状态序列。hmmgenerate在第0步从状态1开始，在第一步转移到状态i1
，并返回i1作为状态的第一个入口。

2. 估计状态序列

给定了转移和输出矩阵TRANS和EMIS，函数hmmviterbi使用Viterbi算法计算模型给定输出序列seq最有可能
通过的状态序列：
likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates是和seq一样长的序列。计算hmmvertibi的精度如下：
sum(states == likelystates) / length(states)
ans =
    0.8680

3. 估计转移和输出矩阵
函数hmmestimate和hmmtrain用于估计给定输出序列seq的转移和输出矩阵TRANS和EMIS。

使用hmmestimate
[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)
TRANS_EST =
    0.9065    0.0935
    0.0406    0.9594
EMIS_EST =
    0.1452    0.1516    0.1581    0.1968    0.1581    0.1903
    0.5841    0.0754    0.0986    0.0812    0.0841    0.0768
由上面使用方式可知，hmmestimate函数需要事先知道了得到输出序列seq，以及得到此结果的状态变化序
列。

使用hmmtrain
如果不知道状态序列，但是知道TRANS和EMIS的初始猜测，那就可以使用hmmtrain来估计TRANS和EMIS。

假设已知如下初始猜测：
TRANS_GUESS = [.85 .15; .1 .9];
EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

TRANS和EMIS的估计如下：
[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)
TRANS_EST2 =
    0.9207    0.0793
    0.0370    0.9630
EMIS_EST2 =
    0.1792    0.1437    0.1436    0.1855    0.1509    0.1971
    0.5774    0.0775    0.1042    0.0840    0.0859    0.0710

hmmtrain使用迭代算法来不断修改TRANS_GUESS和EMIS_GUESS，使得每一步修改得到的矩阵都更加可能产生观测序列seq。当前后两个两次迭代矩阵的变化在一个小的容错范围内时，迭代停止。如果算法无法达到容错的范围，则迭代到达一定次数时就会停止，并返回一个警告提示。默认的最大迭代次数为100。

如果算法达不到目标误差范围，则可以通过增加迭代次数和/或加大容错误差值来使其获得较合适结果：
改变迭代次数maxiter：hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)
改变容错误差tol：hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

影响hmmtrain输出的矩阵可靠性的两点因素：
（1）算法收敛于局部极值，这点可以使用不同的初始猜测矩阵来尝试解决；
（2）序列seq太短而无法很好的训练矩阵，可以尝试使用较长的序列。

4. 估计后验状态概率（不太理解）
一个输出序列seq的后验状态概率是在特定状态下的模型产生在seq中一个输出的条件概率。假定seq已经给出，你可以使用hmmdecode得到后验状态概率。

PSTATES = hmmdecode(seq,TRANS,EMIS)

输出为一个M * N的矩阵。M是状态的个数，L是seq的长度。PSTATES(i, j)是模型在状态i时，产生seq第j个输出的条件概率。

参考文献：
[1] Matlab R2007b帮助文档，http://www.aiseminar.cn/bbs，jink2005简译！

 

 

----------------------------------------------------------------------------------------------------------------

统计工具箱包含5个与隐马尔可夫模型相关的函数：
hmmgenerate 从一个马尔可夫模型产生一个状态序列和输出序列；
hmmestimate 计算转移和输出的极大似然估计；
hmmtrain 从一个输出序列计算转移和输出概率的极大似然估计；
hmmviterbi 计算一个隐马尔可夫模型最可能的状态变化过程；
hmmdecode 计算一个给定输出序列的后验状态概率。

下面部分介绍如何使用这些函数来分析隐马尔可夫模型。

1. 产生一个测试序列

下面代码产生上面简介中模型的转移和输出矩阵：
TRANS = [.9 .1; .05 .95;];
EMIS = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6;...
7/12, 1/12, 1/12, 1/12, 1/12, 1/12];

要从模型产生一个随机的状态序列和输出序列，使用hmmgenerate：
[seq,states] = hmmgenerate(1000,TRANS,EMIS);

输出中，seq是输出序列，states是状态序列。hmmgenerate在第0步从状态1开始，在第一步转移到状态i1
，并返回i1作为状态的第一个入口。

2. 估计状态序列

给定了转移和输出矩阵TRANS和EMIS，函数hmmviterbi使用Viterbi算法计算模型给定输出序列seq最有可能
通过的状态序列：
likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates是和seq一样长的序列。计算hmmvertibi的精度如下：
sum(states == likelystates) / length(states)
ans =
    0.8680

3. 估计转移和输出矩阵
函数hmmestimate和hmmtrain用于估计给定输出序列seq的转移和输出矩阵TRANS和EMIS。

使用hmmestimate
[TRANS_EST, EMIS_EST] = hmmestimate(seq, states)
TRANS_EST =
    0.9065    0.0935
    0.0406    0.9594
EMIS_EST =
    0.1452    0.1516    0.1581    0.1968    0.1581    0.1903
    0.5841    0.0754    0.0986    0.0812    0.0841    0.0768
由上面使用方式可知，hmmestimate函数需要事先知道了得到输出序列seq，以及得到此结果的状态变化序
列。

使用hmmtrain
如果不知道状态序列，但是知道TRANS和EMIS的初始猜测，那就可以使用hmmtrain来估计TRANS和EMIS。

假设已知如下初始猜测：
TRANS_GUESS = [.85 .15; .1 .9];
EMIS_GUESS = [.17 .16 .17 .16 .17 .17;.6 .08 .08 .08 .08 08];

TRANS和EMIS的估计如下：
[TRANS_EST2, EMIS_EST2] = hmmtrain(seq, TRANS_GUESS, EMIS_GUESS)
TRANS_EST2 =
    0.9207    0.0793
    0.0370    0.9630
EMIS_EST2 =
    0.1792    0.1437    0.1436    0.1855    0.1509    0.1971
    0.5774    0.0775    0.1042    0.0840    0.0859    0.0710

hmmtrain使用迭代算法来不断修改TRANS_GUESS和EMIS_GUESS，使得每一步修改得到的矩阵都更加可能产生观测序列seq。当前后两个两次迭代矩阵的变化在一个小的容错范围内时，迭代停止。如果算法无法达到容错的范围，则迭代到达一定次数时就会停止，并返回一个警告提示。默认的最大迭代次数为100。

如果算法达不到目标误差范围，则可以通过增加迭代次数和/或加大容错误差值来使其获得较合适结果：
改变迭代次数maxiter：hmmtrain(seq,TRANS_GUESS,EMIS_GUESS,'maxiterations',maxiter)
改变容错误差tol：hmmtrain(seq, TRANS_GUESS, EMIS_GUESS, 'tolerance', tol)

影响hmmtrain输出的矩阵可靠性的两点因素：
（1）算法收敛于局部极值，这点可以使用不同的初始猜测矩阵来尝试解决；
（2）序列seq太短而无法很好的训练矩阵，可以尝试使用较长的序列。

4. 估计后验状态概率（不太理解）
一个输出序列seq的后验状态概率是在特定状态下的模型产生在seq中一个输出的条件概率。假定seq已经给出，你可以使用hmmdecode得到后验状态概率。

PSTATES = hmmdecode(seq,TRANS,EMIS)

输出为一个M * N的矩阵。M是状态的个数，L是seq的长度。PSTATES(i, j)是模型在状态i时，产生seq第j个输出的条件概率。