动态规划之隐含马尔可夫模型(HMM)和维特比算法(Viterbi Algorithm)

动态规划之(HMM)和(Viterbi Algorithm)

1. 实际问题

HMM－韦小宝的骰子
• 两种骰子，开始以2/5的概率出千。
– 正常A：以1/6的概率出现每个点
– 不正常B： 5,6出现概率为3/10,其它为1/10
• 出千的随机规律

观测到其一次投掷结果

Y={1,3,4,5,5,6,6,3,2,6}

问题是韦小宝何时出千了？

2. 马尔可夫模型(HMM)介绍

隐马尔可夫模型（Hidden Markov Model，HMM）作为一种统计分析模型，创立于20世纪70年代。80年代得到了传播和发展，成为信号处理的一个重要方向，现已成功地用于语音识别，行为识别，文字识别以及故障诊断等领域。模型如下图所示：

隐马尔可夫模型是马尔可夫链的一种，它的状态不能直接观察到，但能通过观测向量序列观察到，每个观测向量都是通过某些概率密度分布表现为各种状态，每一个观测向量是由一个具有相应概率密度分布的状态序列产生。所以，隐马尔可夫模型是一个双重随机过程—-具有一定状态数的隐马尔可夫链和显示随机函数集。自20世纪80年代以来，HMM被应用于语音识别，取得重大成功。到了90年代，HMM还被引入计算机文字识别和移动通信核心技术“多用户的检测”。HMM在生物信息科学、故障诊断等领域也开始得到应用。
隐马尔可夫模型（HMM）可以用五个元素来描述，包括2个状态集合和3个概率矩阵：

2.1 HMM数学模型表达

(1)隐过程为X={X1,X2,X3,...,XT}
(2)观察过程为Y={Y1,Y2,Y3,...,YT}
(3)模型参数{π,A,B}

2.1.1 HMM模型各参数详细说明

(1). 隐含状态 S
这些状态之间满足马尔可夫性质，是马尔可夫模型中实际所隐含的状态。这些状态通常无法通过直接观测而得到。（例如S1、S2、S3等等)
(2). 可观测状态Y
在模型中与隐含状态相关联，可通过直接观测而得到。(例如Y1,Y2,Y3,...,YT等等，可观测状态的数目不一定要和隐含状态的数目一致）
(3). 初始状态概率矩阵π
表示隐含状态在初始时刻t=1或t=0的概率矩阵，(例如t=1时，P(S1)=p1、P(S2)=P2、P(S3)=p3，则初始状态概率矩阵 π=[ p1 p2 p3 ].
(4). 隐含状态转移概率矩阵 A。
描述了HMM模型中各个状态之间的转移概率。
其中Aij=P(Sj|Si),1≤i,,j≤N.表示在 t 时刻、状态为 Si的条件下，在t+1时刻状态是 Sj的概率。
(5). 观测状态转移概率矩阵 B （英文名为Confusion Matrix，直译为混淆矩阵不太易于从字面理解）。
令N代表隐含状态数目，M代表可观测状态数目，则：
Bij=P(Yi|Sj),1≤i≤M,1≤j≤N.
表示在 t 时刻、隐含状态是 Sj 条件下，观察状态为 Oi 的概率。
总结：一般的，可以用λ=(A,B,π)三元组来简洁的表示一个隐马尔可夫模型。隐马尔可夫模型实际上是标准马尔可夫模型的扩展，添加了可观测状态集合和这些状态与隐含状态之间的概率关系。

3. Viterbi算法介绍

维特比算法是一种动态规划算法用于寻找最有可能产生观测事件序列的-维特比路径-隐含状态序列,特别是在马尔可夫信息源上下文和隐马尔可夫模型中。术语“维特比路径”和“维特比算法”也被用于寻找观察结果最有可能解释相关的动态规划算法。例如在统计句法分析中动态规划算法可以被用于发现最可能的上下文无关的派生(解析)的字符串，有时被称为“维特比分析”。维特比算法由安德鲁·维特比(Andrew Viterbi)于1967年提出，用于在数字通信链路中解卷积以消除噪音。 此算法被广泛应用于CDMA和GSM数字蜂窝网络、拨号调制解调器、卫星、深空通信和802.11无线网络中解卷积码。现今也被常常用于语音识别、关键字识别、计算语言学和生物信息学中。例如在语音(语音识别)中，声音信号作为观察到的事件序列,而文本字符串,被看作是隐含的产生声音信号的原因，因此可对声音信号应用维特比算法寻找最有可能的文本字符串。

维特比算法的基础可以概括成下面三点：
(1)如果概率最大的路径p(或者说最短路径)经过某个点，比如途中的X22，那么这条路径上的起始点S到X22的这段子路径Q，一定是S到X22之间的最短路径。否则，用S到X22的最短路径R替代Q，便构成一条比P更短的路径，这显然是矛盾的。证明了满足最优性原理。
(2)从S到E的路径必定经过第i个时刻的某个状态，假定第i个时刻有k个状态，那么如果记录了从S到底i个状态的所有k个节点的最短路径，最终的最短路径必经过其中一条，这样，在任意时刻，只要考虑非常有限的最短路即可。
(3)结合以上两点，假定当我们从状态i进入状态i+1时，从S到状态i上各个节的最短路径已经找到，并且记录在这些节点上，那么在计算从起点S到第i+1状态的某个节点Xi+!的最短路径时，只要考虑从S到前一个状态i所有的k个节点的最短路径，以及从这个节点到Xi+1，j的距离即可。

3.1 Viterbi算法伪代码描述

现在我们回到问题：
构建模型三个参数如下：

接下来进行用Viterbi算法求解，实现过程如下：

/********************************************************Description: Hidden Markov Model(HMM) and Viterbi's decoding Algorithm Author:Robert.TianyiDate:2016.10.28********************************************************/#include<iostream>#include<iomanip>using namespace std;#define NumState 2  //状态数量 A B #define NumViewValue  6//观测值种类数 #define ViewSequnceLength  10//观测到的序列长度 void Viterbi_Algorithm(float A[][NumState],float B[][NumViewValue],float Pi[],int Y[]); int main(){    float A[NumState][NumState]={{0.8,0.2},                                 {0.1,0.9}};//状态转移矩阵    cout<<"输出转移矩阵A:"<<endl;    for(int i=0;i<NumState;i++){        for(int j=0;j<NumState;j++){            cout<<setw(4)<<A[i][j]<<" ";        }        cout<<endl;    }    float B[NumState][NumViewValue]={{1.0/6,1.0/6,1.0/6,1.0/6,1.0/6,1.0/6,},                                    {1.0/10,1.0/10,1.0/10,1.0/10,3.0/10,3.0/10,}};//B矩阵     cout<<"输出矩阵B:"<<endl;    for(int i=0;i<NumState;i++){        for(int j=0;j<NumViewValue;j++){            cout<<setw(10)<<B[i][j]<<" ";        }        cout<<endl;    }    float Pi[NumState]={0.6,0.4};//初始各状态分布矩阵    cout<<"输出初始状态分布："<<endl;    for(int i=0;i<NumState;i++)        cout<<setw(4)<<Pi[i]<<" ";    cout<<endl;     int Y[ViewSequnceLength]={1,3,4,5,5,6,6,3,2,6};//观测序列     Viterbi_Algorithm(A, B,Pi,Y);    return 0;}void Viterbi_Algorithm(float A[][NumState],float B[][NumViewValue],float Pi[],int Y[]){    float delta[ViewSequnceLength][NumState];    delta[0][0]=Pi[0]*B[0][0];    delta[0][1]=Pi[1]*B[1][0];    int path[ViewSequnceLength][NumState];    path[0][0]=0;path[0][1]=0;    //迭代     for(int t=1;t<ViewSequnceLength;t++){        for(int j=0;j<NumState;j++){            float tempMax=0.0;            int flag_i=0;            for(int i=0;i<NumState;i++){                if(delta[t-1][i]*A[i][j]>tempMax){                    tempMax=delta[t-1][i]*A[i][j];                    flag_i=i;                }            }            delta[t][j]=tempMax*B[j][Y[t]-1];            path[t][j]=flag_i;          }    }     //终止 求出最大概率     int X[ViewSequnceLength];    X[ViewSequnceLength-1]=0;    float p=delta[ViewSequnceLength-1][0];    for(int i=1;i<NumState;i++){        if(delta[ViewSequnceLength-1][i]>p){                p=delta[ViewSequnceLength-1][i];                X[ViewSequnceLength-1]=i;        }       }    cout<<"Max_probablity="<<p<<endl;    //回溯    for(int t=ViewSequnceLength-2;t>=0;t--){        X[t]=path[t+1][X[t+1]];    }     //信息输出     cout<<endl;    float result;    cout<<"各时刻具体情况："<<endl;    for(int t=0;t<ViewSequnceLength;t++){        cout<<"时刻t="<<t<<" ";        for(int j=0;j<NumState;j++){                cout<<setw(18)<<delta[t][j]<<" ";            }            cout<<endl;    }    cout<<"输出隐含过程："<<endl;     for(int t=0;t<ViewSequnceLength;t++){            cout<<setw(4)<<X[t]<<" ";    }    char charX[ViewSequnceLength];    for(int t=0;t<ViewSequnceLength;t++){        if(X[t]==0)            charX[t]='A';        else             charX[t]='B';    }    cout<<endl;    for(int t=0;t<ViewSequnceLength;t++){            cout<<setw(4)<<charX[t]<<" ";    }    cout<<endl<<"输出观测序列："<<endl;    for(int t=0;t<ViewSequnceLength;t++){            cout<<setw(4)<<Y[t]<<" ";    }}

运行结果如下：

更形象展示结果：

我们可以知道韦小宝在第四次出千啦！好狡猾的！！！