HMM 隐马尔可夫模型 MATLAB

HMM 隐马尔可夫模型

Intro

维基百科定义：隐马尔可夫模型用来描述一个含有隐含未知参数的马尔可夫过程，而马尔可夫过程是是一个具备了马尔可夫性质的随机过程，当一个随机过程在给定现在状态及所有过去状态情况下，其未来状态的条件概率分布仅依赖于当前状态；换句话说，在给定现在状态时，它与过去状态（即该过程的历史路径）是条件独立的，那么此随机过程即具有马尔可夫性质。

Matlab HMM 定义：

A hidden Markov model (HMM) is one in which you observe a sequence of emissions, but do not know the sequence of states the model went through to generate the emissions. Analyses of hidden Markov models seek to recover the sequence of states from the observed data.

人话（个人理解）：

我们可以做一个假设：世界上的很多随机现象，比如天气（晴天或是下雨天），看似随机其实背后都是有其真正的原因，我们假设这些随机表象背后是由数个真正的原因所决定的，而且每一个状态只与上一个状态相关（今天的天气只与昨天天气背后的隐藏状态有关），我们称满足这样假设的模型为马尔可夫模型。我们在这里称表象（晴天或是下雨天）为emissions（发射）或者the observed data，称背后真正的原因为hidden state（隐藏状态）。之所以将表象称为发射，是因为一个表象是由一个隐藏状态“发射”而来，每一天都是由一个隐藏状态产生，而每一个隐藏状态都可能产生其他的表象，产生表象的概率矩阵被称为发射矩阵\(B\)。每一天的隐藏状态之间也存在转移概率，此概率矩阵被称为转移矩阵\(A\)。初始隐藏状态的概率矩阵被称为\(\pi\)

于是这就是隐马尔可夫模型(HMM)的重要三个矩阵\(A,B,\pi\)

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

• 给定的发射序列，求最可能的状态序列
• 给定的发射序列，你将如何估计模型的转移概率和输出概率？
• 模型生成给定序列 的 先验概率是多少？
• 在给定序列中某个状态时，模型的后验概率是多少？

Hidden Markov models raise the following questions:

• Given a sequence of emissions, what is the most likely state path?
• Given a sequence of emissions, how can you estimate transition and emission probabilities of the model?
• What is the forward probability that the model generates a given sequence?
• What is the posterior probability that the model is in a particular state at any point in the sequence?

模型需要指定一个转移矩阵TRANS和一个发射矩阵EMIS，TRANS(i,j)是隐藏状态i转移到j的概率，EMIS(i,j) 是序列 seq 在隐藏状态i时发射出符号\(j\)的概率。

设\(N\)为隐藏状态的状态个数，\(M\)为发射序列(emissions)中元素(symbol)可能的值的个数，则TRANS矩阵大小为\(N\times N\)，EMIS矩阵大小为\(N\times M\)

时间复杂度，再设\(V\)为发射序列(emissions)中元素(symbol)的维数，则易知隐藏状态\(x(t)\)转移到\(x(t+1)\)复杂度为\(N\times N\)，隐藏状态\(x(t)\)发射一个元素\(y(t)\)的复杂度为\(N\times M \times V\)

实现

MATLAB中提供了相关的工具箱

生成测试数列

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

输出发射序列seq 和序列状态 states 。

测试状态序列

函数 hmmviterbi 使用算法Viterbi计算模型生成序列seq最有可能经过的状态序列

likelystates = hmmviterbi(seq, TRANS, EMIS);

likelystates 序列的长度和 seq相同。

测试hmmviterbi的准确度可以将序列likelystates和原状态序列states比较

sum(states==likelystates)/1000
ans =
   0.8200

此样例可知准确度为82%

估算转移矩阵TRANS和发射矩阵EMIS

函数 hmmestimate和hmmtrain 可以通过发射序列seq估算模型的转移矩阵TRANS 和发射矩阵 EMIS 。

使用hmmestimate函数时需要序列的states （即在生成测试序列时的返回值 states ）

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

如果不知道序列的 states并且对于转移矩阵TRANS 和发射矩阵 EMIS有初步猜想，可以可以使用hmmtrain函数来估算序列的转移矩阵TRANS 和发射矩阵 EMIS

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

hmmtrain 使用Baum-Welch迭代算法，通过改变TRANS_GUESS 和 EMIS_GUESS 以便让模型更有可能生成观察序列 seq。当两个连续迭代中的矩阵彼此公差很小的时候，算法停止。

如果算法在默认值为100的最大迭代次数内未能达到此公差，则算法将停止，hmmtrain将返回TRANS_EST和EMIS_EST的最后一个值，并发出未达到公差的警告。

如果算法未能达到所需的公差，请使用以下命令增加最大迭代次数的默认值：

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

其中，maxiter是算法执行的最大步数。
使用以下命令更改公差的默认值：

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

其中，tol是公差的期望值。增加tol的值会使算法更快地停止，但结果不太准确。

有两个因素会降低hmmtrain输出矩阵的可靠性：

• 该算法收敛到一个局部极大值，该极大值不代表真实的过渡矩阵和发射矩阵。如果您怀疑这一点，请对矩阵TRANS_EST和EMIS_EST使用不同的初始猜测。
• 序列seq可能太短，无法正确训练矩阵。如果您怀疑这一点，请对seq使用更长的序列。

估算先验后验概率

可以通过函数hmmdecode计算已经发射的序列seq在某一历史状态下的后验概率

PSTATES = hmmdecode(seq,TRANS,EMIS)

输出值PSTATES是一个M-by-L（M乘L）矩阵，这里M是状态的数量，L是序列seq的长度。PSTATES(i,j)是模型在状态i时生成序列seq中符号j的条件概率。

更详细的介绍：

PSTATES = hmmdecode(seq,TRANS,EMIS)
[PSTATES,logpseq] = hmmdecode(...)
[PSTATES,logpseq,FORWARD,BACKWARD,S] = hmmdecode(...)
hmmdecode(...,'Symbols',SYMBOLS)

PSTATES = hmmdecode(seq,TRANS,EMIS) 可以根据隐马尔可夫(HMM)模型，计算出序列seq后面状态的可能性PSTATES。 后续状态的可能性是k和i的条件概率（给定观察序列的符号sym）。你应该通过指定一个转移矩阵TRANS和一个发射矩阵EMIS来使用模型。TRANS(i,j)是状态i转移到j的概率，EMIS(k,seq) 是序列 seq 由状态k发射出的概率。

（输出值PSTATES的介绍见上）

PSTATES = hmmdecode(seq,TRANS,EMIS) calculates the posterior state probabilities, PSTATES, of the sequence seq, from a hidden Markov model. The posterior state probabilities are the conditional probabilities of being at state k at step i, given the observed sequence of symbols, sym. You specify the model by a transition probability matrix, TRANS, and an emissions probability matrix, EMIS. TRANS(i,j) is the probability of transition from state i to state j. EMIS(k,seq) is the probability that symbol seq is emitted from state k.

注意：函数hmmdecode在第一次发射之前，状态1从step 0开始。hmmdecode根据模型从状态 1 开始计算PSTATES中的概率。

**Note **The function hmmdecode begins with the model in state 1 at step 0, prior to the first emission. hmmdecode computes the probabilities in PSTATES based on the fact that the model begins in state 1.

[PSTATES,logpseq] = hmmdecode(...) 返回 给定转移矩阵TRANS和发射矩阵EMIS下的对数概率 logpseq

[PSTATES,logpseq] = hmmdecode(...) returns logpseq, the logarithm of the probability of sequence seq, given transition matrix TRANS and emission matrix EMIS.

[PSTATES,logpseq,FORWARD,BACKWARD,S] = hmmdecode(...) 返回由S缩放后的先验概率和后验概率

[PSTATES,logpseq,FORWARD,BACKWARD,S] = hmmdecode(...) returns the forward and backward probabilities of the sequence scaled by S.

hmmdecode(...,'Symbols',SYMBOLS) 指定了发射的符号。SYMBOLS可以是一个数值数组、字符串数组或者是一个包含符号的cell array 。默认的符号是整数1到N，N是可能的发射数量。

hmmdecode(...,'Symbols',SYMBOLS) specifies the symbols that are emitted. SYMBOLS can be a numeric array, a string array, or a cell array of the names of the symbols. The default symbols are integers 1 through N, where N is the number of possible emissions.

后记

评价：不如Python hmmlearn，MATLAB再见

附

其中，关于矩阵的指定（参见）：

基本上来说，对于这个问题我们没有一种简单直接的答案。相反的，经验告诉我们无论是使用随机（subject to the stochastic and the nonzero value constraints）还是通过先验概率和转移矩阵的初始估计，都足以在几乎所有情况下对参数进行重新调整。

但是……（未翻译，参见原文）

可以通过多种方式获得矩阵参数

1. 手动将观测序列分割为状态，并在状态内平均观测值
2. 用平均值对观测值进行最大似然分割
3. 带聚类的分段k-均值分段
4. etc

Basically, there is no simple or straightforward answer to the above question. Instead, experience has shown that either random (subject to the stochastic and the nonzero value constraints) or uniform initial estimates of the prior probabilities and the transition matrix is adequate for giving useful reestimates of these parameters in almost all cases.

However, for the emission matrix, experience has shown that good initial estimates are helpful in the discrete symbol case, and are essential (when dealing with multiple mixtures) in the continuous distribution case.

Such initial estimates can be obtained in a number of ways, including:

1) manual segmentation of the observation sequences into states with averaging observations within states,

2) maximum likelihood segmentation of observations with averaging,

3) segmental k-means segmentation with clustering,

etc.