气候统计方法的新成员——HMM（隐马尔可夫模型）

hillside

        HMM（隐马尔可夫模型）在语音识别、图像识别等领域已经取得了巨大的成功。但在气象气候领域的应用还非常欠缺，就中文文献而言，目前近于空白。
        去年猝然离世的南信大教授、气候学家丁裕国先生已经认识到HMM对于气候研究的意义，并已指导学生开展这方面的研究。
        MATLAB自带隐马尔可夫模型，统计工具箱包含5个与隐马尔可夫模型相关的函数：
                hmmgenerate 从一个马尔可夫模型产生一个状态序列和输出序列；
                hmmestimate 计算转移和输出的极大似然估计；
                hmmtrain 从一个输出序列计算转移和输出概率的极大似然估计；
                hmmviterbi 计算一个隐马尔可夫模型最可能的状态变化过程；
                hmmdecode 计算一个给定输出序列的后验状态概率。
        MATLAB自带的HMM据称有局限性，比如只能计算离散序列、不能设置初始分布等,一个具体的实例性用法见附件。网上有一些单独开发的HMM工具箱，现已上传国外高校与研究机构的多种HMM模型与使用信息。
        本人也是一个初学者，只能做入门性的介绍。
        附件《 东北地区土壤水分预报模型研究(采用HMM法)》（由丁裕国教授协助指导的2012年硕士论文）可称国内将HMM模型运用于气象气候领域乃至地学领域难得的率先性研究。
        至于作为隐马尔可夫过程基础的马尔可夫过程在气候统计研究中的运用，可阅读“气候统计”栏目另一帖《介绍一本原黄委总工、水文专家马秀峰关于游程理论的统计学新著  》中由丁裕国教授协助指导的硕士论文《Markov链模式在逐日降水持续性特征研究中的应用》（2012）。

附：[转载]Matlab2011b的HMM(隐马尔可夫模型)相关函数介绍

                                   

Matlab 2011b  Statistics Toolbox HMM

作者：yuheng666

Email:wuyuheng666@163.com

关键字：HMM，隐马尔科夫模型，Matlab，Statistics Toolbox

声明：本文主要介绍Matlab2011b中 Statistics Toolbox工具箱里与隐马尔科夫模型相关的函数及其用法（请勿与其它HMM工具箱混淆）。本文的主要内容来自Matlab 2011b的帮助文档，为作者自学笔记。水平有限，笔记粗糙，本着“交流探讨，知识分享”的宗旨，希望对HMM感兴趣的同学有些许帮助，欢迎指教，共同进步。

有关马尔科夫模型的基本知识，请参考其他资料。如：

http://www.eecis.udel.edu/~lliao/cis841s06/hmmtutorialpart1.pdf

http://www.eecis.udel.edu/~lliao/cis841s06/hmmtutorialpart2.pdf

http://cs229.stanford.edu/section/cs229-hmm.pdf

http://jedlik.phy.bme.hu/~gerjanos/HMM/node2.html

http://uirvli.ai.uiuc.edu/dugad/hmm_tut.html

.......

变量说明：

设有M个状态，N个符号Markov模型。

TRANS：对应状态转移矩阵，大小为M*M，表示各状态相互转换的概率，TRANS(i,j)表示从状态i转换到状态j的概率。

EMIS：对应符号生成矩阵，又叫混淆矩阵，观察符号概率分布。EMIS(i,j)代表在状态i时，产生符号j的概率。

函数介绍：

hmmgenerate — Generates a sequence of states and emissions from a Markov model

从一个马尔科夫模型产生状态序列和输出序列，该序列具有模型所表达的随机性特征。

A random sequence seq of emission symbols

A random sequence states of states

用法：

[seq,states] = hmmgenerate(len,TRANS,EMIS)

hmmgenerate(...,'Symbols',SYMBOLS)
hmmgenerate(...,'Statenames',STATENAMES)

示例：

trans = [0.95,0.05;    0.10,0.90];

emis = [1/6 1/6 1/6 1/6 1/6 1/6;

    1/10 1/10 1/10 1/10 1/10 1/2];

[seq,states] = hmmgenerate(100,trans,emis)

[seq,states] = hmmgenerate(100,trans,emis,...

    'Symbols',{'one','two','three','four','five','six'},...

'Statenames',{'fair';'loaded'})

估计状态序列：Estimating the State Sequence

hmmviterbi — Calculates the most probable state path for a hidden Markov model

Given the transition and emission matrices TRANS and EMIS, the function hmmviterbi uses the Viterbi algorithm to compute the most likely sequence of states the model would go through to generate a given sequence seq of emissions:

给定状态转移矩阵TRANS和混淆矩阵EMIS，hmmviterbi使用Viterbi算法计算该模型最相似的状态序列。

用法：

STATES = hmmviterbi(seq,TRANS,EMIS)
hmmviterbi(...,'Symbols',SYMBOLS)
hmmviterbi(...,'Statenames',STATENAMES)

示例：

likelystates = hmmviterbi(seq, TRANS, EMIS);

trans = [0.95,0.05;

         0.10,0.90];

emis = [1/6 1/6 1/6 1/6 1/6 1/6;

   1/10 1/10 1/10 1/10 1/10 1/2];

[seq,states] = hmmgenerate(100,trans,emis);

estimatedStates = hmmviterbi(seq,trans,emis);

[seq,states] = hmmgenerate(100,trans,emis,...

               'Statenames',{'fair';'loaded'});

estimatesStates = hmmviterbi(seq,trans,emis,...

              'Statenames',{'fair';'loaded'});

hmmestimate和hmmtrain用于通过给定的输出序列估计转移矩阵TRANS和混淆矩阵EMIS。

hmmestimate — Calculates maximum likelihood estimates of transition and emission probabilities from a sequence of emissions and a known sequence of states

通过一个输出序列和已知的状态序列，计算转移概率和输出概率的最大似然估计。

用法：

[TRANS,EMIS] = hmmestimate(seq,states)
hmmestimate(...,'Symbols',SYMBOLS)
hmmestimate(...,'Statenames',STATENAMES)
hmmestimate(...,'Pseudoemissions',PSEUDOE)
hmmestimate(...,'Pseudotransitions',PSEUDOTR)

示例：

通过已知的输出序列和状态序列估计出转移状态和混淆矩阵。

[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

hmmtrain — Calculates maximum likelihood estimates of transition and emission probabilities from a sequence of emissions.在知道输出序列，不知道状态转移序列，但是对转移矩阵和混淆矩阵有个初始猜测的情况下，可以使用hmmtrain估计转移状态和混淆矩阵，改函数可以选择使用BaumWelch或者Viterbi算法，通过迭代的方式进行估计，可以设置迭代次数Maxiterations和精度Tolerance。

用法：

[ESTTR,ESTEMIT] = hmmtrain(seq,TRGUESS,EMITGUESS)
hmmtrain(...,'Algorithm',algorithm)
hmmtrain(...,'Symbols',SYMBOLS)
hmmtrain(...,'Tolerance',tol)
hmmtrain(...,'Maxiterations',maxiter)
hmmtrain(...,'Verbose',true)
hmmtrain(...,'Pseudoemissions',PSEUDOE)
hmmtrain(...,'Pseudotransitions',PSEUDOTR)

示例：

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

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

You estimate TRANS and EMIS as follows:

[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

hmmdecode — Calculates the posterior state probabilities of a sequence of emissions

用于估计后验状态概率。对于一个输出序列seq的后验状态概率，是指这个模型在确定的状态下产生seq中一个符号的条件概率。

PSTATES = hmmdecode(seq,TRANS,EMIS)

The output PSTATES is an M-by-L matrix, where M is the number of states and L is the length of seq. PSTATES(i,j) is the conditional probability that the model is in state i when it generates the jth symbol of seq, given that seq is emitted.

PSTATES是一个M×L的矩阵，M是状态的个数，L是输出序列seq的长度。PSTATES(i,j)是在状态i时，产生seq的第j个符号的条件概率。

细心的读者可能会发现，上面都没有使用到Markov模型的初始状态概率矩阵。实际上，如果指定初始状态的概率，则上述函数默认从第一个状态开始，即初始状态为第一个状态的概率是1。但是statistics toolbox也提供了修改初始状态概率矩阵的方法。

设实际的状态转移矩阵和混淆矩阵分别是TRANS和EMIS，初始状态的概率分布为p=[p1,p2,p3,...,pm]，则可以通过以下方式重新设置初始状态。

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

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

kongfeng0824

matlab的工具箱很强大，这样写代码也很简洁。感谢分享！期待楼主更多更好的帖子。

kongfeng0824

matlab的工具箱很强大，这样写代码也很简洁。感谢分享！期待楼主更多更好的帖子。

yuanliang110

学习一下新方法

ayzqs

好强大，可以试着研究！！

小Q的梦

好复杂，不过是个好东西

meredithxf

好东西啊，值得学习的新方法

黎明的夕阳

学习一下新方法

ggfiri

好东西值得学学，谢谢分享

张学文

感想介绍了，理解还需要个过程。