HMM示例及Matlab计算

Alice 和Bob是好朋友，但是他们离得比较远，每天都是通过电话了解对方那天作了什么。Bob仅仅对三种活动感兴趣：公园散步,购物以及清理房间。他选择做什么事情只凭当天天气。Alice对于Bob所住的地方的天气情况并不了解，但是知道总的趋势。在Bob告诉Alice每天所做的事情基础上，Alice想要猜测Bob所在地的天气情况。

　　Alice认为天气的运行就像一个马尔可夫链。其有两个状态 “雨”和”晴”，但是无法直接观察它们，也就是说,它们对于Alice是隐藏的。每天，Bob有一定的概率进行下列活动：“散步”，“购物”，或 “清理”。因为Bob会告诉Alice他的活动，所以这些活动就是Alice的观察数据。这整个系统就是一个隐马尔可夫模型HMM。

　　Alice知道这个地区的总的天气趋势，并且平时知道Bob会做的事情。也就是说这个隐马尔可夫模型的参数是已知的。可以用程序语言(Python)写下来：

　　　// 状态数目，两个状态：雨或晴

　　　states = (‘Rainy’, ‘Sunny’)

　　　// 每个状态下可能的观察值

　　　observations = (‘walk’, ‘shop’, ‘clean’)

　　　//初始状态空间的概率分布

　　　start_probability = {‘Rainy’: 0.6, ‘Sunny’: 0.4}

　　　// 与时间无关的状态转移概率矩阵

　　　transition_probability = {

　　　’Rainy’ : {‘Rainy’: 0.7, ‘Sunny’: 0.3},

　　　’Sunny’ : {‘Rainy’: 0.4, ‘Sunny’: 0.6},

　　　}

　　　//给定状态下，观察值概率分布,发射概率

　　　emission_probability = {

　　　’Rainy’ : {‘walk’: 0.1, ‘shop’: 0.4, ‘clean’: 0.5},

　　　’Sunny’ : {‘walk’: 0.6, ‘shop’: 0.3, ‘clean’: 0.1},

　　　}

　　在这些代码中，start_probability代表了Alice对于Bob第一次给她打电话时的天气情况的不确定性(Alice知道的只是那个地方平均起来下雨多些)。在这里,这个特定的概率分布并非平衡的，平衡概率应该接近（在给定变迁概率的情况下）{‘Rainy’: 0.571, ‘Sunny’: 0.429}。 transition_probability 表示马尔可夫链下的天气变迁情况，在这个例子中,如果今天下雨，那么明天天晴的概率只有30%。代码emission_probability 表示了Bob每天作某件事的概率。如果下雨，有 50% 的概率他在清理房间；如果天晴，则有60%的概率他在外头散步。

　　Alice和Bob通了三天电话后发现第一天Bob去散步了，第二天他去购物了，第三天他清理房间了。Alice现在有两个问题：这个观察序列“散步、购物、清理”的总的概率是多少？(注：这个问题对应于HMM的基本问题之一：已知HMM模型λ及观察序列O，如何计算P(O|λ)？) 最能解释这个观察序列的状态序列（晴/雨）又是什么？（注：这个问题对应HMM基本问题之二：给定观察序列O=O1,O2,…OT以及模型λ，如何选择一个对应的状态序列S = q1,q2,…qT，使得S能够最为合理的解释观察序列O？）

　　至于HMM的基本问题之三：如何调整模型参数，使得P(O|λ)最大？这个问题事实上就是给出很多个观察序列值，来训练以上几个参数的问题。

 

For illustration, the section uses the example described in A Concrete Example for Hidden Markov Model.

 

First, create the transition and emission matrices by entering the following commands.

%%%%%%%%%%%%%%%%%%%%%

转移矩阵

TRANS=[0.7 0.3;0.4 0.6];

 

发射概率

EMIS=[0.1 0.4 0.5;0.6 0.3 0.1];

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%

Next, generate a random sequence of emissions from the model, seq, of length 1000, using the function hmmgenerate. You can also return the corresponding random sequence of states in the model as the second output, states.

 

Generating Data

%%%%%%%%%%%%%%%%%%%%%%%%%

随机生成状态和概率

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

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

得到这么一段东西：

seq=[1, 2, 2, 2, 3, 2,...]

说明朋友最近的行踪：【散步，购物，购物，购物，打扫卫生，购物】，是个购物狂哦！

 

Computing the Most Likely Sequence of States

 

Used the Viterbi algorithm to compute the most likely sequence of states that the model would go through to generate the given sequence of emissions.

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

likelystates = hmmviterbi(seq, TRANS, EMIS);

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

To test the accuracy of hmmviterbi, you can compute the percentage of the time that the actual sequence states agrees with the sequence likelystates.

 

%%%%%%%%%%%%%%%%%%%%%%%%%%

sum(states==likelystates)/1000

%%%%%%%%%%%%%%%%%%%%%%%%%

 

I got 0.7710 in my computer. 看来对天气的预测不是很准哈

 

Estimating the Transition and Emission Matrices

Suppose you do not know the transition and emission matrices in the model, and you observe a sequence of emissions, seq. If you have an initial guess as to the values of TRANS and EMIS, you can estimate the transition and emission matrices using the function hmmtrain.

For example, suppose you have the following initial guesses for TRANS and EMIS.

 

 

%%%%%%%%%%%%%%%%%%%%%%%%

TRANS_GUESS = [0.5 0.5;0.5 0.5];

EMIS_GUESS = [0.3 0.4 0.3;0.6 0.3 0.2];

%%%%%%%%%%%%%%%%%%%%%%%%

 

You can estimate TRANS and EMIS with the following command.

 

%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%