GMM的EM算法实现

GMM的全名是Gaussian Mixture Model,即高斯混合模型,这个名字挺直观。

如果给了一组数据点,假设这些数据点符合iid,那么很容就能根据数据拟合出一个最佳的高斯分布。

但是如果给了一堆数据点,而他们却不符合iid,只能做到互相独立,并不能做到同分布。相反的,数据点们有可能来自不同的高斯分布,这也就是『高斯混合』模型名字的意思。这样如何能拟合出最佳的高斯分布们呢?

EM算法提供了一种方法,这种方法的思路有些像牛拉法。先是确定一个初值,然后看看当前值距离目标还有多远,然后逐渐的缩小差距(不断修正membership weight)。

具体原理见这里

这篇小短文大概记录下这次作业的内容。作业题目是生成了2个二维高斯分布混合的数据点,共2000个,然后用这两千个点拟合出2个二维高斯分布的参数。

clear;
%产生2个二维正态数据
MU1    = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2    = [-1 -1];
SIGMA2 = [1 0; 0 1];
X      = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'.');

%initial values
h_mu1 = X(randi(2000),:);
h_sigma1 = cov(X);
h_mu2 = X(randi(2000),:);
h_sigma2 = cov(X);
h_alpha1 = 0.5;
h_alpha2 = 0.5;

%simple iteration
for iter=1:100
    %E step
    w1 = mvnpdf(X,h_mu1,h_sigma1)*h_alpha1./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);
    w2 = mvnpdf(X,h_mu2,h_sigma2)*h_alpha2./...
    (mvnpdf(X,h_mu1,h_sigma1)*h_alpha1 + mvnpdf(X,h_mu2,h_sigma2)*h_alpha2);

    %M step
    h_alpha1 = sum(w1)/2000;
    h_alpha2 = sum(w2)/2000;
    
    h_mu1 = sum([w1.*X(:,1) w1.*X(:,2)])/sum(w1);
    h_mu2 = sum([w2.*X(:,2) w2.*X(:,2)])/sum(w2);

    M1 = X-ones(2000,1)*h_mu1;
    h_sigma1 = ([sum(w1.*M1(:,1).^2) 0;0 sum(w1.*M1(:,2).^2)]+...
        sum(w1.*(M1(:,1).*M1(:,2)))*[0 1;1 0])/sum(w1);
    M2 = X-ones(2000,1)*h_mu2;
    h_sigma2 = ([sum(w2.*M2(:,1).^2) 0;0 sum(w2.*M2(:,2).^2)]+...
        sum(w2.*(M2(:,1).*M2(:,2)))*[0 1;1 0])/sum(w2);
end

figure;
[x,y] = meshgrid(linspace(-5,5,75));
p1 = mvnpdf([x(:),y(:)],h_mu1,h_sigma1);
p2 = mvnpdf([x(:),y(:)],h_mu2,h_sigma2);
surf(x,y,reshape(p1,75,75));hold on
surf(x,y,reshape(p2,75,75));

options = statset('Display','final');
obj = gmdistribution.fit(X,2,'Options',options);
figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);

绝大多数时候都没有什么问题,不过由于初始点是从X里随便选的一个数据点,极偶然会跑偏。

以下分别是数据点的分布、自行拟合出的结果、利用matlab自带拟合函数拟合出的结果。

originalgmm

根据给定的均值、协方差矩阵画出的GMM的图

xscatter

数据点的散点图,保存之后可能有点糊

fitbymyself

自己写的脚本拟合出的GMM的图

fitbymatlabtoolbox

matlab提供的拟合函数拟合出的图

迭代100次之后,高斯分布1的均值分别偏差0.2776%、4.4462%,协方差分别偏差了2.4458%、5.0146%;高斯分布2的均值偏差了2.6981%、1.2911%,协方差偏差了3.8498%、1.5517%。精度凑凑活活把...

对比之下,matlab自带的拟合函数拟合出来的高斯分布1的均值偏差2.2415%、2.2415%,协方差偏差为4.02%、1.32%;高斯分布2的均值偏差了0.4165%、0.4165%,协方差偏差为2.39%、4.62%。精度也没高到哪里去...不过均值偏差一样让我很是在意...

然而看图1、图3和图4的话,明显还是自己写的代码拟合效果更好一些啊,有些不科学。