#include<iostream>
using namespace std;
int test[1111],train[1111];
int i,j,k,t,n,m,T;
double pi[10],a[10][10],b[10][10],a1[1111][10],b1[1111][10],e[1111][10][10],r[1111][10];
double tmp,tmp1,tmp2,ans;
int main(){
//训练数据T为100，状态数n为2，产生值数m为4，初始分布pi和每个状态产生值矩阵b概率随机，状态间转移概率矩阵a为[0.6,0.4;0.3,0.7]
T=100;m=4;n=2;
for(i=1;i<=n;i++)pi[i]=1.0/n;
for(i=1;i<=n;i++)for(j=1;j<=m;j++)b[i][j]=1.0/m;
a[1][1]=0.6;a[1][2]=0.4;a[2][1]=0.3;a[2][2]=0.7;
//输入训练数据
freopen("input1.txt","r",stdin);
for(i=1;i<=T;i++)cin>>train[i];
//迭代1000次后退出
for(k=1;k<=1000;k++){
for(j=1;j<=n;j++)a1[1][j]=pi[j]*b[j][train[1]];
//计算forward函数a1
for(t=1;t<T;t++)
for(j=1;j<=n;j++){
tmp=0;
for(i=1;i<=n;i++)tmp+=a1[t][j]*a[i][j];
a1[t+1][j]=b[j][train[t+1]]*tmp;
}
//计算backward函数b1
for(j=1;j<=n;j++)b1[T][j]=1;
for(t=T-1;t>0;t--)
for(i=1;i<=n;i++){
tmp=0;
for(j=1;j<=n;j++)tmp+=b1[t+1][j]*a[i][j]*b[j][train[t+1]];
b1[t][i]=tmp;
}
//e[t][i][j]表示在已知观测序列和模型参数时，在时刻t处于状态i，时刻t+1处于状态j的概率
for(t=1;t<T;t++){
tmp=0;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)tmp+=a1[t][i]*a[i][j]*b1[t+1][j]*b[j][train[t+1]];
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)e[t][i][j]=a1[t][i]*a[i][j]*b1[t+1][j]*b[j][train[t+1]]/tmp;
}
//r[t][i]表示在已知观测序列和模型参数时，在时刻t处于状态i的概率
for(t=1;t<T;t++)
for(i=1;i<=m;i++){
tmp=0;
for(j=1;j<=n;j++)tmp+=e[t][i][j];
r[t][i]=tmp;
}
//更新初始分布pi
for(i=1;i<=n;i++)pi[i]=r[1][i];
//更新状态转移概率矩阵a
for(i=1;i<=n;i++)
for(j=1;j<=n;j++){
tmp1=0;tmp2=0;
for(t=1;t<T;t++){
tmp1+=e[t][i][j];
tmp2+=r[t][i];
}
a[i][j]=tmp1/tmp2;
}
//更新状态产生值的概率矩阵b
for(i=1;i<=n;i++)
for(j=1;j<=m;j++){
tmp1=0;tmp2=0;
for(t=1;t<T;t++){
if(train[t]==j)tmp1+=r[t][i];
tmp2+=r[t][i];
}
b[i][j]=tmp1/tmp2;
}
}
//输入测试数据
freopen("input2.txt","r",stdin);
for(i=1;i<=T;i++)cin>>test[i];
//用forward函数计算在当前模型下得到观测序列的概率ans
for(j=1;j<=n;j++)a1[1][j]=pi[j]*b[j][test[1]];
for(t=1;t<T;t++)
for(j=1;j<=n;j++){
tmp=0;
for(i=1;i<=n;i++)tmp+=a1[t][j]*a[i][j];
a1[t+1][j]=b[j][test[t+1]]*tmp;
}
for(i=1;i<=n;i++)ans+=a1[T][i];
//输出最终状态转移矩阵a，状态产生值矩阵b，模型产生测试数据的概率ans
freopen("output.txt","w",stdout);
cout<<"状态转移矩阵为："<<endl;
for(i=1;i<=n;i++){for(j=1;j<=n;j++)cout<<a[i][j]<<" ";cout<<endl;}
cout<<"状态产生值矩阵为："<<endl;
for(i=1;i<=n;i++){for(j=1;j<=m;j++)cout<<b[i][j]<<" ";cout<<endl;}
cout<<"模型产生测试数据的概率为：\n"<<ans<<endl;
return 0;
}
