HMM在3D轨迹识别中的应用matlab——适合小白
近日在用HMM做3D轨迹识别，看了许多博客与文章，在此做一个总结，希望对刚刚接触HMM的人有所帮助。
关于HMM原理方面，我爱自然语言处理博客有很细致的介绍，在这里就不妄自菲薄地介绍HMM了。
使用HMM解决问题，我们需要清楚地知道“什么是马尔科夫过程”（即我的问题是否符合马尔科夫过程的条件），隐马尔科夫模型解决的三个问题（评估，解码，学习）的用途（三个问题的介绍请移步传送门），三个问题的输入输出是什么，我们需要解决的问题是属于哪一类，才能按照需求去寻找解决问题的方法。比如以我需要解决的轨迹识别为例，我的训练数据是n个轨迹样本，每个样本是一条3D轨迹时间序列，每个样本的时间帧是不同的，这n个轨迹样本共分为五类（有标签的），我希望输入一条轨迹样本 ，输出这个轨迹的类别。在HMM中我的问题要拆分为两个，首先它是个“学习”问题，我需要分别训练出5个行为的HMM的“状态转移矩阵A”和“混淆矩阵B”；然后它是个“评估问题”，我需要将待识别的序列分别放进这5个HMM中，得到5个loglik,这条序列属于这个HMM的概率（取对数后的值），这5个中哪个最大，这条序列就属于哪个行为。
“Talk is cheap，Show me the code.”用一个简单的matlab代码来举例“训练”和“识别”的过程，用的是matlab自带的hmm代码。有关matlab自带的hmm代码解释，请戳matlab官网和博客。不知道自己说的对不对，根据自己的调试来看，hmmtrain函数好像只接受一维的训练样本，即Trseq中的时间序列【3 3 3 4 5 4 5 5 6 10 ……10 11 11 12 12 1 1】只能是一维的，假如放多维的数据会在运行过程中出现错误。
%简单的离散HMM“学习”示例
statenum=5;%状态数 可以变，状态数不一样 识别性能会些许改变
trans=rand(statenum);
sumTrans=sum(trans,2);
Temp=inv(diag(sumTrans));
trans=(trans'*Temp)';%状态转移矩阵
observationnum=23;%观察状态数
emis = ones(statenum,observationnum)*(1/observationnum);%混合矩阵
%Trseq为某一类型轨迹训练样本集，每一条轨迹的点数可能不一样。
%故Trseq为一个cell结构。
%比如说：
Trseq={[3 3 3 4 5 4 5 5 6 6 7 7 8 8 8 9 10 10 11 12 1 2 22 3 3 3 5 8 8 9 10 10 10 11 11 12 12 1 1]
[4 5 5 5 6 6 7 8 8 9 9 10 11 12 1 2 2 3 3 3 3 7 9 9 10 10 11 12 12 1 1 2 3]
[4 5 6 6 7 7 8 8 9 10 10 11 1 2 2 3 3 3 3 1 9 9 10 10 11 12 12 1 1 3]
[4 5 5 6 6 7 7 8 9 9 10 11 12 2 3 3 3 3 3 10 9 9 10 10 11 12 12 1 2]
[4 4 4 5 5 6 6 7 8 8 9 9 10 10 12 1 2 3 3 3 3 3 4 10 9 9 10 10 10 11 12 12 1 1]
[4 4 5 5 6 7 7 8 8 9 9 10 11 12 1 2 3 3 3 3 3 4 1 8 9 9 10 10 11 11 12 1 1 2]
[4 4 4 5 5 6 6 7 8 8 9 9 10 10 11 1 2 2 3 3 3 3 3 1 8 9 10 10 11 11 12 12 1 12]
[3 4 4 4 5 5 6 7 7 8 9 9 10 10 11 12 1 2 3 3 3 3 3 4 9 9 10 10 10 11 12 12 1 12]
[2 4 4 4 4 4 5 6 7 7 8 9 9 10 10 11 1 2 2 3 3 3 3 4 8 9 10 10 11 12 12 1 1 2]};
%estTrain,estEmis分别指 训练得到的状态转移矩阵，混合转移矩阵（状态到观察值转移矩阵）
[estTrain,estEmis] = hmmtrain(Trseq,trans,emis);%训练
%Ttest为一条即将要识别的类型，或者测试的轨迹编码，比如说
Ttest=[3 4 4 4 4 4 5 6 7 8 8 9 9 10 10 11 12 1 2 3 3 3 12 8 10 10 10 11 11 12 1 2];
%Ttest =[12 11 12 10 3 3 12 8 10 10 10 11 11 12 1 2 3 4 4 4 4 4 5 6 7 8 8 9 9 10 10 11];
%通过hmmdecode计算之后得到logPro。
[~,logPro]= hmmdecode(Ttest,estTrain,estEmis);%测试
%不同的训练样本集，训练得出不同的HMM
%给定的一条轨迹编码，带入不同的HMM，得到相应的logPro，选择tlog最大的值对应的HMM类型即为该类型的轨迹模式。
因为我的输入样本是3D轨迹数据，因此不能使用自带的hmm，参考几种不同程序语言的HMM版本，我使用的是博客里介绍的matlab toolbox，也希望使用过其他语言版本HMM的同学能分享一下使用经验哈~
这一版matlab中虽然m文件很多，看着有些吓人，不过其实逻辑很清晰，主页上有介绍它的文件内容：
This toolbox supports inference and learning for HMMs with discrete outputs (dhmm's), Gaussian outputs (ghmm's), or mixtures of Gaussians output (mhmm's). The Gaussians can be full, diagonal, or spherical (isotropic). It also supports discrete inputs, as in a POMDP. The inference routines support filtering, smoothing, and fixed-lag smoothing.
我用的是ghmm即GMM-HMM，有篇博客（传送门）对它进行了简洁的概括性介绍，toolbox的使用说明中也有相关解释，其实耐心下来看懂英文解释就足够理解是什么意思了（奈何英文阅读速度还需磨砺……）。这个版本的代码刚刚放入matlab中可能会报错,workspacefunc 329,解决方案在这里，就是将与matlab自带的toolbox重名的函数改过名就ok了，记得是KPMtools的”strsplit“
GMM-HMM的代码是HMM文件夹下的带“mhmm”关键字的m文件，具体用法参考mhmm_em_demo.m
%mhmm_em_demo.m
if 1 O = 4; T = 10; nex = 50; M = 2; Q = 3; else O = 8; %Number of coefficients in a vector 向量维度，在我的问题中由于是3D轨迹，所以是3 T = 420; %Number of vectors in a sequence T序列长度 nex = 1; %Number of sequences nex样本数量 M = 1; %Number of mixtures Q = 6; %Number of states end cov_type = 'full'; data = randn(O,T,nex); % initial guess of parameters prior0 = normalise(rand(Q,1)); transmat0 = mk_stochastic(rand(Q,Q)); if 0 Sigma0 = repmat(eye(O), [1 1 Q M]); % Initialize each mean to a random data point indices = randperm(T*nex); mu0 = reshape(data(:,indices(1:(Q*M))), [O Q M]); mixmat0 = mk_stochastic(rand(Q,M)); else [mu0, Sigma0] = mixgauss_init(Q*M, data, cov_type); mu0 = reshape(mu0, [O Q M]); Sigma0 = reshape(Sigma0, [O O Q M]); mixmat0 = mk_stochastic(rand(Q,M)); end [LL, prior1, transmat1, mu1, Sigma1, mixmat1] = ... mhmm_em(data, prior0, transmat0, mu0, Sigma0, mixmat0, 'max_iter', 5); loglik = mhmm_logprob(data, prior1, transmat1, mu1, Sigma1, mixmat1);
最后，祝大家调代码顺利~
