上学期学习了HMM模型,自认为弄明白,嘴上也时不时的也蹦出个HMM,但细细想来竟理不出个条条框框,禁不住心慌起来,于是又拿起HMM的资料再次揣摩,希望能够将其真正的转化为自己的知识。(公式都没有显示,稍后解决)
1、 前向算法
前向算法用于HMM的模型评估问题。假如我们已知HMM的三元组(初始状态概率向量,状态转移矩阵,观测概率矩阵)和观测序列,计算模型产生观测的概率,即
计算所求最直接的思路便是,首先列举出所有可能的长度为T的状态序列,然后求各个状态序列I产生观测序列O的联合概率,最后再对计算出的所有联合概率求和得到结果。假如共有N种状态,那么对于观测序列O,共有N的T次方种路径。这种时间复杂度为指数级的思想显然不可取。对该过程分析可知,该过程之所以产生指数级的时间复杂度是由于做了大量重复的计算,比如N条路径前t-1个时刻所经历的状态完全相同,仅T时刻的状态不同,那么前t-1时刻的概率就重复计算了N次。一个很好的解决办法就是记住已经计算过的路径的概率这就是前向算法的核心。
结合图立体的赶脚一下,横向代表时刻t,纵向代表所有状态I:
于是我们采用前向算法来计算;
前向算法的思想:
(1)初值 t=1,观测为,
注释:在时刻t=1,由初始状态概率矩阵,得到达每个状态的概率Pi(Ii),然后再由Pi乘以P(O1|Ii),得到由各个状态产生观测O1的概率。我们将计算出的结果放在一个N*T的二维数组M中,为第一列的元素。
(2)递推 t = 2,3,..,T
注释:我们要计算时刻t观测为Ot的概率,第一步,需求时刻t到达状态Ii的概率P(Ii|t),第二步,再乘以由状态Ii产生观测Ot的概率。由t=1时刻达到t时刻的状态Ii,有N1*N2*…Nt中可能的路径,然而在我们的二维数组M中的t-1列已将保存了前N1*N2*…Nt-1种路径倒到t-1的概率,因此我们在计算第一步时只需在t-1列的基础上,求出到达li的N种路径的概率之和即。第二步,由第一步的结果乘以由状态Ii产生观测Ot的概率即可。
(3)
注释:我们把M数组的最后一列相加即为所求。
C++代码如下:VS2012中简单测试过。
#include <iostream>
using namespace std;
typedef int StateT;
typedef int ObserveT;
//求一维数组中元素的个数 ,也可做求二维数组的行数。
//不适用于数组指针
template <typename T>
int getArrayLen(T &arr)
{
return sizeof(arr)/sizeof(arr[0]);
}
//前向算法计算给定HMM的条件下产生某观测observeset的概率
void Forward(
int N, //状态个数
int T, //观测的个数
int *observeset,//观测集合 一维数组 T 下标以为时刻t, 值为该观测在观测集合中的序号,用于观测概率的获取
double *initial, //初始概率向量 一维数组 N
double **Transfer, //状态转移概率矩阵 二维数组 N*N
double ** observer, //观测概率矩阵 N*T
double **p //存储计算的二维矩阵 N(状态个数)*T(观测序列长度)
){
for (int t = 0; t < T; t++)//T个观测,循环计算T次
{
for (int j = 0; j < N; j++)//时刻t,处于各个状态时,产生观测Ot的概率
{ //即p[j][t]=jointP*observer[j][Ot]
if (t== 0) //第一步:在初始时刻状态i的概率,由初始状态概率向量得到
{
p[j][t] = initial[j]*( observer[j][observeset[t]]);
}
else//第二步 开始递推求结果矩阵
{
double jointP = 0.0;
for (int k = 0; k < N; k++)
{//计算t时刻处于状态j的概率 = t-1时刻从N个状态转移到状态j的概率之和 jointP
jointP += p[k][t-1]*Transfer[k][j];
}
//t时刻处于状态j的条件下,产生观测Ot的概率
p[j][t] = jointP * observer[j][observeset[t]];
}
}
}
//第三步
double fp = 0.0;
for (int i = 0; i < N; i++)
{
fp += p[i][T-1];
}
cout << "该模型产生观测O的概率为:"<< fp << endl;
}
//测试
void main()
{
int N = 3; //状态个数
//HMM 三元组
//初始概率
double pai[] = {0.2,0.4,0.4};
//状态转移概率矩阵
double A[3][3] = {{0.5,0.2,0.3},{0.3,0.5,0.2},{0.2,0.3,0.5}};
double *pA[3];
pA[0] = A[0];
pA[1] = A[1];
pA[2] = A[2];
//观测概率矩阵
double B[3][2] = {{0.5,0.5},{0.4,0.6},{0.7,0.3}};
double *pB[3];// 为满足形参 double** observe, 我们需将B 转换为pB这种形式传递
pB[0] = B[0];
pB[1] = B[1];
pB[2] = B[2];
//已知的观测状态序列
int O[3] = {0,1,0};//0 红球 1白球
int T = getArrayLen(O);
double **M = new double*[N];//用来存放计算的概率矩阵N*T
for (int i = 0; i < N; i++)
{
M[i] = new double[T];
}
Forward(N,T,O,pai,(double**)pA,(double**)pB,M);
}
参考文献:
HMM学习最佳范例------- 来自52nlp
《统计学习方法》------李航著