隐马尔科夫（HMM）的Matlab实现

说明：

1. 本文实现了PRML一书第13章的隐马尔科夫（HMM）算法，并与K-means聚类、GMM模型聚类进行了对比。当然，HMM的用处远不止是聚类；

2. 非职业码农，代码质量不高，变量命名也不规范，凑合着看吧，不好意思。

隐马尔科夫模型（HMM）是我一直想弄清楚的一种模型。可能是上学时随机过程只考了60分，心有不甘。

模型和算法本身不赘述了，详见PRML第13章。直接看结果，把隐变量z作为类别编号，HMM可用于聚类。这里可以看到，由于HMM利用了序列跳转信息（也就是马尔科夫特性），可以获得比混合高斯模型（GMM，http://blog.csdn.net/foreseerwang/article/details/75222522）更为准确的结果。

Matlab输出：

****************************************************************

K-means聚类计算，结果将作为HMM初始值......

K-means聚类算法结束

****************************************************************

GMM聚类计算中......

GMM聚类算法结束

****************************************************************

HMM EM算法迭代计算中......

loop 1, A error: 7.19e-01

loop 2, A error: 4.20e-02

loop 3, A error: 1.36e-02

loop 4, A error: 5.59e-03

loop 5, A error: 2.96e-03

loop 6, A error: 1.61e-03

loop 7, A error: 8.77e-04

loop 8, A error: 4.87e-04

loop 9, A error: 2.81e-04

loop 10, A error: 1.70e-04

loop 11, A error: 1.08e-04

loop 12, A error: 7.28e-05

HMM EM算法结束

****************************************************************

真实质心为：[1.000,2.000], [-4.000, 2.000], [7.000, 1.000]

Kmeans聚类质心为：[0.647, 1.799],[-4.714, 2.394], [6.713, 0.870]

GMM聚类质心为：[0.977, 2.009],[-4.179, 2.010], [6.945, 0.973]

HMM拟合质心为：[1.011, 1.994],[-4.091, 2.005], [6.984, 1.001]

****************************************************************

Kmeans聚类错误率：10.517%

GMM 聚类错误率： 3.360%

HMM 拟合错误率： 1.673%

****************************************************************

真实转移概率矩阵A：

0.8000 0.1000 0.1000

0.2000 0.7000 0.1000

0 0.5000 0.5000

HMM拟合转移概率矩阵A：

0.7978 0.0981 0.1041

0.2038 0.6963 0.0999

0.0006 0.4987 0.5007

****************************************************************

聚类结果的图示：

代码：

%% by foreseer wang
% web: http://blog.csdn.net/foreseerwang
% QQ: 50834
%% 新买了个机械键盘，茶轴，码字很爽clear all;
close all;
rng('default');
fprintf('****************************************************************\n');
fprintf('****************************************************************\n');%% 初始化并生成测试数据
dim=[2,30000];   % k*N，k为数据维度，N为数据数目
Nclst=3;         % 簇的数量，也就是隐变量z可能取值的个数
k=dim(1);
N=dim(2);
Aerr=1e-4;       % 用于迭代终止A0real = [0.2, 0.3, 0.5];                   % z1点分布概率
Areal = [0.8, 0.1, 0.1;...0.2, 0.7, 0.1;...0.0, 0.5, 0.5];                    % z序列转移概率矩阵mureal = [1 2; -4 2; 7 1]';                 % p(x|z)的均值
sigreal=zeros(k,k,Nclst);                   % p(x|z)的方差
sigreal(:,:,1)=[2 -1.5; -1.5 2];
sigreal(:,:,2)=[5 -2.; -2. 3];
sigreal(:,:,3)=[1 0.1; 0.1 2];zreal=zeros(1,dim(2));                      % 隐变量z
zreal(1)=randsrc(1,1,[1:Nclst;A0real]);
x=zeros(dim);                               % 可观测到的变量x
x(:,1)=funGaussSample(mureal(:,zreal(1)),squeeze(sigreal(:,:,zreal(1))),1);
for ii=2:dim(2),zreal(ii)=randsrc(1,1,[1:Nclst;Areal(zreal(ii-1),:)]);x(:,ii)=funGaussSample(mureal(:,zreal(ii)),squeeze(sigreal(:,:,zreal(ii))),1);
end;%% 利用EM算法求解HMM问题
% 只有序列x是已知量
% 通过x求解mu（p(x|z)的均值）、sig（p(x|z)的协方差）、
%   A0（z1的分布概率）、A（z序列的状态跳转概率矩阵）% 其中涉及的中间变量包括（详见PRML第13章）：
%   alpha(zn) = p(x1,..., xn, zn)（未用到）
%   alp_caret(zn) = p(zn|x1,..., xn) = alpha(zn)/p(x1, ..., xn)
%   beta(zn) = p(x(n+1), ..., xN|zn)（未用到）
%   bet_caret(zn) = p(x(n+1), ..., xN|zn)/p(x(n+1), ..., xN|x1,...xn)
%               = beta(zn)/p(x(n+1), ..., xN|x1,...xn)
%   c(n) = p(xn|x1, ..., x(n-1))%   gamma(zn) = p(zn|X)=alpha(zn)*beta(zn)/p(X) = alp_caret(zn)*bet_caret(zn)
%   xi(z(n-1), zn) = alpha(z(n-1))*p(xn|zn)*p(zn|z(n-1))*beta(zn)/p(X)
%                  = alp_caret(z(n-1))*p(xn|zn)*p(zn|z(n-1))*bet_caret(zn)/c(n)% 算法步骤：
% 1. 变量随机初始化，包括：mu（k*Nclst）、sig（k*k*Nclst）、A0（1*Nclst）、
%    A（k*k）、z(Nclst*N)
% 2. E步和M步轮换迭代：
%    2.1 使用forward-backward算法alp_caret和bet_caret，进而得出gama和xi
%    2.2 M步：重新计算mu、sig、A0、A
% 3. 使用Viterbi算法求得最可能的z序列% 1. 初始化 ---------------------------------------------------------------
mu = zeros(k,Nclst);
sig = zeros(k,k,Nclst);
z = zeros(Nclst,N);% K-means聚类，结果作为HMM的初始值
fprintf('K-means聚类计算，结果将作为HMM初始值......\n');
k_idx=kmeans(x',Nclst);         % k_idx即为聚类后的z
tmp_k_idx1=zeros(1,Nclst);
mu_kmeans=zeros(k, Nclst);
for ii=1:Nclst,idx=(k_idx==ii);tmp_mu=mean(x(:,idx),2);    % 旋转操作，确保类别编号与真实数据一致dist=sum((repmat(tmp_mu,1,Nclst)-mureal).^2);[dummy,tmp_idx]=min(dist);mu_kmeans(:,tmp_idx)=tmp_mu;mu(:,tmp_idx)=tmp_mu;sig(:,:,tmp_idx)=cov(x(:,idx)');tmp_k_idx1(ii)=tmp_idx;
end;tmp_k_idx2=k_idx;
for ii=1:Nclst,k_idx(tmp_k_idx2==ii)=tmp_k_idx1(ii);
end;
errrate_Clst=sum(abs(k_idx'-zreal)>0.5)*1.0/N; % kmeans聚类错误率
fprintf('K-means聚类算法结束\n');
fprintf('****************************************************************\n');% GMM聚类 -----------------------------------------------------------------
% 此处仅为对比，相对独立。删除此处横线间的代码及最后的相关打印代码，对结果无影响
fprintf('GMM聚类计算中......\n');
k_idx_GMM=funGMM(x',Nclst);     % k_idx_GMM为GMM聚类后的z
mu_GMM = zeros(k,Nclst);
sig_GMM = zeros(k,k,Nclst);
for ii=1:Nclst,idx=(k_idx_GMM==ii);tmp_mu=mean(x(:,idx),2);    % 旋转dist=sum((repmat(tmp_mu,1,Nclst)-mureal).^2);[dummy,tmp_idx]=min(dist);mu_GMM(:,tmp_idx)=tmp_mu;sig_GMM(:,:,tmp_idx)=cov(x(:,idx)');tmp_k_idx1(ii)=tmp_idx;
end;tmp_k_idx2=k_idx_GMM;
for ii=1:Nclst,k_idx_GMM(tmp_k_idx2==ii)=tmp_k_idx1(ii);
end;errrate_GMM=sum(abs(k_idx_GMM'-zreal)>0.5)*1.0/N; % GMM聚类错误率
fprintf('GMM聚类算法结束\n');
fprintf('****************************************************************\n');
% GMM聚类结束 --------------------------------------------------------------% K-means聚类结果赋给隐变量z
for ii=1:N,z(k_idx(ii),ii)=1;
end;% 随机初始化A0和A
tmp = rand(1,Nclst);
A0 = tmp/sum(tmp);
tmp = rand(Nclst, Nclst);
A = tmp./repmat(sum(tmp,2),1,Nclst);% 2. EM迭代
alp_caret = zeros(Nclst, N);
c=zeros(1,N);
Pzx=zeros(Nclst,N);         % p(x|z)bet_caret = zeros(Nclst, N);
bet_caret(:,N) = ones(Nclst, 1);gama = zeros(Nclst,N);
xi = zeros(Nclst,Nclst,N-1);
Aold=ones(size(A))/Nclst;
fprintf('HMM EM算法迭代计算中......\n');for ii=1:100,for jj=1:Nclst,Pzx(jj,:)=mvnpdf(x',mu(:,jj)',squeeze(sig(:,:,jj)))';end;% 2.1 E step ----------------------------------------------------------tmp1=A0'.*Pzx(:,1);    % 实际上是alpha(z1)c(1)=sum(tmp1);alp_caret(:,1)=tmp1/c(1);% 迭代计算PRML式13.59for kk=2:N,tmp1=alp_caret(:,kk-1)'*A;tmp2=tmp1'.*Pzx(:,kk);c(kk)=sum(tmp2);alp_caret(:,kk)=tmp2/c(kk);end;% end of forward. alpha_caret got% 迭代计算PRML式13.62for kk=N:-1:2,tmp2=A*(Pzx(:,kk).*bet_caret(:,kk));bet_caret(:,kk-1)=tmp2/c(kk);end;% end of backward. beta_caret gotgama = alp_caret.*bet_caret;for jj=1:N-1,xi(:,:,jj)=alp_caret(:,jj)/c(jj+1)*...(Pzx(:,jj+1).*bet_caret(:,jj+1))'.*A;end;% gama and xi got% 2.2 M step ----------------------------------------------------------A0=gama(:,1)'/sum(gama(:,1));A=sum(xi,3)./repmat(sum(sum(xi,3),2),1,Nclst);mu=x*gama'./repmat(transpose(sum(gama,2)),k,1);for kk=1:Nclst,sig(:,:,kk)=(repmat(gama(kk,:),k,1).*(x-repmat(mu(:,kk),1,N)))*...(x-repmat(mu(:,kk),1,N))'/sum(gama(kk,:));end;% M step completediter_err=sum(sum(abs(A-Aold)))/Nclst;Aold=A;fprintf('loop %2d, A error: %4.2e \n', ii, iter_err);if iter_err<Aerr,break;end;
end;
fprintf('HMM EM算法结束\n');
fprintf('****************************************************************\n');for jj=1:Nclst,Pzx(jj,:)=mvnpdf(x',mu(:,jj)',squeeze(sig(:,:,jj)))';
end;% 3. Viterbi算法确定最可能的序列 -------------------------------------------
w = zeros(Nclst,N);
psi = zeros(Nclst, N-1);
zpred=zeros(1,N);w(:,1)=log(A0')+log(Pzx(:,1));
% 迭代计算PRML式13.70
for ii=1:N-1,tmp_1=log(A)+repmat(w(:,ii),1,Nclst);[tmp_max, tmp_idx]=max(tmp_1, [], 1);w(:,ii+1)=log(Pzx(:,ii+1))+tmp_max';psi(:,ii)=tmp_idx';
end;[dummy, tmp_idx]=max(w(:,N));
% 反向计算PRML式13.71
zpred(N)=tmp_idx;
for ii=N-1:-1:1,zpred(ii)=psi(zpred(ii+1),ii);
end;errrate_HMM=sum(abs(zpred-zreal)>0.5)*1.0/N;
for ii=1:Nclst,mu(:,ii)=mean(x(:,zpred==ii),2);
end;%% 结果打印
% 打印聚类质心
fprintf('真实质心为      ：[%5.3f, %5.3f], [%5.3f, %5.3f], [%5.3f, %5.3f]\n',...mureal(1,1),mureal(2,1),mureal(1,2),mureal(2,2),mureal(1,3),mureal(2,3));
fprintf('Kmeans聚类质心为：[%5.3f, %5.3f], [%5.3f, %5.3f], [%5.3f, %5.3f]\n',...mu_kmeans(1,1),mu_kmeans(2,1),mu_kmeans(1,2),...mu_kmeans(2,2),mu_kmeans(1,3),mu_kmeans(2,3));
fprintf('GMM聚类质心为   ：[%5.3f, %5.3f], [%5.3f, %5.3f], [%5.3f, %5.3f]\n',...mu_GMM(1,1),mu_GMM(2,1),mu_GMM(1,2),mu_GMM(2,2),mu_GMM(1,3),mu_GMM(2,3));
fprintf('HMM拟合质心为   ：[%5.3f, %5.3f], [%5.3f, %5.3f], [%5.3f, %5.3f]\n',...mu(1,1),mu(2,1),mu(1,2),mu(2,2),mu(1,3),mu(2,3));% 打印散点图，直观观看效果
% 可以看到，HMM算法中用到了序列信息，结果更为准确
figure(1);
subplot(2,2,1); hold on;
scatter(x(1,zreal==1),x(2,zreal==1),'ko');
scatter(x(1,zreal==2),x(2,zreal==2),'go');
scatter(x(1,zreal==3),x(2,zreal==3),'bo');
xlabel('x1');
ylabel('x2');
axis([-15,15,-5,10]);
title('1. 原始数据', 'FontSize', 20);subplot(2,2,2); hold on;
scatter(x(1,k_idx==1),x(2,k_idx==1),'ko');
scatter(x(1,k_idx==2),x(2,k_idx==2),'go');
scatter(x(1,k_idx==3),x(2,k_idx==3),'bo');
xlabel('x1');
ylabel('x2');
axis([-15,15,-5,10]);
title('2. K-means聚类结果', 'FontSize', 20);subplot(2,2,3); hold on;
scatter(x(1,k_idx_GMM==1),x(2,k_idx_GMM==1),'ko');
scatter(x(1,k_idx_GMM==2),x(2,k_idx_GMM==2),'go');
scatter(x(1,k_idx_GMM==3),x(2,k_idx_GMM==3),'bo');
xlabel('x1');
ylabel('x2');
axis([-15,15,-5,10]);
title('3. GMM聚类结果', 'FontSize', 20);subplot(2,2,4); hold on;
scatter(x(1,zpred==1),x(2,zpred==1),'ko');
scatter(x(1,zpred==2),x(2,zpred==2),'go');
scatter(x(1,zpred==3),x(2,zpred==3),'bo');
xlabel('x1');
ylabel('x2');
axis([-15,15,-5,10]);
title('4. HMM结果', 'FontSize', 20);% 打印聚类错误率，其实就是隐变量z的错误率
fprintf('****************************************************************\n');
fprintf('Kmeans聚类错误率：%6.3f%%\n', errrate_Clst*100);
fprintf('GMM   聚类错误率：%6.3f%%\n', errrate_GMM*100);
fprintf('HMM   拟合错误率：%6.3f%%\n', errrate_HMM*100);
fprintf('****************************************************************\n');% 打印z序列转移概率矩阵
fprintf('真实转移概率矩阵A：\n');
disp(Areal);
fprintf('HMM拟合转移概率矩阵A：\n');
disp(A);
fprintf('****************************************************************\n');

子程序1（就一句话，好像有点多此一举）：

function z = funGaussSample( mu, sigma, dim )
%GAUSSAMPLE Summary of this function goes here
%   Detailed explanation goes here%R = chol(sigma);
%z = repmat(mu,dim(1),1) + randn(dim)*R;
z = mvnrnd(mu, sigma, dim(1));end

子程序2（混合高斯模型聚类算法）：

function clst_idx = funGMM( z, Nclst )
%% by foreseer wang
% web: http://blog.csdn.net/foreseerwang
% QQ: 50834[len, k]=size(z);k_idx=kmeans(z,Nclst);     % 用K-means聚类结果做GMM聚类的初始值
mu=zeros(Nclst,k);
sigma=zeros(k,k,Nclst);
for ii=1:Nclst,mu(ii,:)=mean(z(k_idx==ii,:));sigma(:,:,ii)=cov(z(k_idx==ii,:));
end;
Pw=ones(Nclst,1)*1.0/Nclst;px=zeros(len,Nclst);
r=zeros(len,Nclst);
for jj=1:1000, % 简单起见，直接循环，不做结束判断for ii=1:Nclst,px(:,ii)=mvnpdf(z,mu(ii,:),squeeze(sigma(:,:,ii)));end;% E steptemp=px.*repmat(Pw',len,1);r=temp./repmat(sum(temp,2),1,Nclst);% M steprk=sum(r);pw=rk/len;mu=r'*z./repmat(rk',1,k);for ii=1:Nclstsigma(:,:,ii)=z'*(repmat(r(:,ii),1,k).*z)/rk(ii)-mu(ii,:)'*mu(ii,:);end;
end;% display
[dummy,clst_idx]=max(px,[],2);end

隐马尔科夫（HMM）的Matlab实现相关推荐

隐马尔科夫(HMM)模型
隐马尔科夫(Hidden Markov model)模型是一类基于概率统计的模型,是一种结构最简单的动态贝叶斯网,是一种重要的有向图模型.自上世纪80年代发展起来,在时序数据建模,例如:语音识别.文字 ...
【深度剖析HMM（附Python代码）】1.前言及隐马尔科夫链HMM的背景
1. 前言隐马尔科夫HMM模型是一类重要的机器学习方法,其主要用于序列数据的分析,广泛应用于语音识别.文本翻译.序列预测.中文分词等多个领域.虽然近年来,由于RNN等深度学习方法的发展,HMM模型逐 ...
python做空矩阵_【手把手教你】Python实现基于隐马尔科夫的多空策略
前言我们通常使用股市的一手数据来创建一个策略模型,预测下一时刻价格的多少.走势的判断或其他. 今天,我们想结合多样的市场条件(波动性,交易量,价格变化等等)和结合隐马尔科夫(HMM)来构建我们的交易 ...
HMM隐马尔科夫时间序列预测 Markov马尔科夫时间序列预测（Matlab）
HMM隐马尔科夫时间序列预测 Markov马尔科夫时间序列预测(Matlab) 1.所有程序经过验证,保证可以运行 2.程序包括源码(主程序一个,子函数两个)和数据集: 3.程序适用于单变量时间序列预 ...
m基于隐马尔科夫模型(HMM)的手机用户行为预测(MMUB)算法matlab仿真
目录 1.算法描述 2.仿真效果预览 3.MATLAB核心程序 4.完整MATLAB 1.算法描述隐马尔可夫模型(Hidden Markov Model,HMM)是一种统计模型,广泛应用在语音识别, ...
隐马尔科夫模型matlab工具箱说明
转自 http://blog.csdn.net/whiteinblue/article/details/40625291 隐马尔可夫模型(HiddenMarkov Model,HMM)是统计模型,它用 ...
隐马尔科夫模型HMM（一）HMM模型
2019独角兽企业重金招聘Python工程师标准>>> 隐马尔科夫模型(Hidden Markov Model,以下简称HMM)是比较经典的机器学习模型了,它在语言识别,自然语言处理 ...
机器学习知识点(二十五)Java实现隐马尔科夫模型HMM之jahmm库
1.隐马尔可夫模型HMM的应用场景,关乎于序列和状态变化的都可以. 发现java有可库,专为开发HMM,可惜只能在CSDN上有得下载. 2.jahmm是java开发隐马尔科夫模型的一个j ...
隐马尔科夫模型HMM自学（3）
Viterbi Algorithm 本来想明天再把后面的部分写好,可是睡觉今天是节日呢?一时情不自禁就有打开电脑.......... 找到可能性最大的隐含状态序列崔晓源翻译多数情况下,我们都希望 ...

隐马尔科夫（HMM）的Matlab实现

隐马尔科夫（HMM）的Matlab实现相关推荐

最新文章

热门文章