EM_GMM.rar_EM估计GMM参数_EM改进算法_GMM参数估计_改进EM_改进的EM

clear all;clc; X=[]; N=1000; mix_means=[6 0 -6]; mix_sigm2=[4,1,1]; mix_priors=[0.5,0.3,0.2]; for n=1:N comp=find(rand<cumsum(mix_priors)); comp=comp(1); X(n)=mvnrnd(mix_means(comp),mix_sigm2(comp)); end X=X'; figure(1); scatter(1:N,X,'fill'); %% 初始化参数 K = 3; % 隐变量的取值有3种可能情况，此处可理解为观测数据来自由三个高斯分布混合而成的混合高斯模型 MaxIts=50;%允许的最大迭代次数 means = rand(K,1); sigm2 = rand(K,1);%更为可取的初始化高斯混合模型的方式是用K-means算法，会减少迭代次数 priors = repmat(1/K,K,1); %% 运行算法 q = zeros(N,K);%初始化E步所需计算的隐变量的条件概率，即伽马函数 B(1) = -inf;%初始化对数似然函数的下界，用来判断收敛的 converged = 0; it = 0;%it是进行迭代（iterate）的次数 tol = 1e-1;%判断收敛与否设定的阀值 while it<MaxIts it = it + 1; % E步，计算隐变量的后验概率q for k=1:K temp(:,k)=priors(k,it)*(1/(2*pi*sigm2(k,it))^0.5)*exp(-((X-repmat(means(k,it),N,1)).^2)/(2*sigm2(k,it))); end q=temp./repmat(sum(temp,2),1,K); % 计算似然函数的下界并判断是否收敛 % if it>1 % B(it) = sum(sum(q.*log(temp)))-sum(sum(q.*log(q))); % if abs(B(it)-B(it-1))<tol % converged = 1; % end % if converged == 1 % break % end % end %%M步 % 更新隐含变量的先验概率 a = mean(q,1);%对N*K维矩阵q的每一列求平均，返回一个1*K维的行向量 priors(:,it+1)=a'; % 更新均值 for k = 1:K means(k,it+1) = sum(X.*q(:,k),1)/sum(q(:,k),1); end % 更新方差 for k = 1:K b=0 for n=1:N b=q(n,k)*(X(n)-means(k,it+1))^2+ b; end sigm2(k,it+1) = b/sum(q(:,k),1); end end figure(2); for k=1:K subplot(1,3,k); plot(0:(it-1),means(k,1:it),'-k.'); ti=sprintf('mean %g',k); title(ti); end; figure(3); for k=1:K subplot(1,3,k); plot(0:(it-1),sigm2(k,1:it),'-k.'); ti=sprintf('sigma %g',k); title(ti); end; figure(4); for k=1:K subplot(1,3,k); plot(0:(it-1),priors(k,1:it),'-k.'); ti=sprintf('priors %g',k); title(ti); end; figure(5); plot(2:length(B),B(2:end),'k'); xlabel('Iterations'); ylabel('Bound');