Forward-Backward算法做HMM的Inference

马尔可夫链的三个基本问题

1.已知一个序列，他的likelihood是什么样的。（使用前向算法）
2.求解一个最好的状态链（Viterbi算法）
3.优化或者重新估计这个HMM，比如重新估计发散矩阵和转换矩阵（Baum-Welch算法）

F/B算法是一个求解HMM的重要算法，它是动态规划（Dynamic programming, 这个名字的翻译有点意思）的重要一种。

F/B前提：假设发散概率矩阵（Emission Probablity Matrix）, 转换矩阵（Transition Probabilty Matrix) , 和初始概率（Initial probability）已知。

F/B目的：求在已知的观测数据X下的某个状态Z的概率。

F/B为什么要分成两部分：前向算法（Forward）和后向算法（Backward）

首先来看我们要求的概率p(Zk|X):

p(Zk|X)∝p(Zk,X)=p(Xk+1:n|Zk,X1:k)p(Zk,X1:k)(1)

在观测数据X固定的情况下，p(Zk|X)会等比例于p(Zk,X)。
根据Chain Rule，可以将p(Zk,X)分解，这里n的意思是X的数量。

根据图论的D-separation,我们可以将p(Xk+1:n|Zk,X1:k)中的X1:k去掉。因为如果我们依赖于Zk的话，Zk跟X1:k是条件不相关的，所以可以删掉。
那么p(Zk|X)就变成了：

p(Zk|X)∝p(Zk,X1:k)p(Xk+1:n|Zk)(2)

(2)式中右边的前半部分就是前向算法，后半部分是后向算法。

前向算法（Forward）

前向算法可以计算给定了观测数据后，这些数据跟目前HMM的参数的相似程度。或者说来计算观测数据在这些参数下的似然（likelihood）。前向算法要求p(Zk,X1:k),这里要明确的是状态的序号是跟观测序列的结尾是同步的。这个过程可以叫做filtering，就是在求这个序列的后验分布。我们可以用D-separation来继续化解。

p(Zk,X1:k)=∑Zk−1=1mp(Zk,Zk−1,X1:k)

m代表了状态Z的数量，根据边际概率，我们得到了这个式子。然后根据Chain rule:

=∑Zk−1=1mp(Xk|Zk,Zk−1,X1:k−1)p(Zk|Zk−1,X1:k−1)p(Zk−1,X1:k−1)

使用D-separation

=∑Zk−1=1mp(Xk|Zk)p(Zk|Zk−1)p(Zk−1,X1:k−1)

这个式子可以解释成为对于每个Zk来说先计算它的发散概率和转换概率乘以上一个Zk的对应的式子，然后把它们全部加起来。这样就构成了一个递归的式子。动态规划的最基本的一个思想就是要用递归的方法来解决问题。下面贴一段C语言的前向算法的代码：

typedef struct  {  int N; /* 隐藏状态数目;Q={1,2,…,N} */  int M; /* 观察符号数目; V={1,2,…,M}*/  double **A; /* 状态转移矩阵A[1..N][1..N]. a[i][j] 是从t时刻状态i到t+1时刻状态j的转移概率 */  double **B; /* 混淆矩阵B[1..N][1..M]. b[j][k]在状态j时观察到符合k的概率。*/  double *pi; /* 初始向量pi[1..N]，pi[i] 是初始状态概率分布 */  } HMM;  前向算法程序示例如下：  /* 　函数参数说明： 　*phmm：已知的HMM模型；T：观察符号序列长度； 　*O：观察序列；**alpha：局部概率（到目前状态为止所有概率的和）；*pprob：最终的观察概率 */  void Forward(HMM *phmm, int T, int *O, double **alpha, double *pprob)  {  　　int i, j; 　　/* 状态索引 */  　　int t; 　　 /* 时间索引 */  　　double sum; /*求局部概率时的中间值 */  　　/* 1. 初始化：计算t=1时刻所有状态的局部概率： */  　　for (i = 1; i <= phmm->N; i++)  　　　　alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]];  　　  　　/* 2. 归纳：递归计算每个时间点，t=2，… ，T时的局部概率 */  　　for (t = 1; t < T; t++)  　　{  　　　　for (j = 1; j <= phmm->N; j++)  　　　　{  　　　　　　sum = 0.0;  　　　　　　for (i = 1; i <= phmm->N; i++) /*将所有元素加起来*/ 　　　　　　　　sum += alpha[t][i]* (phmm->A[i][j]);  　　　　　　alpha[t+1][j] = sum*(phmm->B[j][O[t+1]]);  　　　　}  　　}  　　/* 3. 终止：观察序列的概率等于T时刻所有局部概率之和*/  　　*pprob = 0.0;  　　for (i = 1; i <= phmm->N; i++)  　　　　*pprob += alpha[T][i];  }  

后向算法

后向算法的意义没有前向算法的那么强，但是是构成F/B算法的一部分。我们要来计算这一部分：

p(Xk+1:n|Zk)=∑zk+1mp(Xk+1:n,Zk+1|Zk)

根据chain rule分解：

=∑zk+1mp(Xk+2:n|Zk+1,Zk,Xk+1)p(Xk+1|Zk+1,Zk)p(Zk+1|Zk)

根据D-separation:

=∑zk+1mp(Xk+2:n|Zk+1)p(Xk+1|Zk+1)p(Zk+1|Zk)

此时可以看出来上面式子中的第一项可以作为递归项，后面的两项依次是发散概率和转换概率。
这样我们就可以进行后向算法的计算了。

前后向算法

我们根据上述两个算法可以计算p(Zk|X)，即给定转换矩阵，发散矩阵，观测数据和一系列状态Z的情况下，可以计算某个状态Zk的概率。