关于Vertebi算法的理解以及程序实现

HMM解码问题

      给定一个观察序列O=O₁O₂...O_T,和模型μ=（A,B,π），如何快速有效地选择在一定意义下“最优”的状态序列Q=q₁q₂...q_T，使该状态最好地解释观察序列。

           

      一种想法是求出每个状态的概率r_t(i)最大(r_t(i)=P(qt=si,O|μ）)，记q'_t(i)=arg_Qmax(r_t(i))，但是这样做，忽略了状态之间的关系，很可能两个状态之间的概率为0，即a_{q't(i)q't+1(i)}=0,这样求得的“最优”状态序列是不合法的。

      为防止状态之间转移概率为0（断续问题），换一种思路，不是求单个状态求得最大值，而是求得整个状态序列最大值，即求

                                   Q'= arg_QmaxP(Q|O,μ）

      此时用维特比算法，先定义下维特比变量δt(i):在时间t，HMM沿着一条路径到达状态si，并输出观察序列O=O₁O₂...O_t的最大概率:

　　      δ_t(i)=max P(q₁q₂...q_t=s_i,O₁O₂...O_t|μ)

           

                            t                      t+1

      上图中，对于从t时刻三个到 t+1时刻的状态1，到底取状态1,2还是3，不是看单独状态1,2还是3的概率，而是看在状态1,2,3各自的维特比变量值乘以相应的状态转换概率，从中选出最大值，假设2时最大，那么记下t+1时刻状态1之前的路径是t时刻的状态2，以此类推。

      δ_t(i)的递归关系式: δ_t+1(i)=max_j δ_t(j)*a_ji*b_i(O_t+1),为了记忆路径，定义路径变量ψ_t(i)，记录该路径上的状态s_i的前一个状态。

      维特比算法：

      step1 初始化：

              δ_t(i) = π_i*b_i(O₁), 1≤i≤N

              ψ_t(i) = 0

      step2 归纳计算：

　　　　  δ_t(i)=max_1≤j≤N δ_t-1(j)*a_ji*b_i(O_t),2≤t≤T;1≤i≤N

             记忆路径   ψ_t(i) = arg [max_1≤j≤Nδ_t-1(j)*a_ji*b_i(O_t)]

      step3 终结:

            Q_T' = arg max_1≤i≤N [δ_T(i)]

            P'(Q_T') = max_1≤i≤N [δ_T(i)]

　　 step4 路径回溯:

             q_t'=ψ_t+1(q_t+1') , t=T-1,T-2...1

      时间复杂度

      计算某时刻的某个状态的前向变量需要比较前一时刻的N个状态，此时时间复杂度为O(N),每个时刻有N个状态，此时时间复杂度为N*O(N)=O(N²),又有T个时刻，所以时间复杂度为T*O(N²)=O(N²T)。

      程序例证

             

       step1 初始化：δ₁(1) = 0.2*0.5=0.1 ，δ₁(2) = 0.4*0.4=0.16， δ₁(3) = 0.4*0.7=0.21

       step2 归纳计算：δ₂(1) =max[0.1*0.5,0.16*0.3,0.21*0.2]*0.6

       ...

      step3 终结：最佳路径是δ₄(1)δ₄(2)δ₄(3)最大的一个对应的状态

      step4 回溯：从最后一个状态往回返

     程序代码

package Vertebi;


public class vertebi {


/**
* @param args
*/
public static void main(String[] args) {
// TODO Auto-generated method stub
vertebi vb=new vertebi();
}

public vertebi(){
double[][] a=new double[][]{{0.5,0.2,0.3},{0.3,0.5,0.2},{0.2,0.3,0.5}};
double[][] b=new double[][]{{0.5,0.5},{0.4,0.6},{0.7,0.3}};
double[] Pi=new double[]{0.2,0.4,0.4};
double[][]deta=new double[4][4];
int[] list=new int[]{0,1,0,1};
int count=1;
int[] record=new int[4];
//初始化
deta[0][0]=0.2*0.5;
deta[0][1]=0.4*0.4;
deta[0][2]=0.4*0.7;
deta[0][3]=0;
//初始化得出最大概率
deta[0][3]=0.28;
record[0]=2;

//迭代
for(int i=1;i<=3;i++){
for(int j=0;j<=2;j++){
deta[i][j]=0;
deta[i][j]+=deta[i-1][record[i-1]]*a[record[i-1]][j]*b[j][list[count]];
}
for(int m=2;m>0;m--){
if(deta[i][m]<deta[i][m-1])
{deta[i][3]=deta[i][m-1];
record[i]=m-1;}
else {
deta[i][3]=deta[i][m];
record[i]=m;
}
}
count+=1;
}

//比较大小
for(int k=0;k<=3;k++){
System.out.println(deta[k][3]+"  "+record[k]);
}
}
}

0.28  2
0.05039999999999999  1
0.010079999999999999  1
0.0030239999999999993  1