EM算法matlab程序

X=zeros(600,2);X(1:200,:)=normrnd(0,1,200,2);X(201:400,:)=normrnd(0,2,200,2);X(401:600,:)=normrnd(0,3,200,2);[W,M,V,L]=EM_GM(X,3,[],[],1,[])下面是程序源码：查看源代码打印帮助function[W,M,V,L]=EM_GM(X,k,ltol,maxiter,pflag,Init)%[W,M,V,L]=EM_GM(X,k,ltol,maxiter,pflag,Init)%%EMalgorithmforkmultidimensionalGaussianmixtureestimation%%Inputs:%X(n,d)-inputdata,n=numberofobservations,d=dimensionofvariable%k-maximumnumberofGaussiancomponentsallowed%ltol-percentageoftheloglikelihooddifferencebetween2iterations([]fornone)%maxiter-maximumnumberofiterationallowed([]fornone)%pflag-1forplottingGMfor1Dor2Dcasesonly,0otherwise([]fornone)%Init-structureofinitialW,M,V:Init.W,Init.M,Init.V([]fornone)%%Ouputs:%W(1,k)-estimatedweightsofGM%M(d,k)-estimatedmeanvectorsofGM%V(d,d,k)-estimatedcovariancematricesofGM%L-loglikelihoodofestimates%%Writtenby%PatrickP.C.Tsui,%PAMIresearchgroup%DepartmentofElectricalandComputerEngineering%UniversityofWaterloo,%March,2006%%%%%Validateinputs%%%%ifnargin=1,disp('EM_GMmusthaveatleast2inputs:X,k!/n')returnelseifnargin==2,ltol=0.1;maxiter=1000;pflag=0;Init=[];err_X=Verify_X(X);err_k=Verify_k(k);iferr_X|err_k,return;endelseifnargin==3,maxiter=1000;pflag=0;Init=[];err_X=Verify_X(X);err_k=Verify_k(k);[ltol,err_ltol]=Verify_ltol(ltol);iferr_X|err_k|err_ltol,return;endelseifnargin==4,pflag=0;Init=[];err_X=Verify_X(X);err_k=Verify_k(k);[ltol,err_ltol]=Verify_ltol(ltol);[maxiter,err_maxiter]=Verify_maxiter(maxiter);iferr_X|err_k|err_ltol|err_maxiter,return;endelseifnargin==5,Init=[];err_X=Verify_X(X);err_k=Verify_k(k);[ltol,err_ltol]=Verify_ltol(ltol);[maxiter,err_maxiter]=Verify_maxiter(maxiter);[pflag,err_pflag]=Verify_pflag(pflag);iferr_X|err_k|err_ltol|err_maxiter|err_pflag,return;endelseifnargin==6,err_X=Verify_X(X);err_k=Verify_k(k);[ltol,err_ltol]=Verify_ltol(ltol);[maxiter,err_maxiter]=Verify_maxiter(maxiter);[pflag,err_pflag]=Verify_pflag(pflag);[Init,err_Init]=Verify_Init(Init);iferr_X|err_k|err_ltol|err_maxiter|err_pflag|err_Init,return;endelsedisp('EM_GMmusthave2to6inputs!');returnend%%%%InitializeW,M,V,L%%%%t=cputime;ifisempty(Init),[W,M,V]=Init_EM(X,k);L=0;elseW=Init.W;M=Init.M;V=Init.V;endLn=Likelihood(X,k,W,M,V);%InitializeloglikelihoodLo=2*Ln;%%%%EMalgorithm%%%%niter=0;while(abs(100*(Ln-Lo)/Lo)ltol)&(niter=maxiter),E=Expectation(X,k,W,M,V);%E-step[W,M,V]=Maximization(X,k,E);%M-stepLo=Ln;Ln=Likelihood(X,k,W,M,V);niter=niter+1;endL=Ln;%%%%Plot1Dor2D%%%%ifpflag==1,[n,d]=size(X);ifd2,disp('Canonlyplot1or2dimensionalapplications!/n');elsePlot_GM(X,k,W,M,V);endelapsed_time=sprintf('CPUtimeusedforEM_GM:%5.2fs',cputime-t);disp(elapsed_time);disp(sprintf('Numberofiterations:%d',niter-1));end%%%%%%%%%%%%%%%%%%%%%%%%%%EndofEM_GM%%%%%%%%%%%%%%%%%%%%%%%%%%functionE=Expectation(X,k,W,M,V)[n,d]=size(X);a=(2*pi)^(0.5*d);S=zeros(1,k);iV=zeros(d,d,k);forj=1:k,ifV(:,:,j)==zeros(d,d),V(:,:,j)=ones(d,d)*eps;endS(j)=sqrt(det(V(:,:,j)));iV(:,:,j)=inv(V(:,:,j));endE=zeros(n,k);fori=1:n,forj=1:k,dXM=X(i,:)'-M(:,j);pl=exp(-0.5*dXM'*iV(:,:,j)*dXM)/(a*S(j));E(i,j)=W(j)*pl;endE(i,:)=E(i,:)/sum(E(i,:));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofExpectation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[W,M,V]=Maximization(X,k,E)[n,d]=size(X);W=zeros(1,k);M=zeros(d,k);V=zeros(d,d,k);fori=1:k,%Computeweightsforj=1:n,W(i)=W(i)+E(j,i);M(:,i)=M(:,i)+E(j,i)*X(j,:)';endM(:,i)=M(:,i)/W(i);endfori=1:k,forj=1:n,dXM=X(j,:)'-M(:,i);V(:,:,i)=V(:,:,i)+E(j,i)*dXM*dXM';endV(:,:,i)=V(:,:,i)/W(i);endW=W/n;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofMaximization%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionL=Likelihood(X,k,W,M,V)%ComputeLbasedonK.V.Mardia,MultivariateAnalysis,AcademicPress,1979,PP.96-97%toenchancecomputationalspeed[n,d]=size(X);U=mean(X)';S=cov(X);L=0;fori=1:k,iV=inv(V(:,:,i));L=L+W(i)*(-0.5*n*log(det(2*pi*V(:,:,i)))...-0.5*(n-1)*(trace(iV*S)+(U-M(:,i))'*iV*(U-M(:,i))));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofLikelihood%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionerr_X=Verify_X(X)err_X=1;[n,d]=size(X);ifnd,disp('Inputdatamustbenxd!/n');returnenderr_X=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofVerify_X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionerr_k=Verify_k(k)err_k=1;if~isnumeric(k)|~isreal(k)|k1,disp('kmustbearealinteger=1!/n');returnenderr_k=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofVerify_k%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[ltol,err_ltol]=Verify_ltol(ltol)err_ltol=1;ifisempty(ltol),ltol=0.1;elseif~isreal(ltol)|ltol=0,disp('ltolmustbeapositiverealnumber!');returnenderr_ltol=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofVerify_ltol%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[maxiter,err_maxiter]=Verify_maxiter(maxiter)err_maxiter=1;ifisempty(maxiter),maxiter=1000;elseif~isreal(maxiter)|maxiter=0,disp('ltolmustbeapositiverealnumber!');returnenderr_maxiter=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofVerify_maxiter%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[pflag,err_pflag]=Verify_pflag(pflag)err_pflag=1;ifisempty(pflag),pflag=0;elseifpflag~=0&pflag~=1,disp('Plotflagmustbeeither0or1!/n');returnenderr_pflag=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EndofVerify_pflag%%%%%%%%%%%%%%%%%%%%%%%%%%%