算法特点

GMM聚类与Keams聚类很相似。K-means是将每个数据点分配到某一个类别中，而GMM则是给出这些数据点被分配到每个类别的概率，又叫做soft assignment。其除了被用在clustering上，还经常被用于density estimation上。

得出一个概率有很多的好处，因为它的信息量比简单的一个结果要多。这个概率可以被看做是得出这个结果的把握。比如说，在诸如疾病诊断时，算法得出49%的概率患病，为51%的概率正常，如果仅仅用50%的阈值将患者诊断为“正常”的话，风险是非常大的，也就是说得出正常的结果的把握是非常小的。这个情况下，就会拒绝发表评论，而把这个任务留给有经验的医生去解决。

学习的过程可以看做是一种“归纳”的过程，在归纳的时候需要有一个提前假设，即“归纳偏置”。比如说线性回归，它的归纳偏执就是“数据是符合线性函数的，满足要求的函数必须是线性函数”。如果去掉“线性函数”这个归纳偏执，那么就会产生过拟合的问题，我们就可以构造n-1次多项式函数，让它完美的穿过所有的N的点，甚至可以构造出无数个满足条件的函数。然而，我们学习的目的通常不是解释现有的事物，而是从中归纳知识，并且能够应用到新的事物上。

没有归纳偏执或者归纳偏执太宽泛会导致过拟合，而限制过大的归纳偏执也是有问题的；比如说，如果数据不是非线性分布的，强行的使用非线性函数去做分类通常并不能够得到好的结果。难点就是在这之间找到一个平衡点。

混合高斯模型的归纳偏执很明显，就是假设数据符合混合高斯分布，也就是看作数据是从多个高斯分布中生成而来。混合模型其实是可以变的任意复杂的，通过增加模型的个数，我们可以任意的逼近任何连续的概率密度分布。

主程序：

clear;

clc;

%从实验结果可以看出先使用Keams方法计算初始点，再使用GMM进行聚类得到的结果

%明显优于直接使用GMM方法进行聚类得到的结果

%数据X1和X2的真实均值为[1,1;4,4;8,1]

%协方差矩阵均为[2,0;0,2];

%数据集X1的混合系数为[1/3,1/3,1/3]

%数据集X2的混合系数为[0.6,0.3,0.1]

[X1,X2]=generateData();

K=3;

%数据集X1

[a_init,mu_init,sigma_init]=initPara(X1);

[a_GMM,mu_GMM,sigma_GMM]=GMM(X1,a_init,mu_init,sigma_init);

[centroids,~] = Kmeans( X1,mu_init,K);

disp('1.Keams算法执行后数据集X1的均值为：')

disp(centroids);

disp('2.GMM后高斯混合模型数据集X1的均值，协方差矩阵为：');

disp('均值为：');

disp(mu_GMM);

disp('协方差矩阵为：')

disp(sigma_GMM);

disp(['混合系数为：',num2str(a_GMM)]);

[a_Kea_GMM,mu_Kea_GMM,sigma_Kea_GMM]=GMM(X1,a_init,centroids,sigma_init);

disp('3.使用Keams计算初始均值再使用GMM进行聚类后的数据集X1的均值，协方差矩阵为：');

disp('均值为：');

disp(mu_Kea_GMM);

disp('协方差矩阵为：')

disp(sigma_Kea_GMM);

disp(['混合系数为：',num2str(a_Kea_GMM)]);

%数据集X2

[a_init,mu_init,sigma_init]=initPara(X2);

[a_fin,mu_fin,sigma_GMM]=GMM(X2,a_init,mu_init,sigma_init);

[centroids,~] = Kmeans( X2,mu_init,K);

disp('1.Keams算法执行后数据集X2的均值为：')

disp(centroids);

disp('2.GMM后高斯混合模型数据集X2的均值，协方差矩阵为：');

disp('均值为：');

disp(mu_GMM);

disp('协方差矩阵为：')

disp(sigma_GMM);

disp(['混合系数为：',num2str(a_GMM)]);

[a_Kea_GMM,mu_Kea_GMM,sigma_Kea_GMM]=GMM(X2,a_init,centroids,sigma_init);

disp('3.使用Keams计算初始均值再使用GMM进行聚类后的数据集X2的均值，协方差矩阵为：');

disp('均值为：');

disp(mu_Kea_GMM);

disp('协方差矩阵为：')

disp(sigma_Kea_GMM);

disp(['混合系数为：',num2str(a_Kea_GMM)]);

首先根据均值协方差矩阵混合系数生成高斯混合分布数据：

function [X1,X2]=generateData()

mu=[1,1;4,4;8,1];

sigma=[2,0;0,2];

a1=[1/3,1/3,1/3];

a2=[0.6,0.3,0.1];

N=1000;

rand1=randsrc(N,1,[[1,2,3];a1]);

rand2=randsrc(N,1,[[1,2,3];a2]);

X1=[];

X2=[];

mean1=[];

for i=1:size(a1,2)

X1_temp=mvnrnd(mu(i,:),sigma,length(find(rand1==i)));

X1=[X1;X1_temp];

subplot(1,2,1);

plot(X1_temp(:,1),X1_temp(:,2),'+');

title('X1');

legend('模型1','模型2','模型3');

xlabel('x');ylabel('y');

hold on;

mean1=[mean1;mean(X1_temp)];

cov1(:,:,i)=cov(X1_temp);

X2_temp=mvnrnd(mu(i,:),sigma,length(find(rand2==i)));

X2=[X2;X2_temp];

subplot(1,2,2);

plot(X2_temp(:,1),X2_temp(:,2),'*');

title('X2');

legend('模型1','模型2','模型3');

xlabel('x');ylabel('y');

hold on;

end

disp(['直接使用数据集X1来估计的属于三个模型随机矢量的均值分别为：']);

disp(mean1);

disp(['直接使用数据集X1来的属于三个模型随机矢量的协方差矩阵分别为：']);

disp(cov1);

end

初始化参数：

function [ a,mu,sigma ] = initPara(X)

[m,n]=size(X);

r=randperm(m);

mu=X(r(1:3),:);

a=[1/3,1/3,1/3];

sigma=[1,0;0,1];

end

GMM算法：

function [a_fin,mu_fin,sigma_fin]=GMM(X1,a_init,mu_init,sigma_init)

%UNTITLED2 此处显示有关此函数的摘要

% 此处显示详细说明

K=size(a_init,2);

[M,N]=size(X1);

a=a_init;

mu=mu_init;

sigma=sigma_init;

px=zeros(M,K);

thre=1e-5;

LLD_pro=inf;

while true

for i=1:K

px(:,i)=mvnpdf(X1,mu(i,:),sigma);

end

pGramm=repmat(a,M,1).*px;

pGramm=pGramm./repmat(sum(pGramm,2),1,K);

Nk=sum(pGramm,1); %1*K

mu=diag(1./Nk)*pGramm'*X1;

for kk=1:K

Xshift=X1-repmat(mu(kk,:),M,1);

sigma(:,:,kk)=(Xshift'*diag(pGramm(:,kk))*Xshift)/Nk(:,kk);

end

a=Nk/M;

LLD=sum(log(a*px'));

if (LLD-LLD_pro)

break;

end

LLD_pro=LLD;

end

a_fin=a;

mu_fin=mu;

sigma_fin=sigma;

end

Keams：

function [ centroids,centroidsIndex ] = Kmeans( X,mu_init,K)

centroids=mu_init;

thre=1e-5;

while true

pro_centroids=centroids;

[centroidsIndex]=findCloestCentroids(X,centroids,K);

[centroids]=computeCentroids(centroidsIndex,X,K);

if norm(centroids-pro_centroids)

break;

end

end

function [centroidsIndex]=findCloestCentroids(X,centroids,K)

[M,~]=size(X);

distMat=zeros(M,K);

centroidsIndex=zeros(M,1);

for i=1:M

for j=1:K

distMat(i,j)=norm(X(i,:)-centroids(j,:));

end

centroidsIndex(i,:)=find(distMat(i,:)==min(distMat(i,:)));

end

end

function [centroids]=computeCentroids(centroidsIndex,X,K)

[M,N]=size(X);

centroids=zeros(K,N);

for i=1:K

centroidsNum=sum(centroidsIndex==i);

centroids(i,:)=sum(X(centroidsIndex==i,:))./centroidsNum;

end

end

算法问题

GMM和K-means有同样的问题，如果初始值选择的不好，那么可能会得到局部最优解，得不到全局最优解。我们可以先使用K-means方法得到的结果作为初始值传给GMM模型，再使用GMM进行细致迭代。 GMM本身只是一个模型，EM算法并不是唯一的求解方法。