隐马尔科夫模型(Hidden Markov Model,HMM)是机器学习中一个极其重要又有效的方法，说到隐马尔科夫模型，就要先说一下马尔科夫过程，马尔科夫过程指的是一个东西，它有N种状态，而在时刻T的状态只与前面的k个状态有关，k≥1，称之为k阶马尔科夫模型，我们一般讨论的都是k=1的情况。这里要说明的是，假设k=1，下一个状态只与上一个状态有关，但是并不是上一个状态确定了，下一个状态就确定了，而是下一个状态有多种情况，而具体出现哪一种情况是一个概率事件，这个概率决定于上一个状态是什么。比如说天气问题，假设第二天的天气只跟前一天的天气有关，但是今天下雨，明天的天气可以是90%下雨，10%晴天，而今天晴天，明天可能20%下雨，80%晴天这样。

我们可以知道，一个一阶马尔科夫模型主要由初始状态概率分布π，状态转移矩阵A和状态数目N组成。初始概率分布就是t=0的时候各个状态出现的概率，状态转移矩阵就是每个状态下之间转换的概率矩阵，状态数目就是有几个状态。

那什么是HMM呢？一个很经典的例子，我们有很多个不同的盒子，每个盒子了里面都有很多不同颜色的小球，这些盒子我们称之为状态，每次我们只能选定一个盒子，t=0的时候每个盒子都有一定概率被选中，假设随机选了某一个盒子，然后我们就从盒子里面随机挑出一个小球，然后我们根据盒子的状态概率转移矩阵挑选出下一个盒子，然后在下一个盒子里面再随机挑出一个小球，然后不断重复，隐马尔科夫模型研究的就是这么一个过程。只是啊，在这个过程中我们只能看到每次被挑出来的小球的颜色，而不能知道小球每次具体是从哪个盒子里面挑出来的，盒子（状态）之间的转换过程是隐藏的！！

所以隐马尔科夫模型组成的参数就变成了：初始状态概率分布π(就是t=0的时候每个盒子被选中的概率)，盒子(状态)之间的概率转移矩阵A，每个盒子挑出每种小球的概率B，观测到的小球的状态序列O（就是输出结果），状态数目N，小球的种类M。一般符号λ=(A,B,π)；

可以看出，HMM是一个双随机过程，其中一个过程不可见。

随便找本讲HMM的教程或课本都会提到HMM主要就是解决三种问题：

1. 评估：给定一个隐马尔科夫模型λ，和输出结果O，计算出这个模型下出现这种结果的概率。
2. 解码：给定一个隐马尔科夫模型λ，和输出结果O，请问最有可能出现这种输出结果的状态转移序列是怎么样的？
3. 学习：给定输出结果O，M,N,计算出λ来。

首先是出现观测序列的概率的计算问题，即P(O|λ),可以用计算机暴力计算所有路径下的概率，然后相加，也就是所谓的穷举法，但是不会有人推荐这么做的，计算机都表示压力很大，尤其是状态数目比较多的情况下小心计算机罢工！

所以呢，就有两种动态规划的方法可以相对高速地计算出这个结果来，一般称之为向前算法(Forward Algorithm)和向后算法(Backward Algorithm)；

所谓向前算法，就是我们引进一个变量α_t(i)，它的计算方式如下：

α_t(i)=P(o₁o₂o₃…o_t,q_t= i|λ)

上式中o_t表示t时刻的观测结果，q_t表示时刻t的状态，很明显α_t(i)的含义就是给定HMM模型λ的情况下，在1到t时刻观测结果分别是o₁o₂o₃…o_t,t时刻的状态是i的概率。显然在t=1的时候有下式成立：

α₁(i)=π_ib_i(o₁) (1 ≤i≤N )

上式中b_i(o₁)状态i出现观测结果o₁的概率，其实就是矩阵B中第i行第o₁个数。

假设有了α_t(i)(1 ≤i≤N )，就可以计算出α_t+1(i)了：

\(\alpha_{t+1}(j)=[\sum_{t=1}^N\alpha_t(j)a_{ij}]b_j(o_{t+1}),1\le t\le T-1,1\le j\le N\)

上式中a_ij指的就是状态i向状态j转移的概率。有了α的初始值和递归公式，我们就可以计算出α_T(i)来，而我们想要的结果P(O|λ)其实就是α_T(i)对所有i求和。

而向后算法思路其实很像，引进变量β_t(i)，定义为：

β_t(i)=P(o_t+1o_t+2o_t+3…o_T|q_t= i,λ)

显然意义就是在给定HMM模型λ，在t时刻状态是i的情况下，后面的观测序列分别是o_t+1o_t+2o_t+3…o_T的概率。

向后算法是从后往前算的，我们知道β_T(i)=1，然后由βt+1(i)计算βt(i)得方法如下：

\(\beta_t(i)=\sum_{j=1}^{N}\beta_{t+1}(j)a_{ij}b_j(o_{t+1})1\le t\le T-1,1\le j\le N\)

最后要计算的结果P(O|λ)如下：

\(P(O|\lambda)=\sum_{i=1}^N\pi_ib_i(o_1)\beta_1(i)\)

解码问题其实就是求解状态转换序列q，使得P(O,q|λ)可以达到最大。同样，我们可以选择暴力计算所有的转化路径，然后找出最大的那个，还是同理，我们不推荐。

通信原理里面对卷积码的解码用到的是维特比算法(Viterbi Algorithm)，我们知道卷积码每个码都是跟前面的若干个有关的，而我们在解码的过程中寻求的就是一条使得误码率最低的路径来解决的。和我们这个问题很像，不是么？

好吧，我懒得慢慢敲了，流程如下，大家自己YY去吧。。

好吧，还是瞎扯两句，维特比译码的思想就是我们要计算全局最优路径之前，我们先计算局部最优路径。什么叫做局部最优路径呢？就是说我们计算出从0时刻开始到在t时刻出现状态i的最有可能的路径，这个概率就是我们上式中的δ_t(i)，含义就是序列中在t时刻以状态i终止，并且观测结果为o_t的概率。

有了δ_t(i)，我们在计算δ_t+1(j)的时候，就可以把在t时刻各个状态向t+1时刻状态j转换的最大概率给挑出来就可以了。计算方法就是上图第2步第一个计算公式，就是寻找出转换过来后概率最大的那个概率，然后乘以j状态下出现o_t+1的概率。好吧，有点拗口，看公式理解就好。。

有次基础就可以理解第一部初始化的意义了吧。。

然后呢，上图中还有另外一个变量Ψ，我们虽然可以计算出出现最后结果的那条路径的概率，但是，我们不知道那条路径是啥呀！！这才是我们真正要的，概率那个顶个鬼用。所以呢，我们在计算局部路径最大概率的时候就要把到达这里的最优路径也给顺便记录下来。看一下公式就能懂得，这里就不解释了。。

以上就是维特比译码求解最佳状态转换路径的方法了~

还有一个问题，额，我有一个很好的方法，但是这张纸不够了(哥德巴赫和费马从棺材里跳出来怒煽我两巴掌！！),好吧，其实是这样的，研究了一下那个方法，就是传说中的Baum-Welch Algorithm，然后我自己设定一个HMM模型，输出一串观测序列，再把这个序列送进这个Baum-Welch Algorithm里面训练，但是出来的结果不怎么好，所以容我有空再研究研究，出了好结果在来这里填坑吧。。

下面是向前向后算法，还有维特比译码的Matlab代码，如果读懂算法但是写代码有难度的孩子可以看看。但是那个代码用了很多矩阵相乘的运算，所以可读性不是那么高，但是尽量不用for语句写代码是我的原则，so…

code:

function HMM_Viterbi()
clc;
clear;
A = [0.5  0.375 0.125
     0.25 0.125 0.625
     0.25 0.375 0.375]';
B = [0.6  0.25 0.05
     0.2  0.25 0.1
     0.15 0.25 0.35
     0.05 0.25 0.5];
pai = [0.63 0.17 0.2];
output = [1 2 4 1 3 2];
N = size(A,1);
T = size(output,2); 

%%向前算法
alpha = zeros(T,N);
alpha(1,:) = pai .* B(output(1),:);
for i = 2:T
    alpha(i,:) = (alpha(i-1,:) * A) .* B(output(i),:);
end
% alpha'
P1 = sum(alpha(T,:))

%%向后算法
beta = zeros(T,N);
beta(T,:) = ones(1,N);
for i = T-1:-1:1
    beta(i,:) = (A * (beta(i+1,:) .* B(output(i+1),:))')';
end
% beta'
P2 = sum(pai .* B(output(1),:) .* beta(1,:))

%%viterbi
dirac = zeros(T,N);
psai = zeros(T,N);
dirac(1,:) = pai .* B(output(1),:);
for i = 2 : T
    dirac(i,:) = max([dirac(i-1,:)' dirac(i-1,:)' dirac(i-1,:)'] .* A).*B(output(i),:);
    [temp psai(i,:)] = max([dirac(i-1,:)' dirac(i-1,:)' dirac(i-1,:)'] .* A,[],1);
end
% dirac'
% psai'
[dir_max dir_max_idx] = max(dirac');

final_path = zeros(1,T);
for i = 1 : T-1
    final_path(i) = psai(i+1,dir_max_idx(i+1));
end
final_path(T) = dir_max_idx(T);

final_path

result:

P1 =

    0.0047

P2 =

    0.0047

final_path =

     1     3     3     2

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com