隐马尔科夫模型

大致概念：隐马尔科夫模型可以表示为一个五元组（S,V,A,B,P）,S是状态集合，V是输出符号组成的集合，A是状态转移矩阵，B是输出符号在某状态下的概率分布，P是初始状态概率分布。

问题1：估算观察序列的概率。给定隐马尔科夫模型的A、B、P和观察序列O=(o1,o2……)，求出出现此观察序列的概率。下面的输入假设初始状态的概率分布默认为落在每个状态上的概率相同。

思路：动态规划。可以使用向前算法和向后算法。a、向前算法为：dp[i][j]表示在时刻i，处在状态j，并且部分观察序列符合(o1,o2...oi)的概率。

输入数据解释：

//数据组数
1

//状态数，输出符号数，观察序列长度
3 2 3

//以下是转移矩阵
0.9 0.05 0.05
0.45 0.1 0.45
0.45 0.45 0.1

//以下每行表示某个状态的不同输出的概率
0.5 0.5
0.75 0.25
0.25 0.75

//输出序列
1 1 2

#include <stdio.h>
#define N 50
double s[N][N],t[N][N];//s表示状态转移矩阵，t[i][j]表示状态为i输出为j的概率
int out[N];//输出序列
double dp[N][N];
int T,n,m,p;//n表示状态数，m表示输出符号数量,p表示输出序列的数量
int main(){
	freopen("a.txt","r",stdin);
	scanf("%d",&T);
	while(T--){
		int i,j,k;
		double res = 0;
		scanf("%d %d %d",&n,&m,&p);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=n;j++)
				scanf("%lf",&s[i][j]);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=m;j++)
				scanf("%lf",&t[i][j]);
		for(i = 1;i<=p;i++)
			scanf("%d",&out[i]);
		
		for(i = 1;i<=n;i++)//初始化
			dp[1][i] = 1./n * t[i][out[1]];
		for(i = 2;i<=p;i++)
			for(j = 1;j<=n;j++){
				double sum = 0;
				for(k = 1;k<=n;k++)//上一步在状态k，这一步转到状态j
					sum += dp[i-1][k] * s[k][j];
				dp[i][j] = sum * t[j][out[i]];
			}
		for(i = 1;i<=n;i++)
			res += dp[p][i];
		printf("%lf\n",res);
	}
	return 0;
}

b、后向算法。此时的dp[i][j]表示在时刻i，状态为j下，观察序列为{oi+1...on}的概率。

#include <stdio.h>
#define max(a,b) ((a)>(b)?(a):(b))
#define N 50
double s[N][N],t[N][N];//s表示状态转移矩阵，t[i][j]表示状态为i输出为j的概率
int out[N];//输出序列
double dp[N][N];
int T,n,m,p;//n表示状态数，m表示输出符号数量,p表示输出序列的数量
int main(){
	freopen("a.txt","r",stdin);
	scanf("%d",&T);
	while(T--){
		int i,j,k;
		double res = 0;
		scanf("%d %d %d",&n,&m,&p);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=n;j++)
				scanf("%lf",&s[i][j]);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=m;j++)
				scanf("%lf",&t[i][j]);
		for(i = 1;i<=p;i++)
			scanf("%d",&out[i]);
		
		for(i = 1;i<=n;i++)//初始化
			dp[p][i] = 1;
		for(i = p-1;i>=1;i--)
			for(j = 1;j<=n;j++){
				double sum = 0;
				for(k = 1;k<=n;k++)//上一步在状态k，这一步转到状态j
					sum += dp[i+1][k] * s[j][k] * t[k][out[i+1]];
				dp[i][j] = sum;
			}
		for(i = 1;i<=n;i++)
			dp[0][i] = 1./n*t[i][out[1]]*dp[1][i];
		for(i = 1;i<=n;i++)
			res += dp[0][i];
		printf("%lf\n",res);
	}
	return 0;
}

问题2：求解最佳状态转移序列。输入内容同问题1，只是问题变成寻找一个状态转换序列，使得该状态转换序列最有可能产生输入所示的观察序列。

思路：Viterbi算法。动态规划思路与问题1的前向算法如出一辙。区别在于1、中间的迭代过程由求和变为求最大值。2、需要一个二维数组path[i][j]来记录时刻i到达状态j的最佳状态转换序列中，其在i-1时刻的最佳状态，然后递归输出。

#include <stdio.h>
#define max(a,b) ((a)>(b)?(a):(b))
#define N 50
double s[N][N],t[N][N];//s表示状态转移矩阵，t[i][j]表示状态为i输出为j的概率
int out[N];//输出序列
double dp[N][N];
int path[N][N];//记录最佳路径
int T,n,m,p;//n表示状态数，m表示输出符号数量,p表示输出序列的数量
void print(int p,int i){
	if(p != 1)
		print(p-1,path[p][i]);
	printf("%d ",i);
}
int main(){
	freopen("a.txt","r",stdin);
	scanf("%d",&T);
	while(T--){
		int i,j,k,now;
		double res = 0;
		scanf("%d %d %d",&n,&m,&p);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=n;j++)
				scanf("%lf",&s[i][j]);
		for(i = 1;i<=n;i++)
			for(j = 1;j<=m;j++)
				scanf("%lf",&t[i][j]);
		for(i = 1;i<=p;i++)
			scanf("%d",&out[i]);
		
		for(i = 1;i<=n;i++){//初始化
			dp[1][i] = 1./n * t[i][out[1]];
			path[1][i] = 0;
		}
		for(i = 2;i<=p;i++)
			for(j = 1;j<=n;j++){
				now = 1;
				for(k = 2;k<=n;k++)//上一步在状态k，这一步转到状态j
					if(dp[i-1][k]*s[k][j] > dp[i-1][now]*s[now][j])
						now = k;
				dp[i][j] = dp[i-1][now] * s[now][j] * t[j][out[i]];
				path[i][j] = now;
			}
		for(i = 1;i<=n;i++)
			if(res < dp[p][i]){
				res = dp[p][i];
				now = i;//求出最佳路径的最后一个状态
			}
		printf("%lf\n",res);
		print(p,now);
	}
	return 0;
}