数据挖掘算法总结-EM算法

 数据挖掘算法总结-EM算法

 

作者：刘伟民

毕业于：中科院计算所

职业：搜索引擎爱好者

2011-05-07 08:54:52

 

 

1.  什么是EM算法？

EM算法是机器学习中一个很重要的算法，即期望最大化算法，主要包括以下两个步骤：

E步骤：estimate the expected values

M步骤：re-estimate parameters

这个算法的主要作用在于对参数的估计上。虽然EM算法也可以进行数据聚类，并且基于混合高斯分布进行数据拟合，但是由于EM算法进行迭代速度很慢，比kmeans性能差很多，并且KMEANS算法 聚类效果没有比EM差多少，所以一般用kmeans进行聚类，而不是EM。

2.   EM算法简述

可以有一些比较形象的比喻说法把这个算法讲清楚。比如说食堂的大师傅炒了一份菜，要等分成两份给两个人吃，显然没有必要拿来天平一点一点的精确的去称分量，最简单的办法是先随意的把菜分到两个碗中，然后观察是否一样多，把比较多的那一份取出一点放到另一个碗中，这个过程一直迭代地执行下去，直到大家看不出两个碗所容纳的菜有什么分量上的不同为止。EM算法就是这样，假设我们估计知道A和B两个参数，在开始状态下二者都是未知的，并且知道了A的信息就可以得到B的信息，反过来知道了B也就得到了A。可以考虑首先赋予A某种初值，以此得到B的估计值，然后从B的当前值出发，重新估计A的取值，这个过程一直持续到收敛为止。

EM 算法是 Dempster，Laind，Rubin 于 1977 年提出的求参数极大似然估计的一种方法，它可以从非完整数据集中对参数进行 MLE 估计，是一种非常简单实用的学习算法。这种方法可以广泛地应用于处理缺损数据，截尾数据，带有噪声等所谓的不完全数据(incomplete data)。 假定集合Z = (X,Y)由观测数据 X 和未观测数据Y 组成，Z = (X,Y)和 X 分别称为完整数据和不完整数据。假设Z的联合概率密度被参数化地定义为P(X，Y|Θ)，其中Θ 表示要被估计的参数。Θ 的最大似然估计是求不完整数据的对数似然函数L(X;Θ)的最大值而得到的： L(Θ; X )= log p(X |Θ) = ∫log p(X ,Y |Θ)dY ； EM算法包括两个步骤：由E步和M步组成，它是通过迭代地最大化完整数据的对数似然函数Lc( X;Θ )的期望来最大化不完整数据的对数似然函数，其中Lc(X;Θ) =log p(X，Y |Θ) ； 假设在算法第t次迭代后Θ 获得的估计记为Θ(t ) ，则在（t+1）次迭代时， E-步：计算完整数据的对数似然函数的期望，记为： 　Q(Θ |Θ (t) ) = E{Lc(Θ;Z)|X;Θ(t) }； M-步：通过最大化Q(Θ |Θ(t) ) 来获得新的Θ 。

通过交替使用这两个步骤，EM算法逐步改进模型的参数，使参数和训练样本的似然概率逐渐增大，最后终止于一个极大点。直观地理解EM算法，它也可被看作为一个逐次逼近算法：事先并不知道模型的参数，可以随机的选择一套参数或者事先粗略地给定某个初始参数λ0 ，确定出对应于这组参数的最可能的状态，计算每个训练样本的可能结果的概率，在当前的状态下再由样本对参数修正，重新估计参数λ ，并在新的参数下重新确定模型的状态，这样，通过多次的迭代，循环直至某个收敛条件满足为止，就可以使得模型的参数逐渐逼近真实参数。

EM算法的主要目的是提供一个简单的迭代算法计算后验密度函数，它的最大优点是简单和稳定，但容易陷入局部最优。

3.  EM算法的缺点

EM算法迭代速度慢，次数多。

EM算法对高位数据处理速度相当慢，并且在估计协方差的时候不够准确。

（话说EM算法主要用于参数估计，协方差只是在基于混合高斯分布的EM算法中存在）

4.  下面附我用matlab实现的EM代码

  close all; clear; clc; % 作者：刘伟民 % 单位：中科院计算所 % 邮箱：liuwm@ics.ict.ac.cn % 功能：最经典的二维两个高斯类的EM算法实现 % 参考书籍Pattern.Recognition.and.Machine.Learning.pdf % 时间20091118 %% 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 cn=1:N*a_real(1) while (~((x(1,cn)>0)&&(x(2,cn)>0)&&(x(1,cn)<10)&&(x(2,cn)<10))) x(:,cn)=mvnrnd(mu_real(:,1),cov_real(:,:,1),1)'; % num=num+1; end end for cn=N*a_real(1)+1:N while (~((x(1,cn)>0)&&(x(2,cn)>0)&&(x(1,cn)<10)&&(x(2,cn)<10))) x(:,cn)=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);