1.软件版本

matlab2013b

2.本算法理论知识

EM算法可以解决由不完全观测数据估计概率模型参数的问题。HMT模型中的EM算法思想是构造一个递推算法，使每一步估计的参数都使似然函数更大，使估计的参数逐渐接近真实值。其中，隐马尔科夫树的观测树的似然函数可以表示为：

隐马尔科夫树模型参数的E步算法

隐马尔科夫树模型参数的M步算法

3.部分源码

function [ESN,PSN,MUN,SIN]=emlht(w,ES,PS,MU,SI,zm)%
% Input :
% w : data PxP matrix (wavelet transform of an image)
% ES : state transition matrix (MxMxPxP)
% PS : state probability matrix (MxPxP)
% MU : mean matrix (MxPxP)
% SI : variance matrix (MxPxP)
% zm : type of density functions
%  zm = 1 : zero mean (do not update MU)
%  zm = 0 : nonzero mean (update MU)M=size(ES,1);
P=size(w,1);
level=log2(P);
BE=zeros(M,P,P);
BEP=zeros(M,P,P);
BER=zeros(M,P,P);AL=zeros(M,P,P);
P1=zeros(M,P,P);
P2=zeros(M,M,P,P);%UP step
wtmp = repmat(w,[1 1 M]);
wtmp = shiftdim(wtmp,2);
si=1; ei=2^(level-1); sj=2^(level-1)+1; ej=P;
gtmp = gauss(wtmp,MU,SI);
scale = repmat(mean(gtmp,1),[M 1 1]);
BE(:,si:ei,sj:ej) = gtmp(:,si:ei,sj:ej)./scale(:,si:ei,sj:ej);%clear MUtmp SItmp;for k=level:-1:2J=2^(k-1);J2=J*J; si = 1; ei = J; sj = J+1; ej = 2*J;EStmp = reshape(ES(:,:,si:ei,sj:ej),M,M*J2);if M == 2%%%%%% For M=2 the following is fasterBEtmp = zeros(M,M*J2);BEtmp(:,1:M:(M*J2))=reshape(BE(:,si:ei,sj:ej),M,J2);BEtmp(:,2:M:(M*J2))=BEtmp(:,1:M:(M*J2)); else% For general M (not equal to 2) use the following%BEtmp = zeros(M,M*4^(k-1));for m=1:MBEtmp(:,m:M:(M*4^(k-1)))=reshape(BE(:,si:ei,sj:ej,:),M,4^(k-1));end;end;BEtmp = reshape(EStmp.*BEtmp,[M M J J]);BEP(:,si:ei,sj:ej) = squeeze(sum(BEtmp,1));sni = 1; eni = J/2; snj = J/2+1; enj = J;%construct betachild matrix hereBCtmp = BEP(:,si:2:ei,sj:2:ej);BCtmp = BCtmp.*BEP(:,si+1:2:ei,sj:2:ej);BCtmp = BCtmp.*BEP(:,si:2:ei,sj+1:2:ej);BCtmp = BCtmp.*BEP(:,si+1:2:ei,sj+1:2:ej); scaletmp = repmat(mean(BCtmp,1),[M 1 1]);scale(:,sni:eni,snj:enj) = scale(:,sni:eni,snj:enj).*scaletmp;BE(:,sni:eni,snj:enj)=gtmp(:,sni:eni,snj:enj)./scale(:,sni:eni,snj:enj).*BCtmp;%construct BE(:,pai(i),paj(j),dindex) matrixBtmp=zeros(M,J,J);Btmp(:,1:2:J,1:2:J)=BE(:,sni:eni,snj:enj);Btmp(:,2:2:J,1:2:J)=BE(:,sni:eni,snj:enj);Btmp(:,1:2:J,2:2:J)=BE(:,sni:eni,snj:enj);Btmp(:,2:2:J,2:2:J)=BE(:,sni:eni,snj:enj);BER(:,si:ei,sj:ej)=Btmp./BEP(:,si:ei,sj:ej);
end;clear EStmp BEtmp BCtmp Btmp;%DOWN step%initializeAL(:,1,2) = PS(:,1,2);for k=2:levelJ = 2^(k-1); J2=J*J;si=1; ei=J; sj=J+1; ej=2*J;sni = 1; eni = J/2; snj = J/2+1; enj = J;Atmp=zeros(M,J,J);Atmp(:,1:2:J,1:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,2:2:J,1:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,1:2:J,2:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,2:2:J,2:2:J)=AL(:,sni:eni,snj:enj);Atmp = repmat(reshape(Atmp.*BER(:,si:ei,sj:ej),1,M*J2),[M 1]);EStmp = reshape(ES(:,:,si:ei,sj:ej),M,M*J2);ALtmp = reshape(EStmp.*Atmp,[M M J J]);AL(:,si:ei,sj:ej) = squeeze(sum(ALtmp,2));end;clear Atmp EStmp ALtmp;%compute probabilities
for k=2:levelJ=2^(k-1); J2=J*J;si=1; ei=J; sj=J+1; ej=2*J;sni = 1; eni = J/2; snj = J/2+1; enj = J;temp = repmat(sum(AL(:,si:ei,sj:ej).*BE(:,si:ei,sj:ej), 1),[M 1]);P1(:,si:ei,sj:ej) = AL(:,si:ei,sj:ej).*BE(:,si:ei,sj:ej)./temp;%compute P2if M == 2 % For M=2 the following may be fasterBEtmp = zeros(M,M*J2);BEtmp(:,1:M:(M*J2))=reshape(BE(:,si:ei,sj:ej),M,J2);BEtmp(:,2:M:(M*J2))=BEtmp(:,1:M:(M*J2));else% For general M (not equal to 2) use the followingBEtmp = zeros(M,M*J2);for m=1:MBEtmp(:,m:M:(M*J2))=reshape(BE(:,si:ei,sj:ej,:),M,J2);end;end;BEtmp = reshape(BEtmp,[M M J J]);EStmp = ES(:,:,si:ei,sj:ej);Atmp=zeros(M,J,J);Atmp(:,1:2:J,1:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,2:2:J,1:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,1:2:J,2:2:J)=AL(:,sni:eni,snj:enj);Atmp(:,2:2:J,2:2:J)=AL(:,sni:eni,snj:enj);Atmp = repmat(reshape(Atmp,1,M*J2),[M 1]);Atmp = reshape(Atmp,[M M J J]);BERtmp = repmat(reshape(BER(:,si:ei,sj:ej),1,M*J2),[M 1]);BERtmp = reshape(BERtmp,[M M J J]);temp = repmat(reshape(temp,1,M*J2),[M 1]);temp = reshape(temp, [M M J J]);P2(:,:,si:ei,sj:ej)=BEtmp.*EStmp.*Atmp.*BERtmp./temp;
end;
P1(:,1,2)=AL(:,1,2).*BE(:,1,2)./repmat(sum(AL(:,1,2).*BE(:,1,2),1),[M 1 1]);clear temp BEtmp EStmp Atmp BERtmp;%PSP=PS; ESP=ES; MUP=MU; SIP=SI;%M step
PS(:,1,2)=P1(:,1,2);for k=2:levelJ=2^(k-1); J2=J*J;si=1; ei=J; sj=J+1; ej=2*J;sni = 1; eni = J/2; snj = J/2+1; enj = J;pstmp = sum(sum(P1(:,si:ei,sj:ej),3),2)/J2;pstmp = pstmp.*(pstmp>1e-4)+1e-4*(pstmp<=1e-4);PS(:,si:ei,sj:ej) = repmat(pstmp,[1 J J]);if zm == 0% do not update MU if zero mean densitiesmutmp = sum(sum(wtmp(:,si:ei,sj:ej).*P1(:,si:ei,sj:ej),3),2)/J2;MU(:,si:ei,sj:ej) = repmat(mutmp,[1 J J])./PS(:,si:ei,sj:ej);end;sitmp = sum(sum((wtmp(:,si:ei,sj:ej)-MU(:,si:ei,sj:ej)).^2.*P1(:,si:ei,sj:ej),3),2)/J2;SI(:,si:ei,sj:ej) = repmat(sitmp,[1 J J])./PS(:,si:ei,sj:ej);estmp =sum(sum(P2(:,:,si:ei,sj:ej),4),3)/J2;ptmp = [PS(:,sni,snj)'; PS(:,sni,snj)'];ES(:,:,si:ei,sj:ej)= repmat(estmp,[1 1 J J])./repmat(ptmp,[1 1 J J]);
end; %kESN=ES; PSN=PS; MUN=MU; SIN=SI;

4.仿真结论

5.参考文献

[1]陈静, 范文兵, 甄姬娜. 基于小波变换及隐式马尔科夫树模型的图像信号去噪[J]. 现代电子技术, 2006.A09-06

【EM算法】小波域隐马尔科夫树模型参数的EM算法估计MATLAB仿真相关推荐

【深度】从朴素贝叶斯到维特比算法：详解隐马尔科夫模型
详解隐马尔科夫模型作者:David S. Batista 选自:机器之心本文首先简要介绍朴素贝叶斯,再将其扩展到隐马尔科夫模型.我们不仅会讨论隐马尔科夫模型的基本原理,同时还从朴素贝叶斯的角度讨论 ...
隐马尔科夫(HMM)模型
隐马尔科夫(Hidden Markov model)模型是一类基于概率统计的模型,是一种结构最简单的动态贝叶斯网,是一种重要的有向图模型.自上世纪80年代发展起来,在时序数据建模,例如:语音识别.文字 ...
隐马尔科夫模型(前向后向算法、鲍姆-韦尔奇算法、维特比算法)
隐马尔科夫模型(前向后向算法.鲍姆-韦尔奇算法.维特比算法) 概率图模型是一类用图来表达变量相关关系的概率模型.它以图为表示工具,最常见的是用一个结点表示一个或一组随机变量,结点之间的变表是变量间的概 ...
机器学习算法总结(七)——隐马尔科夫模型(前向后向算法、鲍姆-韦尔奇算法、维特比算法)...
概率图模型是一类用图来表达变量相关关系的概率模型.它以图为表示工具,最常见的是用一个结点表示一个或一组随机变量,结点之间的变表是变量间的概率相关关系.根据边的性质不同,可以将概率图模型分为两类:一类是 ...
【深度剖析HMM（附Python代码）】2.隐马尔科夫链HMM的EM训练过程
隐马尔科夫链HMM的参数θ的EM训练过程现在回到前一节最后提出的参数θ的最大似然函数上来,先对其做个对数变换,做对数变换是考虑到序列X的概率计算公式中包含了连乘,为了方便计算同时避免序列X的概率过小 ...
隐马尔科夫模型(HMMs)之五：维特比算法及前向后向算法
维特比算法(Viterbi Algorithm) 找到可能性最大的隐藏序列通常我们都有一个特定的HMM,然后根据一个可观察序列去找到最可能生成这个可观察序列的隐藏序列. 1.穷举搜索我们可以在下图 ...
视频教程-隐马尔科夫算法：中文分词神器-深度学习
隐马尔科夫算法:中文分词神器在中国知网从事自然语言处理和知识图谱的开发,并负责带领团队完成项目,对深度学习和机器学习算法有深入研究. 吕强 ¥49.00 立即订阅扫码下载「CSDN程序员学院APP ...
学习笔记四：HMM（隐马尔科夫模型）
第四章 HMM 1.隐马尔科夫模型的基本概念引入例子:天气和冰激凌问题想象你是2799年的气候学家,正在研究全球变暖的历史你无法找到任何关于北京2020年夏天的天气记录但你找到了小米的日记本 ...
隐马尔科夫模型(Hidden Markov Models) 系列之五
隐马尔科夫模型(Hidden Markov Models) 系列之五介绍(introduction) 生成模式(Generating Patterns) 隐含模式(Hidden Patterns) ...
一、隐马尔科夫模型HMM
隐马尔科夫模型HMM(一)HMM模型基础隐马尔科夫模型(Hidden Markov Model,以下简称HMM)是比较经典的机器学习模型了,它在语言识别,自然语言处理,模式识别等领域得到广泛的应用. ...

【EM算法】小波域隐马尔科夫树模型参数的EM算法估计MATLAB仿真