EM算法(期望最大算法)实现

close all; clear; clc; %% 这个是已经改进了的EM算法通过统计找到中心点再进行迭代。 % 参考书籍Pattern.Recognition.and.Machine.Learning.pdf % http://www.pr-ml.cn % [email protected] % 2009/10/15 %% M=2; % number of Gaussian N=50000; % total number of data samples th=0.1; % convergent threshold Nit=200; % maximal iteration Nrep=10; % number of repetation to find global maximal jd=0.1; % the step K=2; % demention of output signal ptime=0.01; pi=3.141592653589793; % in case it is overwriten by smae name variable cond_num =100; % prevent the singular covariance matrix in simulation data %% paramethers for random signal genrator a_real =[1/2;1/2];%这里需要手工设置 mu_real=[3 7; 7 3]; cov_real(:,:,1)=[0.5 0; 0 0.5]; cov_real(:,:,2)=[0.5 0; 0 0.5]; %% generate the data%这里生成的数据全部符合标准 x=[mvnrnd(mu_real(:,1),cov_real(:,:,1),N*a_real(1))',mvnrnd(mu_real(:,2),cov_real(:,:,2),N*a_real(2))']; % num=0; for i=1:N*a_real(1) while (~((x(1,i)>0)&&(x(2,i)>0)&&(x(1,i)<10)&&(x(2,i)<10))) x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'; % num=num+1; end end for i=N*a_real(1)+1:N while (~((x(1,i)>0)&&(x(2,i)>0)&&(x(1,i)<10)&&(x(2,i)<10))) x(:,i)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'; % num=num+1; end end %这里生成的数据全部符合标准 %% %%%%%%%%%%%%%%%% 参数初始化 a=[1/3,2/3]; mu=[2 4;4 7];%均值初始化完毕 cov(:,:,1)=[5 0; 0 0.5]; cov(:,:,2)=[5 0; 0 0.5];%协方差初始化 %% EM Algorothm % loop f_best=-inf; for crep=1:Nrep fprintf('第 %04d 个局部最大值\n',crep); while 1 a_old = a; mu_old = mu; cov_old= cov; rznk_p=zeros(M,N); for cm=1:M mu_cm=mu(:,cm); cov_cm=cov(:,:,cm); for cn=1:N p_cm=exp(-0.5*(x(:,cn)-mu_cm)'/cov_cm*(x(:,cn)-mu_cm)); rznk_p(cm,cn)=p_cm; end rznk_p(cm,:)=rznk_p(cm,:)/sqrt(det(cov_cm)); end rznk_p=rznk_p*(2*pi)^(-K/2); %E step %开始求rznk rznk=zeros(M,N);%r(Z pikn=zeros(1,M);%r(Z pikn_sum=0; for cn=1:N for cm=1:M pikn(1,cm)=a(cm)*rznk_p(cm,cn); % pikn_sum=pikn_sum+pikn(1,cm); end for cm=1:M rznk(cm,cn)=pikn(1,cm)/sum(pikn); end end %求rank结束 % M step nk=zeros(1,M); for cm=1:M for cn=1:N nk(1,cm)=nk(1,cm)+rznk(cm,cn); end end a=nk/N; rznk_sum_mu=zeros(M,1); % 求均值MU for cm=1:M rznk_sum_mu=0;%开始的时候就是错在这里，这里要置零。 for cn=1:N rznk_sum_mu=rznk_sum_mu+rznk(cm,cn)*x(:,cn); end mu(:,cm)=rznk_sum_mu/nk(cm); end % 求 协方差COV for cm=1:M rznk_sum_cov=zeros(K,M); for cn=1:N rznk_sum_cov=rznk_sum_cov+rznk(cm,cn)*(x(:,cn)-mu(:,cm))*(x(:,cn)-mu(:,cm))'; end cov(:,:,cm)=rznk_sum_cov/nk(cm); end t=max([norm(a_old(:)-a(:))/norm(a_old(:));norm(mu_old(:)-mu(:))/norm(mu_old(:));norm(cov_old(:)-cov(:))/norm(cov_old(:))]); disp(t); if t<th break; end end %while 1 f=sum(log(sum(pikn))); if f>f_best a_best=a; mu_best=mu; cov_best=cov; f_best=f; end end % for crep=1:Nrep %% 输出结果 disp('a_best='); disp(a_best); disp('mu_best='); disp(mu_best); disp('cov_best='); disp(cov_best);