隐马尔科夫算法通俗解释.doc

HMM介绍(通俗易懂版)
Hidden Markov Models是一种统计信号处理方法,模型中包含2个序列和3个矩阵:状态序列S、观察序列O、初始状态矩阵P、状态转移矩阵A、混淆矩阵B。举个例子来说明。
Example1：
你一个异地的朋友只做三种活动：散步、看书、做清洁。每天只做一种活动。假设天气只有两种状态：晴和兩。每天只有一种天气。你的朋友每天告诉你他做了什么，但是不告诉你他那里的天气。
某一周从周一到周五每天的活动分别是{读书，做清洁，散步，做清洁，散步}----这就是观察序列O，因为你可以观察得到。
从周一到周五的天气依次是{晴，兩，晴，晴，晴}----这就是状态序列S，状态序列是隐藏的，你不知道。
根据长期统计，某天晴的概率是0.6，兩的概率是0.4。则。
从晴转晴的概率是0.7，从晴转兩的概率是0.3。从兩转晴的概率是0.4，从兩转兩的概率是0.6。则。
天气晴时，散步的概率是0.4，看书的概率是0.3，做清洁的概率是0.3。天气兩时，散步的概率是0.1，看书的概率是0.4，做清洁的概率是0.5。则。
该模型和实际情况有明显不符的地方：用一简单的状态转移矩阵A来表示状态的转移概率的前提是t时刻的状态只跟t-1时刻的状态有关，而实际上今天的天气跟过去几天的天气都有关系，而且跟过去几天的晴朗程度、兩量大小都有关系；混淆矩阵B认为今天的活动只跟今天的天气有关系，实际上今天的活动跟过去几天的活动也有关系，比如过去一周都没有打扫房间，那今天做清洁的概率就大大增加。
模型介绍完了。
评估问题
隐马尔可夫模型中包含一个评估问题：已知模型参数，计算某一特定输出序列的概率。通常使用forward算法解决。 
比如计算活动序列{读书，做清洁，散步，做清洁，散步}出现的概率，就属于评估问题。
如果穷举的话，观察序列会有2^5种，需要分别计算它们出现的概率，然后找出概率最大的。
穷举法中有很多重复计算，向前算法就是利用已有的结果，减少重复计算。
算法借助于一个矩阵Q[LEN][M]，其中M是所有状态的种数，Q[i][j]表示从第0天到第i天，满足观察序列，且第i天隐藏状态为Sj的所有可能的隐藏序列的概率之和。最终所求结果为Q[LEN-1][0]+...+Q[LEN-1][M-1]，即最后一天，所有隐藏状态下，分别满足观察序列的概率值之和。
比如Q[0][0]=p(第一天做卫生 且 第一天晴)=p(天晴)*p(做卫生|天晴)=P[0]*B[0][2]=0.6*0.3=0.18
Q[0][1]=p(第一天做卫生 且 第一天下雨)=p(下雨)*p(做卫生|下雨)=P[1]*B[1][2]=0.4*0.5=0.2
Q[1][0]=p(第一天做卫生 且 第二天晴 且 第二天做卫生)
=p(第一天做卫生 且 第二天晴)*p(天晴的情况下做卫生)
=p{ p(第一天做卫生 且 第一天晴)*p(从天晴转天晴)+p(第一天做卫生 且 第一天下雨)*p(从下雨转天晴) }*p(天晴的情况下做卫生)
={ Q[0][0]*A[0][0] + Q[0][1]*A[1][0] } * B[0][2]
Q[1][1]= ……
…… ……
可以看到计算Q矩阵的每i行时都用到了第i-1行的结果。
解码问题
解码问题是：已知模型参数，寻找最可能的能产生某一特定输出序列O(LEN)的隐含状态的序列。通常使用Viterbi算法解决。 
观察序列长度为LEN，则隐藏状态序列长度也是LEN，如果采用穷举法，就有M^LEN种可能的隐藏状态序列，我们要计算每一种隐藏状态到指定观察序列的概率，最终选择概率最大的。
穷举法中有很多重复计算，Viterbi算法就是利用已有的结果，减少重复计算。
跟评估问题非常相似，不同点在于评估算的是和，解码算的是最大值。
Viterbi算法主要就是在计算一个矩阵Q[LEN][M]，其中Q[i][j]表示从第0天到第i天，满足观察序列，且第i天隐藏状态为Sj的所有可能的隐藏序列的概率的最大值。另外还要建立一个矩阵Path[LEN][M]，用来记录状态序列中某一状态之前最可能的状态。
举个例子，假如指定观察序列是{读书，做卫生，散步，做卫生，散步}，求出现此观察序列最可能的状态序列是什么。
Q[0][0]=p(第一天读书 且 第一天晴)=p(天晴)*p(读书|天晴)
Path[0][0]=-1;
Q[0][1]=p(第一天读书 且 第一天下雨)=p(下雨)*p(读书|下雨)
Path[0][1]=-1;
关键是从第二天开始，Q[1][0]表示：满足“第一天读书 且 第二天做卫生 且 第二天晴”的所有可能的隐藏序列的概率的最大值。那么满足“第一天读书 且 第二做卫生 且 第二天晴”的所有可能的隐藏序列有哪些呢？第二天是必须满足晴天的，第二天之前的状态可以任意变。则所有可能的隐藏序列就是“晴  晴”和“雨  晴”。实际上考虑第三天（及第三天以后）时，并不需要考虑“所有”可能的隐藏序列，而只需要考虑第二天的不同状态取值，这是因为马氏过程有无后效性--tm时刻所处状态的概率只和tm-1时刻的状态有关，而与tm-1时刻之前的状态无关。
Q[1][0]=
max{ p(第一天晴 且 第一天读书 且 第二天晴 且第二天做卫生) ，p(第一天下雨 且 第一天读书 且 第二天晴 且 第二天做卫生) }
=max{ p(第一天读书 且 第一天晴)*p(天晴转天晴)，p(第一天读书 且 第一天下雨)*p(下雨转天晴) } * p(做卫生|天晴)
=max{ Q[0][0]*A[0][0]，Q[0][1]*A[1][0] } * B[0][2]
假如Q[0][0]*A[0][0] < Q[0][1]*A[1][0]，则Path[1][0]=1；假如Q[0][0]*A[0][0] > Q[0][1]*A[1][0]，则Path[1][0]=0。
Q[1][1]= ……
…… ……
可以看到计算Q矩阵的每i行时都用到了第i-1行的结果。

Example2：（l自行练习用）
这里设状态数目为Ｎ； 每个状态有Ｍ 个观测值； π＝（ π １， π ２，…， πＮ ） 为各状态初始向量； Ａ＝（ ａ ｉ ｊ） Ｎ×Ｎ为状态转移概率矩阵； Ｂ＝（ ｂ ｉ ｊ） Ｎ×Ｍ
为状态观测概率矩阵； Ｏ＝｛ ｏ １， ｏ ２，…， ｏ Ｔ｝ 为观测值序列； Ｑ＝｛ ｑ １， ｑ ２，…， ｑ Ｔ｝ 为状态序列。
要针对状态观测概率矩阵 ． 图７中 在最初始ＨＭＭ 中 和 都采用随机阵列表示 且
满足下式条件
（上传者ps。Bij也满足这个aij的这个条件） 
在假设状态种类数和观测种类数都为3的前提下。
O:a,c,b,b,c,b,a,c,b

根据上述的HMM中viterbi算法可得：
Q：1,2,2,2,2,2,2,2,2
通过结合O（观测）和Q（状态）矩阵可发现B中在状态1，2产生a，b，c的观测值的概率与初试矩阵不符。
O: a,c,b,b,c,b,a,c,b
Q：1, 2, 2, 2, 2, 2 ,2, 2,2
可见B中p（a/1）=1,p（b/1）=0, p（c/1）=0,
p（a/2）=1/8=0.125, p（b/2）=4/8=0.5, p（c/2）=3/8=0.375,
修正得：

将新的矩阵 B结合上述的 O A л.作为初始值.通过Ｂａｕｍ－Ｗｅｌｃｈ多序列训练算法并结合前向－后向算法判断的收敛性 即可得出模型参数 A Bл 的最终估计值.从原理上可看出 该方法要求样本的观测序列要尽量充分．后面的实验会发现在较充分的样本下 该训练方法得到使得对应的接近全局最大的比较精确的 ＨＭＭ ．

下面给出两种算法的Java代码：（尚未测试，请自行验证）
package hmm;

/**
* @author Orisun
* date 2011-10-20
*/
import java.util.ArrayList;
import java.util.Collections;

public class HMM {
ArrayList<String> state;
ArrayList<String> observation;
double[] P;
double[][] A;
double[][] B;
int M;
int N;

HMM() {
// 天气只有两种状态：晴和兩。每天只有一种天气。
state = new ArrayList<String>();
state.add("sunny");
state.add("rain");
M = state.size();
// 活动只有三种：散步、看书、做清洁。每天只做一种活动。
observation = new ArrayList<String>();
observation.add("walk");
observation.add("read");
observation.add("clear");
N = observation.size();
// 初始状态矩阵。某天晴的概率是0.6，兩的概率是0.4。
P = new double[] { 0.6, 0.4 };
// 状态转移矩阵。
A = new double[M][];
// 从晴转晴的概率是0.7，从晴转兩的概率是0.3
A[0] = new double[] { 0.7, 0.3 };
// 从兩转晴的概率是0.4，从兩转兩的概率是0.6
A[1] = new double[] { 0.4, 0.6 };
// 混淆矩阵
B = new double[M][];
// 天气晴时，散步的概率是0.4，看书的概率是0.3，做清洁的概率是0.3。
B[0] = new double[] { 0.4, 0.3, 0.3 };
// 天气兩时，散步的概率是0.1，看书的概率是0.4，做清洁的概率是0.5。
B[1] = new double[] { 0.1, 0.4, 0.5 };
}

public double forward(ArrayList<String> observe) {
double rect = 0.0;
int LEN = observe.size();
double[][] Q = new double[LEN][];
// 第一天计算，状态的初始概率，乘上隐藏状态到观察状态的条件概率。
Q[0] = new double[M];
for (int j = 0; j < M; j++) {
Q[0][j] = P[j] * B[j][observation.indexOf(observe.get(0))];
}
// 第一天以后的计算，首先从前一天的每个状态，转移到当前状态的概率求和，然后乘上隐藏状态到观察状态的条件概率。
for (int i = 1; i < LEN; i++) {
Q[i] = new double[M];
for (int j = 0; j < M; j++) {
double sum = 0.0;
for (int k = 0; k < M; k++) {
sum += Q[i - 1][k] * A[k][j];
}
Q[i][j] = sum * B[j][observation.indexOf(observe.get(i))];
}
}
for (int i = 0; i < M; i++)
rect += Q[LEN - 1][i];
return rect;
}

public ArrayList<String> viterbi(ArrayList<String> observe) {
ArrayList<String> sta = new ArrayList<String>();
int LEN = observe.size();
double[][] Q = new double[LEN][];
int[][] Path = new int[LEN][];
// 第一天计算，状态的初始概率，乘上隐藏状态到观察状态的条件概率。
Q[0] = new double[M];
Path[0] = new int[M];
for (int j = 0; j < M; j++) {
Q[0][j] = P[j] * B[j][observation.indexOf(observe.get(0))];
Path[0][j] = -1;
}
//第一天以后的计算，首先从前一天的每个状态，转移到当前状态的概率的最大值，然后乘上隐藏状态到观察状态的条件概率。
for (int i = 1; i < LEN; i++) {
Q[i] = new double[M];
Path[i] = new int[M];
for (int j = 0; j < M; j++) {
double max = 0.0;
int index = 0;
for (int k = 0; k < M; k++) {
if (Q[i - 1][k] * A[k][j] > max) {
max = Q[i - 1][k] * A[k][j];
index = k;
}
}
Q[i][j] = max * B[j][observation.indexOf(observe.get(i))];
Path[i][j] = index;
}
}
// 找到最后一天天气呈现哪种状态的概率最大
double max = 0;
int index = 0;
for (int i = 0; i < M; i++) {
if (Q[LEN - 1][i] > max) {
max = Q[LEN - 1][i];
index = i;
}
}
System.out.println("出现最可能的隐藏序列的概率是"+max);
sta.add(state.get(index));
// 动态规划，逆推回去各天出现什么天气概率最大
for (int i = LEN - 1; i > 0; i--) {
index = Path[i][index];
sta.add(state.get(index));
}
// 把状态序列再顺过来
Collections.reverse(sta);
return sta;
}

public static void main(String[] args) {
HMM hmm = new HMM();

ArrayList<String> obs_seq = new ArrayList<String>();
obs_seq.add("read");
obs_seq.add("clear");
obs_seq.add("walk");
obs_seq.add("clear");
obs_seq.add("walk");

double ratio = hmm.forward(obs_seq);
System.out.println(ratio);

ArrayList<String> sta_seq = hmm.viterbi(obs_seq);
System.out.println(sta_seq.toString());
}
}