隐马尔科夫(HMM)模型 前向后向(Forward_backward) 维特比 (viterbi)

#include<iostream>
#include<stdio.h>
#include<cstring>
#include<string>
using namespace std;

#define NN 10  //表示模型中的状态数
#define MM 7  //每个状态中不同观察序号的个数
#define l 100   //观测值的长度
int N,M,ll;   //对应于NN MM LL的含义
double t[l][NN+2];  //l时刻处于状态j的概率  for_ward算法是用到的数组
double t2[l][NN+2];  //l时刻处于状态j的概率  back_ward算法中用到的数组
double a[NN+2][NN+2];//状态转移矩阵
double b[NN+2][MM+2];  //每个状态的观察值分布矩阵
double p[NN+2];;//初始状态
int o[l];  //观察值序列
//从文件读入数据初始化矩阵
void input()
{
    freopen("in4.txt","r",stdin);
	cin>>N;
	cout<<"转移矩阵是："<<endl;
	for(int i=1;i<=N;i++)
	{
	    for(int j=1;j<=N;j++)
		{	
			cin>>a[i][j];
	        cout<<a[i][j]<<" ";
		}cout<<endl;
	}
	cin>>M;
	cout<<"状态概率矩阵是："<<endl;
	for(int i=1;i<=N;i++)
	{
	    for(int j=1;j<=M;j++)
		{	
			cin>>b[i][j];
	        cout<<b[i][j]<<" ";
		}cout<<endl;
	}
	//初始值
	cout<<"在该有向图中初始状态是："<<endl;
	for(int i=1;i<=N;i++)
	{	
		cin>>p[i];
	    cout<<p[i]<<" ";
	}cout<<endl;

	string s;
	cin>>s;
	ll=s.size();
	cout<<"观察值是："<<endl;
	for(int i=1;i<=ll;i++)
	{
		o[i]=int(s[i-1]-'0');
		cout<<o[i]<<" ";
	}cout<<endl;
}

void forward()
{
	double res=0;

	for(int i=1;i<=N;i++)
	{
	    t[1][i]=p[i]*b[i][o[1]]; //在 1 时刻观察到 o[1] 的概率
    	cout<<"时刻 "<<1<<" 时，处于状态 "<<i<<" 并且观察值是 "<<o[1]<<" 的概率是："<<t[1][i]<<endl;
	}
	
	for(int i=2;i<=ll;i++)
	{
	    for(int j=1;j<=N;j++)
		{
			for(int k=1;k<=N;k++)
		    {
				t[i][j]+=t[i-1][k]*a[k][j];
			}
			t[i][j]=t[i][j]*b[j][o[i]]; 
			cout<<"时刻 "<<i<<" 时，处于状态 "<<j<<" 并且观察值是 "<<o[i]<<" 的概率是："<<t[i][j]<<endl;
		}
	}
	for(int i=0;i<=N;i++)
		res+=t[ll][i];
	cout<<"forward 算法求得 输出该序列的概率是："<<res<<endl;
}

void backward()
{
	double res2=0;
    for(int i=1;i<=N;i++)
		t2[ll][i]=1;
	for(int i=ll-1;i>=1;i--)
	{
	    for(int j=1;j<=N;j++)
		{
			for(int k=1;k<=N;k++)
		    {
				t2[i][j]+=t2[i+1][k]*a[j][k]*b[k][o[i+1]]; 
			}
		}
	}
	for(int i=1;i<=N;i++)
		res2+=t2[1][i];
}

void for_back()
{
    double res3=0;
	for(int i=1;i<=N;i++)
	{
	    res3+=t[3][i]*t2[3][i];//这里的3 可以换成1-ll之间的任何一个数  表示某一个时刻 见decode()函数中 sum[i]的值总是相等的
	}
	cout<<"forward_backward 算法求得 输出该序列的概率是："<<res3<<endl;
}
void decode()  //看来一个资料说这个方法也能用于解码 但是好像不对
{
	double sum[l+1]={0},max=0,r[NN+2][l+1];
	for(int i=1;i<=ll;i++)
	{
	    for(int j=1;j<=N;j++)
		{
		    sum[i]+=t[i][j]*t2[i][j];
		}
	}
	cout<<endl;
    for(int i=1;i<=ll;i++)
	{   max=0;
		for(int j=1;j<=N;j++)
		{
			r[j][i]=double(t[i][j]*t2[i][j]*1.0/sum[j]*1.0);
			if(r[j][i]>max)
				max=r[j][i];
		}
		cout<<"解码问题的概率是："<<max<<endl;
	}
}

void viterbi()
{
	double max,theta[l+1][NN+2],w[l+1][NN+2];
	for(int i=1;i<=N;i++)
	{
	    theta[1][i]=p[i]*b[i][o[1]];
		w[1][i]=0; //w用来记录胜出的节点
	}

	int x;  //用来记录每一步的得胜节点
    for(int i=2;i<=ll;i++)
	{
		for(int j=1;j<=N;j++)
		{
			max=0;x=1;
			for(int k=1;k<=N;k++)
			{
		         theta[i][j]=theta[i-1][k]*a[k][j];
		         if(theta[i][j]>max)
				 {
					 max=theta[i][j];
			         x=k;
				 }
			}
			theta[i][j]=max*b[j][o[i]];
			w[i][j]=x;   //时刻i时 处于状态j  他的前一个的得胜节点是k	
		}
	}
	max=0;
	for(int i=1;i<=N;i++)
	{
	    if(theta[ll][i]>max)
		{
			max=theta[ll][i];
	        x=i;  //用来记录最终到达时间ll  即最后一个的得胜节点
		}
	}
	//回溯 反查找到每一步的得胜节点 即求得最后的路径
	cout<<"viterbi路径的反序是："<<x<<"  ";
	for(int i=ll-1;i>=1;i--)
	{
	    x=w[i+1][x];
		cout<<x<<"  ";
	}cout<<endl;
	cout<<"viterbi解码的概率是："<<max<<endl;
}

int main()
{
    input();
	forward();
	backward();
	for_back();
	//decode();
	viterbi();
	system("pause");
	while(1)
	{}
    return 0;

}