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的话,明显还是自己写的代码拟合效果更好一些啊,有些不科学。

算法学习笔记(4):计数排序

在宿舍呆着无聊接着在hackerrank上刷水题。这道题由于排序的依据是一定范围内的整数,所以用计数排序效率应当不错。

说到计数排序,实际上思路也很常见。最常见的一个例子就是比赛或者是考试的时候要进行排名。

比如某班有3名同学考了100分,2名同学考了99分,5名同学考了98分。那么很明显,比99分高的同学一共有3名,所以考99分的同学只能是第4名了。

计数排序大概就是做了个数数每个元素前面有多少元素这个事情。先统计最高分多少分,然后依次统计各个分数的学生个数,然后以此叠加就可以统计出前面有多少个学生了。

而且由于不同分数段的学生是按照原顺序放置到新数组中的,所以在原数组中的顺序并不会被打乱,即这个算法是稳定的。这也是为什么这道题可以用计数排序解,因为并不会打乱文本的顺序。

顺便一提,这道题中的最大值max()的复杂度应该是,所以总运行时间应该是。因为计数排序并不涉及到不同元素之间的比较,即不是比较排序,所以的下界并不适用于计数排序。

最后附上这道题的一个无脑解。

# Enter your code here. Read input from STDIN. Print output to STDOUT
n = int(raw_input())
ar = []
idx = []
out = [None]*n
for _ in range(n):
    ar.append(raw_input().split())
    ar[-1][0] = int(ar[-1][0])
    idx.append(ar[-1][0])
    if _<n/2:
        ar[-1][1] = "-"
countlist = [0]*(max(idx)+1)
for item in ar:
    countlist[item[0]]+=1
for i in range(len(countlist)-1):
    countlist[i+1]+=countlist[i]
for item in reversed(ar):
    out[countlist[item[0]]-1]=item
    countlist[item[0]]-=1
for item in out:
    print item[1],