[zz] 混合高斯模型 Gaussian Mixture Model

聚类(1)——混合高斯模型 Gaussian Mixture Model

http://blog.csdn.net/jwh_bupt/article/details/7663885

聚类系列：

--------------------------------

    聚类的方法有很多种，k-means要数最简单的一种聚类方法了，其大致思想就是把数据分为多个堆，每个堆就是一类。每个堆都有一个聚类中心（学习的结果就是获得这k个聚类中心），这个中心就是这个类中所有数据的均值，而这个堆中所有的点到该类的聚类中心都小于到其他类的聚类中心（分类的过程就是将未知数据对这k个聚类中心进行比较的过程，离谁近就是谁）。其实k-means算的上最直观、最方便理解的一种聚类方式了，原则就是把最像的数据分在一起，而“像”这个定义由我们来完成，比如说欧式距离的最小，等等。想对k-means的具体算法过程了解的话，请看这里。而在这篇博文里，我要介绍的是另外一种比较流行的聚类方法----GMM（Gaussian Mixture Model）。

    GMM和k-means其实是十分相似的，区别仅仅在于对GMM来说，我们引入了概率。说到这里，我想先补充一点东西。统计学习的模型有两种，一种是概率模型，一种是非概率模型。所谓概率模型，就是指我们要学习的模型的形式是P(Y|X)，这样在分类的过程中，我们通过未知数据X可以获得Y取值的一个概率分布，也就是训练后模型得到的输出不是一个具体的值，而是一系列值的概率（对应于分类问题来说，就是对应于各个不同的类的概率），然后我们可以选取概率最大的那个类作为判决对象（算软分类soft assignment）。而非概率模型，就是指我们学习的模型是一个决策函数Y=f(X)，输入数据X是多少就可以投影得到唯一的一个Y，就是判决结果（算硬分类hard assignment）。回到GMM，学习的过程就是训练出几个概率分布，所谓混合高斯模型就是指对样本的概率密度分布进行估计，而估计的模型是几个高斯模型加权之和（具体是几个要在模型训练前建立好）。每个高斯模型就代表了一个类（一个Cluster）。对样本中的数据分别在几个高斯模型上投影，就会分别得到在各个类上的概率。然后我们可以选取概率最大的类所为判决结果。

    得到概率有什么好处呢？我们知道人很聪明，就是在于我们会用各种不同的模型对观察到的事物和现象做判决和分析。当你在路上发现一条狗的时候，你可能光看外形好像邻居家的狗，又更像一点点女朋友家的狗，你很难判断，所以从外形上看，用软分类的方法，是女朋友家的狗概率51%，是邻居家的狗的概率是49%，属于一个易混淆的区域内，这时你可以再用其它办法进行区分到底是谁家的狗。而如果是硬分类的话，你所判断的就是女朋友家的狗，没有“多像”这个概念，所以不方便多模型的融合。

    从中心极限定理的角度上看，把混合模型假设为高斯的是比较合理的，当然也可以根据实际数据定义成任何分布的Mixture Model,不过定义为高斯的在计算上有一些方便之处，另外，理论上可以通过增加Model的个数，用GMM近似任何概率分布。

    混合高斯模型的定义为：

   

    其中K为模型的个数，πk为第k个高斯的权重，则为第k个高斯的概率密度函数，其均值为μk，方差为σk。我们对此概率密度的估计就是要求πk、μk和σk各个变量。当求出的表达式后，求和式的各项的结果就分别代表样本x属于各个类的概率。

    在做参数估计的时候，常采用的方法是最大似然。最大似然法就是使样本点在估计的概率密度函数上的概率值最大。由于概率值一般都很小，N很大的时候这个连乘的结果非常小，容易造成浮点数下溢。所以我们通常取log，将目标改写成：

  

    也就是最大化log-likelyhood function，完整形式则为：

    一般用来做参数估计的时候，我们都是通过对待求变量进行求导来求极值，在上式中，log函数中又有求和，你想用求导的方法算的话方程组将会非常复杂，所以我们不好考虑用该方法求解（没有闭合解）。可以采用的求解方法是EM算法——将求解分为两步：第一步是假设我们知道各个高斯模型的参数（可以初始化一个，或者基于上一步迭代结果），去估计每个高斯模型的权值；第二步是基于估计的权值，回过头再去确定高斯模型的参数。重复这两个步骤，直到波动很小，近似达到极值（注意这里是个极值不是最值，EM算法会陷入局部最优）。具体表达如下：

  

    1、对于第i个样本xi来说，它由第k个model生成的概率为：

   

    在这一步，我们假设高斯模型的参数和是已知的（由上一步迭代而来或由初始值决定）。

   （E step）

   

    （M step）

    3、重复上述两步骤直到算法收敛（这个算法一定是收敛的，至于具体的证明请回溯到EM算法中去，而我也没有具体关注，以后补上）。

    最后总结一下，用GMM的优点是投影后样本点不是得到一个确定的分类标记，而是得到每个类的概率，这是一个重要信息。GMM每一步迭代的计算量比较大，大于k-means。GMM的求解办法基于EM算法，因此有可能陷入局部极值，这和初始值的选取十分相关了。GMM不仅可以用在聚类上，也可以用在概率密度估计上。

http://www.cnblogs.com/Key-Ky/p/3605779.html

混合高斯模型

玩了混合高斯模型，先转几个参考资料，曾经试过自己写代码，结果发现混合高斯模型矩阵运算对我的计算能力要求很高，结果失败了，上网找了代码学习一下大牛们的编程思想，事实证明数学写出来的公式虽然很美，但是现实写代码的时候要考虑各种问题~~~

1.http://www.cnblogs.com/cfantaisie/archive/2011/08/20/2147075.html （主要实现代码）

2.http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/ （对于奇异矩阵的正则化）

3.参考的一本书，Computer Vision Models, Learning, and Inference.（数学推导）

4.斯坦福机器学习课程（数学推导）

利用EM算法：

E-step：

M-step:

matlab实现代码：

复制代码
function prob = MOF_guassPdf(Data,Mu,Sigma)
%

% 根据高斯分布函数计算每组数据的概率密度 Probability Density Function (PDF) 

% 输入 -----------------------------------------------------------------

%   o Data:  D x N ，N个D维数据

%   o Mu:    D x 1 ，M个Gauss模型的中心初始值

%   o Sigma: D x D ，每个Gauss模型的方差（假设每个方差矩阵都是对角阵，

%                                   即一个数和单位矩阵的乘积）

% Outputs ----------------------------------------------------------------

%   o prob:  N x 1 array representing the probabilities for the 

%            N datapoints.     

[dim,N] = size(Data);

Data = Data' - repmat(Mu',N,1);

prob = sum((Data/Sigma).*Data, 2);

prob = exp(-0.5*prob) / sqrt((2*pi)^dim * (abs(det(Sigma))+realmin));
复制代码

EM迭代：

复制代码
function [Alpha, Mu, Sigma] = MOF_EM(Data, Alpha0, Mu0, Sigma0)
%% EM 迭代停止条件
loglik_threshold = 1e-10;

%% 初始化参数

[dim, N] = size(Data); 

M = size(Mu0,2);

loglik_old = -realmax;

nbStep = 0;
                %Data是D*N的矩阵

Mu = Mu0;       %Mu是D*M的矩阵

Sigma = Sigma0; %sigma是D*D*M的三维矩阵

Alpha = Alpha0; %alpha是1*M大小的向量，意思是列代表第M个高斯模型。

Epsilon = 0.0001; 

while (nbStep < 1200)

  nbStep = nbStep+1;

  %% E-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % PDF of each point

    Pxi(:,i) = MOF_guassPdf(Data, Mu(:,i), Sigma(:,:,i));          

  end

  

  % 计算后验概率 beta(i|x)

  Pix_tmp = repmat(Alpha,[N 1]).*Pxi; %Pxi应该是N*M的高斯概率矩阵，意思是第N个数据计算第M个高斯函数的概率

  Pix = Pix_tmp ./ (repmat(sum(Pix_tmp,2),[1 M])+realmin);

  Beta = sum(Pix);

  %% M-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % 更新权值

    Alpha(i) = Beta(i) / N;

    % 更新均值

    Mu(:,i) = Data*Pix(:,i) / Beta(i);

    % 更新方差

    Data_tmp1 = Data - repmat(Mu(:,i),1,N);

    Sigma(:,:,i) = (repmat(Pix(:,i)',dim, 1) .* Data_tmp1 * Data_tmp1') / Beta(i);

    %% Add a tiny variance to avoid numerical instability

    Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(dim,1));

  end

 

%  %% Stopping criterion 1 %%%%%%%%%%%%%%%%%%%%

  for i=1:M

    %Compute the new probability p(x|i)

    Pxi(:,i) = MOF_guassPdf(Data, Mu(:,i), Sigma(i));

  end

  %Compute the log likelihood

  F = Pxi*Alpha';

  F(find(F<realmin)) = realmin;

  loglik = mean(log(F));

  %Stop the process depending on the increase of the log likelihood 

  if abs((loglik/loglik_old)-1) < loglik_threshold

    break;

  end

  loglik_old = loglik;

 

  %% Stopping criterion 2 %%%%%%%%%%%%%%%%%%%%
%{
  v = [sum(abs(Mu - Mu0)), abs(Alpha - Alpha0)];

  s = abs(Sigma-Sigma0);

  v2 = 0;

  for i=1:M

    v2 = v2 + det(s(:,:,i));

  end

 

  if ((sum(v) + v2) < Epsilon)

    break;

  end

  Mu0 = Mu;

  Sigma0 = Sigma;

  Alpha0 = Alpha;
%}
end
复制代码

测试结果

编写代码随机生成3个高斯分布的数据，参数也随机生成（注意sigma要半正定对称）：

这里只画了一个经过EM算法迭代所得参数的高斯函数图..（有点丑，不知道怎么将mesh弄透明，求搜索关键词）。

高斯混合模型（GMM）参数优化及实现  

http://www.cnblogs.com/cfantaisie/archive/2011/08/20/2147075.html

高斯混合模型（GMM）参数优化及实现 （< xmlnamespace prefix ="st1" ns ="urn:schemas-microsoft-com:office:smarttags" />2010-11-13）

1 高斯混合模型概述< xmlnamespace prefix ="o" ns ="urn:schemas-microsoft-com:office:office" />

高斯密度函数估计是一种参数化模型。有单高斯模型（Single Gaussian Model, SGM）和高斯混合模型（Gaussian mixture model，GMM）两类。类似于聚类，根据高斯概率密度函数（PDF，见公式1）参数的不同，每一个高斯模型可以看作一种类别，输入一个样本< xmlnamespace prefix ="v" ns ="urn:schemas-microsoft-com:vml" /> ，即可通过PDF计算其值，然后通过一个阈值来判断该样本是否属于高斯模型。很明显，SGM适合于仅有两类别问题的划分，而GMM由于具有多个模型，划分更为精细，适用于多类别的划分，可以应用于复杂对象建模。

下面以视频前景分割应用场景为例，说明SGM与GMM在应用上的优劣比较：

l        SGM需要进行初始化，如在进行视频背景分割时，这意味着如果人体在前几帧就出现在摄像头前，人体将会被初始化为背景，而使模型无法使用；

l        SGM只能进行微小性渐变，而不可突变。如户外亮度随时间的渐变是可以适应的，如果在明亮的室内突然关灯，单高斯模型就会将整个室内全部判断为前景。又如，若在监控范围内开了一辆车，并在摄像头下开始停留。由于与模型无法匹配，车会一直被视为前景。当车过很长时间离去时，由于车停留点的亮度发生了很大的变化，因此已经无法与先前的背景模型相匹配；

l        SGM无法适应背景有多个状态，如窗帘，风吹的树叶。单高斯模型无法表示这种情况，而使得前背景检测混乱，而GMM能够很好地描述不同状态；

l        相对于单高斯模型的自适应变化，混合高斯模型的自适应变化要健壮的多。它能解决单高斯模型很多不能解决的问题。如无法解决同一样本点的多种状态，无法进行模型状态转化等。

2 理论说明部分

因博客中无法编辑公式，故详细文档见这里。代码如下：

3 源码

3.1 单高斯模型

下面代码实现了SGM，并实现了人脸肤色检测。其中图像处理、矩阵运算采用了openCV库函数

/*****************************************************************************

       Single Gaussian Model for skin color extraction

       Param:

              img -- input image to extract the face region

              skinImg -- result

*****************************************************************************/

void CSkinColor::RunSGM(IplImage *img, IplImage **skinImg)

{

       if (img == NULL) return -1;

       //////////////////////////////////////////////////////////////////////////

       // 以下参数一组(117.4361,156.5599)来自源码 light2，与文章《王航宇:基于 YCbCr 高斯肤色模型的

       // 人脸检测技术研究》相同，另一组来自源码“肤色检测正式版”（103.0056, 140.1309）

       double M[]={103.0056, 140.1309}/*{117.4361,156.5599}*/;//M 为肤色在 YCbCr 颜色空间的样本均值(Cb, Cr)，经验值

       double C[2][2]={{160.1301,12.1430},//C 为肤色相似度模型的协方差矩阵，同上为经验值

              {12.1430,299.4574}};// 注：因为运算仅需要该矩阵的逆矩阵值，故该值没有使用，仅作参考

       double invC[2][2]={0.0077 ,-0.0041,-0.0041 ,0.0047

       };//Ct 为C的逆矩阵值，由matlab计算而得

       //////////////////////////////////////////////////////////////////////////

       IplImage* pImg = img;

       double CrMean=0,CbMean=0,YMean=0;

       // 1 颜色转换：BGR->YCrCb

       IplImage*imgYCrCb=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,3);// YCrCb图像

       cvCvtColor(pImg, imgYCrCb, CV_BGR2YCrCb);// 第0,1,2层分别为Y,Cr,Cb

      

       IplImage *imgY = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);// YCrCb图像

       IplImage *imgCr = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);// YCrCb图像

       IplImage *imgCb = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U,1);// YCrCb图像

       IplImage *imgY32 = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);// YCrCb图像

       IplImage *imgCr32 = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);// YCrCb图像

       IplImage *imgCb32 = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_32F,1);// YCrCb图像

       cvSplit(imgYCrCb, imgY, imgCr, imgCb, NULL);

       cvConvert(imgY, imgY32);

       cvConvert(imgCr, imgCr32);

       cvConvert(imgCb, imgCb32);

       //////////////////////////////////////////////////////////////////////////

       // 2 根据Sigle Gaussian Model计算颜色模型

       IplImage *PCbCr=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型

       IplImage *tempA=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型

       IplImage *tempB=cvCreateImage(cvGetSize(pImg), IPL_DEPTH_32F, 1);//YCrCb颜色模型

       cvSubS(imgCb32, cvScalar(M[0]), imgCb32);// x-m

       cvSubS(imgCr32, cvScalar(M[1]), imgCr32);// x-m

       cvAddWeighted(imgCb32, invC[0][0], imgCr32, invC[1][0], 0, tempA);

       cvAddWeighted(imgCb32, invC[0][1], imgCr32, invC[1][1], 0, tempB);

       cvMul(imgCb32, tempA, tempA, -0.5);

       cvMul(imgCr32, tempB, tempB, -0.5);

       cvAdd(tempA, tempB, PCbCr);

       cvExp(PCbCr, PCbCr);

       double max_val=0,min_val=0;

       cvMinMaxLoc(PCbCr,&min_val,&max_val);

       IplImage *proImg=cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U, 1);//YCrCb颜色模型

       double a=255/(max_val);

       cvConvertScaleAbs(PCbCr,proImg,a,0);

       m_proimg = cvCloneImage(proImg);

      

       if ((*skinImg)!=NULL) cvReleaseImage(skinImg);

       *skinImg = cvCreateImage(cvGetSize(pImg),IPL_DEPTH_8U, 1);//肤色结果

       // 释放内存

       cvReleaseImage(&proImg);

       cvReleaseImage(&imgYCrCb);

       cvReleaseImage(&imgY);

       cvReleaseImage(&imgCr);

       cvReleaseImage(&imgCb);

       cvReleaseImage(&imgY32);

       cvReleaseImage(&imgCr32);

       cvReleaseImage(&imgCb32);

       cvReleaseImage(&PCbCr);

       cvReleaseImage(&tempA);

       cvReleaseImage(&tempB);

}

 

3.1高斯混合模型

（1）以下matlab代码实现了高斯混合模型：

function [Alpha, Mu, Sigma] = GMM_EM(Data, Alpha0, Mu0, Sigma0)

%% EM 迭代停止条件

loglik_threshold = 1e-10;

%% 初始化参数

[dim, N] = size(Data);

M = size(Mu0,2);

loglik_old = -realmax;

nbStep = 0;

 

Mu = Mu0;

Sigma = Sigma0;

Alpha = Alpha0;

Epsilon = 0.0001;

while (nbStep < 1200)

  nbStep = nbStep+1;

  %% E-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % PDF of each point

    Pxi(:,i) = GaussPDF(Data, Mu(:,i), Sigma(:,:,i));         

  end

 

  % 计算后验概率 beta(i|x)

  Pix_tmp = repmat(Alpha,[N 1]).*Pxi;

  Pix = Pix_tmp ./ (repmat(sum(Pix_tmp,2),[1 M])+realmin);

  Beta = sum(Pix);

  %% M-步骤 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % 更新权值

    Alpha(i) = Beta(i) / N;

    % 更新均值

    Mu(:,i) = Data*Pix(:,i) / Beta(i);

    % 更新方差

    Data_tmp1 = Data - repmat(Mu(:,i),1,N);

    Sigma(:,:,i) = (repmat(Pix(:,i)',dim, 1) .* Data_tmp1*Data_tmp1') / Beta(i);

    %% Add a tiny variance to avoid numerical instability

    Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(dim,1));

  end

 

%  %% Stopping criterion 1 %%%%%%%%%%%%%%%%%%%%

%  for i=1:M

    %Compute the new probability p(x|i)

%    Pxi(:,i) = GaussPDF(Data, Mu(:,i), Sigma(i));

%  end

  %Compute the log likelihood

%  F = Pxi*Alpha';

%  F(find(F<realmin)) = realmin;

%  loglik = mean(log(F));

  %Stop the process depending on the increase of the log likelihood

%  if abs((loglik/loglik_old)-1) < loglik_threshold

%    break;

%  end

%  loglik_old = loglik;

 

  %% Stopping criterion 2 %%%%%%%%%%%%%%%%%%%%

  v = [sum(abs(Mu - Mu0)), abs(Alpha - Alpha0)];

  s = abs(Sigma-Sigma0);

  v2 = 0;

  for i=1:M

    v2 = v2 + det(s(:,:,i));

  end

 

  if ((sum(v) + v2) < Epsilon)

    break;

  end

  Mu0 = Mu;

  Sigma0 = Sigma;

  Alpha0 = Alpha;

end

nbStep

 

（2）以下代码根据高斯分布函数计算每组数据的概率密度，被GMM_EM函数所调用

function prob = GaussPDF(Data, Mu, Sigma)

%

% 根据高斯分布函数计算每组数据的概率密度 Probability Density Function (PDF)

% 输入 -----------------------------------------------------------------

%   o Data:  D x N ，N个D维数据

%   o Mu:    D x 1 ，M个Gauss模型的中心初始值

%   o Sigma: M x M ，每个Gauss模型的方差（假设每个方差矩阵都是对角阵，

%                                   即一个数和单位矩阵的乘积）

% Outputs ----------------------------------------------------------------

%   o prob:  1 x N array representing the probabilities for the

%            N datapoints.    

[dim,N] = size(Data);

Data = Data' - repmat(Mu',N,1);

prob = sum((Data*inv(Sigma)).*Data, 2);

prob = exp(-0.5*prob) / sqrt((2*pi)^dim * (abs(det(Sigma))+realmin));

 

（3）以下是演示代码demo1.m

% 高斯混合模型参数估计示例 （基于 EM 算法）

% 2010 年 11 月 9 日

[data, mu, var, weight] = CreateSample(M, dim, N);  // 生成测试数据

[Alpha, Mu, Sigma] = GMM_EM(Data, Priors, Mu, Sigma)

 

（4）以下是测试数据生成函数，为demo1.m所调用：

function [data, mu, var, weight] = CreateSample(M, dim, N)

% 生成实验样本集，由M组正态分布的数据构成

% % GMM模型的原理就是仅根据数据估计参数：每组正态分布的均值、方差，

% 以及每个正态分布函数在GMM的权重alpha。

% 在本函数中，这些参数均为随机生成，

%

% 输入

%   M    : 高斯函数个数

%   dim  : 数据维数

%   N    : 数据总个数

% 返回值

%   data : dim-by-N, 每列为一个数据

%   miu  : dim-by-M, 每组样本的均值，由本函数随机生成

%   var  : 1-by-M, 均方差，由本函数随机生成

%   weight: 1-by-M, 每组的权值，由本函数随机生成

% ----------------------------------------------------

%

% 随机生成不同组的方差、均值及权值

weight = rand(1,M);

weight = weight / norm(weight, 1); % 归一化，保证总合为1

var = double(mod(int16(rand(1,M)*100),10) + 1);  % 均方差，取1~10之间，采用对角矩阵

mu = double(round(randn(dim,M)*100));            % 均值，可以有负数

 

for(i = 1: M)

  if (i ~= M)

    n(i) = floor(N*weight(i));

  else

    n(i) = N - sum(n);

  end

end

 

% 以标准高斯分布生成样本值，并平移到各组相应均值和方差

start = 0;

for (i=1:M)

  X = randn(dim, n(i));

  X = X.* var(i) + repmat(mu(:,i),1,n(i));

  data(:,(start+1):start+n(i)) = X;

  start = start + n(i);

end

save('d:data.mat', 'data');

http://blog.sina.com.cn/s/blog_6a35405d01011s98.html

混合高斯模型的matlab代码剖析

什么是混合高斯模型？可以参见这篇博文，它同时给出了matlab实现代码，编者将数学表达式转化为矩阵运算的代码技巧很值得我们学习和借鉴，以下是个人对代码的剖析与理解。

 function varargout = gmm(X, K_or_centroids)
% ============================================================
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
%  - X: N-by-D data matrix.
%  - K_OR_CENTROIDS: either K indicating the number of
%       components or a K-by-D matrix indicating the
%       choosing of the initial K centroids.
%
%  - PX: N-by-K matrix indicating the probability of each
%       component generating each point.
%  - MODEL: a structure containing the parameters for a GMM:
%       MODEL.Miu: a K-by-D matrix.
%       MODEL.Sigma: a D-by-D-by-K matrix.
%       MODEL.Pi: a 1-by-K vector.
% ============================================================
 
    threshold = 1e-15;
    [N, D] = size(X);
 
    if isscalar(K_or_centroids)
        K = K_or_centroids;
        % randomly pick centroids
        rndp = randperm(N);
        centroids = X(rndp(1:K), :);
    else
        K = size(K_or_centroids, 1);
        centroids = K_or_centroids;
    end
 
    % initial values
    [pMiu pPi pSigma] = init_params(); %初始化
 
    Lprev = -inf; %inf表示正无究大，-inf表示为负无究大
    while true
        Px = calc_prob();
 
        % new value for pGamma
        pGamma = Px .* repmat(pPi, N, 1);
        pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K); %求每个样本由第K个聚类，也叫“component“生成的概率
 
        % new value for parameters of each Component
        Nk = sum(pGamma, 1);
        pMiu = diag(1./Nk) * pGamma' * X; %重新计算每个component的均值
        pPi = Nk/N; %更新混合高斯的加权系数
        for kk = 1:K %重新计算每个component的协方差
            Xshift = X-repmat(pMiu(kk, :), N, 1);
            pSigma(:, :, kk) = (Xshift' * ...
                (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);
        end
 
        % check for convergence
        L = sum(log(Px*pPi')); %求混合高斯分布的似然函数
        if L-Lprev < threshold %随着迭代次数的增加，似然函数越来越大，直至不变
            break; %似然函数收敛则退出
        end
        Lprev = L;
    end
 
    if nargout == 1 %如果返回是一个参数的话，那么varargout=Px;
        varargout = {Px};
    else %否则，返回[Px model],其中model是结构体
        model = [];
        model.Miu = pMiu;
        model.Sigma = pSigma;
        model.Pi = pPi;
        varargout = {Px, model};
    end
 
    function [pMiu pPi pSigma] = init_params()
        pMiu = centroids;
        pPi = zeros(1, K);
        pSigma = zeros(D, D, K);
 
        % hard assign x to each centroids
        distmat = repmat(sum(X.*X, 2), 1, K) + ... %distmat第j行的第i个元素表示第j个数据与第i个聚类点的距离，如果数据有4个，聚类2个，那么distmat就是4*2矩阵
            repmat(sum(pMiu.*pMiu, 2)', N, 1) - 2*X*pMiu'; %sum(A，2)结果为列向量，第i个元素是第i行的求和
        [dummy labels] = min(distmat, [], 2); %返回列向量dummy和labels，dummy向量记录distmat的每行的最小值，labels向量记录每行最小值的列号，即是第几个聚类，labels是N×1列向量，N为样本数
 
        for k=1:K
            Xk = X(labels == k, :); %把标志为同一个聚类的样本组合起来
            pPi(k) = size(Xk, 1)/N; %求混合高斯模型的加权系数，pPi为1*K的向量
            pSigma(:, :, k) = cov(Xk); %分别求单个高斯模型或聚类样本的协方差矩阵,pSigma为D*D*K的矩阵
        end
    end
 
    function Px = calc_prob()
        Px = zeros(N, K);
        for k = 1:K
            Xshift = X-repmat(pMiu(k, :), N, 1); %Xshift表示为样本矩阵-Uk，第i行表示xi-uk
            inv_pSigma = inv(pSigma(:, :, k)); %求协方差的逆
            tmp = sum((Xshift*inv_pSigma) .* Xshift, 2); %tmp为N*1矩阵，第i行表示（xi-uk)^T*Sigma^-1*(xi-uk)
            coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma)); %求多维正态分布中指数前面的系数
            Px(:, k) = coef * exp(-0.5*tmp); %求单独一个正态分布生成样本的概率或贡献
        end
    end
end

Free Mind

漫谈 Clustering (3): Gaussian Mixture Model

 http://blog.pluskid.org/?p=39

本文是“漫谈 Clustering 系列”中的第 4 篇，参见本系列的其他文章。

上一次我们谈到了用 k-means 进行聚类的方法，这次我们来说一下另一个很流行的算法：Gaussian Mixture Model (GMM)。事实上，GMM 和 k-means 很像，不过 GMM 是学习出一些概率密度函数来（所以 GMM 除了用在 clustering 上之外，还经常被用于 density estimation ），简单地说，k-means 的结果是每个数据点被 assign 到其中某一个 cluster 了，而 GMM 则给出这些数据点被 assign 到每个 cluster 的概率，又称作 soft assignment 。

得出一个概率有很多好处，因为它的信息量比简单的一个结果要多，比如，我可以把这个概率转换为一个 score ，表示算法对自己得出的这个结果的把握。也许我可以对同一个任务，用多个方法得到结果，最后选取“把握”最大的那个结果；另一个很常见的方法是在诸如疾病诊断之类的场所，机器对于那些很容易分辨的情况（患病或者不患病的概率很高）可以自动区分，而对于那种很难分辨的情况，比如，49% 的概率患病，51% 的概率正常，如果仅仅简单地使用 50% 的阈值将患者诊断为“正常”的话，风险是非常大的，因此，在机器对自己的结果把握很小的情况下，会“拒绝发表评论”，而把这个任务留给有经验的医生去解决。

废话说了一堆，不过，在回到 GMM 之前，我们再稍微扯几句。我们知道，不管是机器还是人，学习的过程都可以看作是一种“归纳”的过程，在归纳的时候你需要有一些假设的前提条件，例如，当你被告知水里游的那个家伙是鱼之后，你使用“在同样的地方生活的是同一种东西”这类似的假设，归纳出“在水里游的都是鱼”这样一个结论。当然这个过程是完全“本能”的，如果不仔细去想，你也不会了解自己是怎样“认识鱼”的。另一个值得注意的地方是这样的假设并不总是完全正确的，甚至可以说总是会有这样那样的缺陷的，因此你有可能会把虾、龟、甚至是潜水员当做鱼。也许你觉得可以通过修改前提假设来解决这个问题，例如，基于“生活在同样的地方并且穿着同样衣服的是同一种东西”这个假设，你得出结论：在水里有并且身上长有鳞片的是鱼。可是这样还是有问题，因为有些没有长鳞片的鱼现在又被你排除在外了。

在这个问题上，机器学习面临着和人一样的问题，在机器学习中，一个学习算法也会有一个前提假设，这里被称作“归纳偏执 (bias)”（bias 这个英文词在机器学习和统计里还有其他许多的意思）。例如线性回归，目的是要找一个函数尽可能好地拟合给定的数据点，它的归纳偏执就是“满足要求的函数必须是线性函数”。一个没有归纳偏执的学习算法从某种意义上来说毫无用处，就像一个完全没有归纳能力的人一样，在第一次看到鱼的时候有人告诉他那是鱼，下次看到另一条鱼了，他并不知道那也是鱼，因为两条鱼总有一些地方不一样的，或者就算是同一条鱼，在河里不同的地方看到，或者只是看到的时间不一样，也会被他认为是不同的，因为他无法归纳，无法提取主要矛盾、忽略次要因素，只好要求所有的条件都完全一样──然而哲学家已经告诉过我们了：世界上不会有任何样东西是完全一样的，所以这个人即使是有无比强悍的记忆力，也绝学不到任何一点知识。

这个问题在机器学习中称作“过拟合 (Overfitting)”，例如前面的回归的问题，如果去掉“线性函数”这个归纳偏执，因为对于 N 个点，我们总是可以构造一个 N-1 次多项式函数，让它完美地穿过所有的这 N 个点，或者如果我用任何大于 N-1 次的多项式函数的话，我甚至可以构造出无穷多个满足条件的函数出来。如果假定特定领域里的问题所给定的数据个数总是有个上限的话，我可以取一个足够大的 N ，从而得到一个（或者无穷多个）“超级函数”，能够 fit 这个领域内所有的问题。然而这个（或者这无穷多个）“超级函数”有用吗？只要我们注意到学习的目的（通常）不是解释现有的事物，而是从中归纳出知识，并能应用到新的事物上，结果就显而易见了。

没有归纳偏执或者归纳偏执太宽泛会导致 Overfitting ，然而另一个极端──限制过大的归纳偏执也是有问题的：如果数据本身并不是线性的，强行用线性函数去做回归通常并不能得到好结果。难点正在于在这之间寻找一个平衡点。不过人在这里相对于（现在的）机器来说有一个很大的优势：人通常不会孤立地用某一个独立的系统和模型去处理问题，一个人每天都会从各个来源获取大量的信息，并且通过各种手段进行整合处理，归纳所得的所有知识最终得以统一地存储起来，并能有机地组合起来去解决特定的问题。这里的“有机”这个词很有意思，搞理论的人总能提出各种各样的模型，并且这些模型都有严格的理论基础保证能达到期望的目的，然而绝大多数模型都会有那么一些“参数”（例如 K-means 中的 k ），通常没有理论来说明参数取哪个值更好，而模型实际的效果却通常和参数是否取到最优值有很大的关系，我觉得，在这里“有机”不妨看作是所有模型的参数已经自动地取到了最优值。另外，虽然进展不大，但是人们也一直都期望在计算机领域也建立起一个统一的知识系统（例如语意网就是这样一个尝试）。

废话终于说完了，回到 GMM 。按照我们前面的讨论，作为一个流行的算法，GMM 肯定有它自己的一个相当体面的归纳偏执了。其实它的假设非常简单，顾名思义，Gaussian Mixture Model ，就是假设数据服从 Mixture Gaussian Distribution ，换句话说，数据可以看作是从数个 Gaussian Distribution 中生成出来的。实际上，我们在 K-means 和 K-medoids 两篇文章中用到的那个例子就是由三个 Gaussian 分布从随机选取出来的。实际上，从中心极限定理可以看出，Gaussian 分布（也叫做正态 (Normal) 分布）这个假设其实是比较合理的，除此之外，Gaussian 分布在计算上也有一些很好的性质，所以，虽然我们可以用不同的分布来随意地构造 XX Mixture Model ，但是还是 GMM 最为流行。另外，Mixture Model 本身其实也是可以变得任意复杂的，通过增加 Model 的个数，我们可以任意地逼近任何连续的概率密分布。

每个 GMM 由 K 个 Gaussian 分布组成，每个 Gaussian 称为一个“Component”，这些 Component 线性加成在一起就组成了 GMM 的概率密度函数：

displaystyle egin{aligned} p(x) & = sum_{k=1}^K p(k)p(x|k) \ & = sum_{k=1}^K pi_k mathcal{N}(x|mu_k, Sigma_k) end{aligned}

根据上面的式子，如果我们要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K 个 Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数 pi_k ，选中了 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

那么如何用 GMM 来做 clustering 呢？其实很简单，现在我们有了数据，假定它们是由 GMM 生成出来的，那么我们只要根据数据推出 GMM 的概率分布来就可以了，然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算概率密度通常被称作 density estimation ，特别地，当我们在已知（或假定）了概率密度函数的形式，而要估计其中的参数的过程被称作“参数估计”。

现在假设我们有 N 个数据点，并假设它们服从某个分布（记作 p(x) ），现在要确定里面的一些参数的值，例如，在 GMM 中，我们就需要确定 pi_k、mu_k 和 Sigma_k 这些参数。 我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于 prod_{i=1}^N p(x_i) ，我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和 sum_{i=1}^N log p(x_i)，得到 log-likelihood function 。接下来我们只要将这个函数最大化（通常的做法是求导并令导数等于零，然后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，我们就认为这是最合适的参数，这样就完成了参数估计的过程。

下面让我们来看一看 GMM 的 log-likelihood function ：

displaystyle sum_{i=1}^N log left{sum_{k=1}^K pi_k mathcal{N}(x_i|mu_k, Sigma_k) ight}

由于在对数函数里面又有加和，我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们采取之前从 GMM 中随机选点的办法：分成两步，实际上也就类似于 K-means的两步。

  1. 估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率）：对于每个数据 x_i 来说，它由第 k 个 Component 生成的概率为displaystyle gamma(i, k) = frac{pi_k mathcal{N}(x_i|mu_k, Sigma_k)}{sum_{j=1}^K pi_jmathcal{N}(x_i|mu_j, Sigma_j)}

    由于式子里的 mu_k 和 Sigma_k 也是需要我们估计的值，我们采用迭代法，在计算 gamma(i, k) 的时候我们假定 mu_k 和 Sigma_k 均已知，我们将取上一次迭代所得的值（或者初始值）。

  2. 估计每个 Component 的参数：现在我们假设上一步中得到的 gamma(i, k) 就是正确的“数据 x_i 由 Component k 生成的概率”，亦可以当做该 Component 在生成这个数据上所做的贡献，或者说，我们可以看作 x_i 这个值其中有 gamma(i, k)x_i 这部分是由 Component k 所生成的。集中考虑所有的数据点，现在实际上可以看作 Component 生成了 gamma(1, k)x_1, ldots, gamma(N, k)x_N 这些点。由于每个 Component 都是一个标准的 Gaussian 分布，可以很容易分布求出最大似然所对应的参数值：displaystyle egin{aligned} mu_k & = frac{1}{N_k}sum_{i=1}^Ngamma(i, k)x_i \ Sigma_k & = frac{1}{N_k}sum_{i=1}^Ngamma(i, k)(x_i-mu_k)(x_i-mu_k)^T end{aligned}

    其中 N_k = sum_{i=1}^N gamma(i, k) ，并且 pi_k 也顺理成章地可以估计为 N_k/N 。

  3. 重复迭代前面两步，直到似然函数的值收敛为止。

当然，上面给出的只是比较“直观”的解释，想看严格的推到过程的话，可以参考 Pattern Recognition and Machine Learning 这本书的第九章。有了实际的步骤，再实现起来就很简单了。Matlab 代码如下：

（Update 2012.07.03：如果你直接把下面的代码拿去运行了，碰到 covariance 矩阵 singular 的情况，可以参见这篇文章。）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

function varargout = gmm(X, K_or_centroids)
% ============================================================
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
%  - X: N-by-D data matrix.
%  - K_OR_CENTROIDS: either K indicating the number of
%       components or a K-by-D matrix indicating the
%       choosing of the initial K centroids.
%
%  - PX: N-by-K matrix indicating the probability of each
%       component generating each point.
%  - MODEL: a structure containing the parameters for a GMM:
%       MODEL.Miu: a K-by-D matrix.
%       MODEL.Sigma: a D-by-D-by-K matrix.
%       MODEL.Pi: a 1-by-K vector.
% ============================================================
 
    threshold = 1e-15;
    [N, D] = size(X);
 
    if isscalar(K_or_centroids)
        K = K_or_centroids;
        % randomly pick centroids
        rndp = randperm(N);
        centroids = X(rndp(1:K), :);
    else
        K = size(K_or_centroids, 1);
        centroids = K_or_centroids;
    end
 
    % initial values
    [pMiu pPi pSigma] = init_params();
 
    Lprev = -inf;
    while true
        Px = calc_prob();
 
        % new value for pGamma
        pGamma = Px .* repmat(pPi, N, 1);
        pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);
 
        % new value for parameters of each Component
        Nk = sum(pGamma, 1);
        pMiu = diag(1./Nk) * pGamma' * X;
        pPi = Nk/N;
        for kk = 1:K
            Xshift = X-repmat(pMiu(kk, :), N, 1);
            pSigma(:, :, kk) = (Xshift' * ...
                (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);
        end
 
        % check for convergence
        L = sum(log(Px*pPi'));
        if L-Lprev < threshold
            break;
        end
        Lprev = L;
    end
 
    if nargout == 1
        varargout = {Px};
    else
        model = [];
        model.Miu = pMiu;
        model.Sigma = pSigma;
        model.Pi = pPi;
        varargout = {Px, model};
    end
 
    function [pMiu pPi pSigma] = init_params()
        pMiu = centroids;
        pPi = zeros(1, K);
        pSigma = zeros(D, D, K);
 
        % hard assign x to each centroids
        distmat = repmat(sum(X.*X, 2), 1, K) + ...
            repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...
            2*X*pMiu';
        [dummy labels] = min(distmat, [], 2);
 
        for k=1:K
            Xk = X(labels == k, :);
            pPi(k) = size(Xk, 1)/N;
            pSigma(:, :, k) = cov(Xk);
        end
    end
 
    function Px = calc_prob()
        Px = zeros(N, K);
        for k = 1:K
            Xshift = X-repmat(pMiu(k, :), N, 1);
            inv_pSigma = inv(pSigma(:, :, k));
            tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
            coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));
            Px(:, k) = coef * exp(-0.5*tmp);
        end
    end
end

函数返回的 Px 是一个 N imes K 的矩阵，对于每一个 x_i ，我们只要取该矩阵第 i 行中最大的那个概率值所对应的那个 Component 为 x_i 所属的 cluster 就可以实现一个完整的聚类方法了。对于最开始的那个例子，GMM 给出的结果如下：

gmm

相对于之前 K-means 给出的结果，这里的结果更好一些，左下角的比较稀疏的那个 cluster 有一些点跑得比较远了。当然，因为这个问题原本就是完全有 Mixture Gaussian Distribution 生成的数据，GMM （如果能求得全局最优解的话）显然是可以对这个问题做到的最好的建模。

另外，从上面的分析中我们可以看到 GMM 和 K-means 的迭代求解法其实非常相似（都可以追溯到 EM 算法，下一次会详细介绍），因此也有和 K-means 同样的问题──并不能保证总是能取到全局最优，如果运气比较差，取到不好的初始值，就有可能得到很差的结果。对于 K-means 的情况，我们通常是重复一定次数然后取最好的结果，不过 GMM 每一次迭代的计算量比 K-means 要大许多，一个更流行的做法是先用 K-means （已经重复并取最优值了）得到一个粗略的结果，然后将其作为初值（只要将 K-means 所得的 centroids 传入 gmm 函数即可），再用 GMM 进行细致迭代。

如我们最开始所讨论的，GMM 所得的结果（Px）不仅仅是数据点的 label ，而包含了数据点标记为每个 label 的概率，很多时候这实际上是非常有用的信息。最后，需要指出的是，GMM 本身只是一个模型，我们这里给出的迭代的办法并不是唯一的求解方法。感兴趣的同学可以自行查找相关资料。

Free Mind

Regularized Gaussian Covariance Estimation

http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/

我之前写过一篇介绍 Gaussian Mixture Model (GMM) 的文章，并在文章里贴了一段 GMM 实现的 Matlab 示例代码，然后就不断地有人来问我关于那段代码的问题，问得最多的就是大家经常发现在跑那段代码的时候估计出来的 Covariance Matrix 是 singular 的，所以在第 96 行求逆的时候会挂掉。这是今天要介绍的主要话题，我会讲得罗嗦一点，把关于那篇文章我被问到的一些其他问题也都提到一下，不过，在步入正题之前，不妨先来小小地闲扯一下。

我自己以前就非常在痛恨书里看到伪代码，又不能编译运行，还搞得像模像样的像代码一般，所以我自己在写 blog 的时候就尝试给出直接可运行的代码，但是现在我渐渐理解为什么书的作者喜欢用伪代码而不是可运行的代码了（除非是专门讲某编程语言的书）。原因很简单：示例用的代码和真实项目中的代码往往是差别很大的，直接给出可运行的示例代码，读者就会直接拿去运行（因为包括我在内的大部分人都是很懒的，不是说“懒”是程序员的天性么？），而往往示例代码为了结构清晰并突出算法的主要脉络，对很多琐碎情况都没有处理，都使用最直接的操作，没有做任何优化（例如我的那个 GMM 示例代码里直接用 matlab 的 inv 函数来求 Covariance 矩阵的逆），等等，因此直接运行这样的代码——特别是在实际项目中使用，是非常不明智的。

可是那又为什么不直接给出实际项目级别的代码呢？第一个原因当然是工作量，我想程序员都知道从一个算法到一个健壮、高效的系统之间有多长的路要走，而且很多时候都需要根据不同的项目环境和需要做不同的修改。所以当一个人实际上是在介绍算法的时候让他去弄一套工业级别的代码出来作为示例，实在是太费力了。当然更重要的原因是：工业级别的代码通常经过大量的优化并充斥着大量错误处理、以及为了处理其他地方的 bug 或者整个系统混乱的 API 等而做的各种 workaround 等等……也许根本就看不懂，基本完全失去了示例的作用。所以，伪代码就成为了唯一的选择。

罗嗦了半天，现在我们回到正题，那篇文章讲的是 GMM 的学习，所以关于 GMM 的问题可以直接参考那篇文章，出问题的地方仅在于 GMM 的每个 component (Gaussian Distribution) 的参数估计上，特别地，在于 Covariance Matrix 的估计上。所以，我们今天主要来探讨 Gaussian 分布的协方差矩阵的估计问题。首先来看一下多维 Gaussian 分布的概率密度函数：

(eq: 1)

参数估计的过程是给定一堆数据点 ，并假设他们是由某个真实的 Gaussian 分布独立地 (IID) 采样出来的，然后根据这堆点来估计该 Gaussian 分布的参数。对于一个多维 Gaussian 分布来说，需要估计的就是均值  和协方差矩阵  两“个”参数——用引号是因为  和  都不是标量，假设我们的数据是  维的，那么  实际上是  个参数，而作为一个  的矩阵， 则是由  个参数组成。根据协方差矩阵的对称性，我们可以进一步把  所包含的参数个数降低为  个。

值得一提的是，我们从参数估计的角度并不能一下子断定  就是对称的，需要稍加说明（这实际上是《Pattern Recognition and Machine Learning》书里的习题 2.17）。记 ，假设我们的  不是对称的，那么我们把他写成一个对称矩阵和一个反对称阵之和（任何一个矩阵都可以这样表示）：

将起代入式 (eq: 1) 的指数部分中，注意到对于反对称部分（书写方便起见我们暂时另 ）：

所以指数部分只剩下对称阵  了。也就是说，如果我们找出一个符合条件的非对称阵 ，用它的对称部分  代替也是可以的，所以我们就可以“不失一般性”地假设  是对称的（因为对称阵的逆矩阵也是对称的）。注意我们不用考虑式 (eq: 1) 指数前面的系数（里的那个  的行列式），因为前面的系数是起 normalization 作用的，使得全概率积分等于 1，指数部分变动之后只要重新 normalize 一下即可（可以想像一下，如果  真的是非对称的，那么 normalization 系数那里应该也不会是  那么简单的形式吧？有兴趣的同学可以自己研究一下）。

好，现在回到参数估计上。对于给定的数据 ，将其代入式 (eq: 1) 中可以得到其对应的 likelihood ，注意这个并不是  的概率。和离散概率分布不同的是，这里谈论单个数据点的概率是没有多大意义的，因为单点集的测度为零，它的概率也总是为零。所以，这里的 likelihood 出现大于 1 的值也是很正常的——这也是我经常被问到的一个问题。虽然不是概率值，但是我想从 likelihood 或者概率密度的角度去理解，应该也是可以得到直观的认识的。

参数估计常用的方法是最大似然估计（Maximum Likelihood Estimation, MLE）。就是将所有给定数据点的似然全部乘起来，得到一个关于参数（这里是  和 ）的函数，然后最大化该函数，所对应的参数值就是最大似然估计的值。通常为了避免数值下溢，会使用单调函数  将相乘变成相加在计算。

对于最大似然估计，我的理解是寻找最可能生成给定的数据的模型参数。这对于离散型的分布来说是比较好解释的，因为此时一点的似然实际上就是该点的概率，不过对于像 Gaussian 这样的针对连续数据的分布来说，就不是那么自然了。不过我们其实有另一种解释，这是我在 Coursera 的 Probabilistic Graphical Model 课上看到的：

对于一个真实的分布 ，假设我们得到了它的一个估计 ，那么如何来衡量这个估计的好坏呢？比较两个概率分布的一个常用的办法是用 KL-Divergence，又叫 Relative Entropy（以及其他若干外号），定义为：

不过由于我们通常是不知道真实的  的，所以这个值也没法直接算，于是我们把它简单地变形一下：

其中蓝色部分是分布  自己的 Entropy，由于我们比较不同的模型估计好坏的时候这个量是不变的，所以我们可以忽略它，只考虑红色部分。为了使得 KL-Divergence 最小，我们需要将红色部分（注意不包括前面的负号）最大化。不过由于仍然依赖于真实分布 ，所以红色部分还是不能直接计算，但是这是一个我们熟悉的形式：某个函数关于真实分布的期望值，而我们手里面又有从真实分布里 IID 采样出来的  个数据点，所以可以直接使用标准的近似方法：用数据的经验分布期望来代替原始的期望：

于是我们很自然地就得到了最大似然估计的目标函数。这也从另一个角度解释了最大似然估计的合理性：同时对离散型和连续型的数据都适用。

接下来我们再回到 Gaussian 分布，将 Gaussian 的概率密度函数 (eq: 1) 代入最大似然估计的目标函数，可以得到：

(eq: 2)

由于本文主要讨论协方差矩阵的估计，我们这里就不管均值矩阵  的估计了，方法上是类似的，下面式子中的  将表示已知的或者同样通过最大似然估计计算出来的均值。为了找到使得上式最大化的协方差矩阵，我们将它对于  求导：

令导数等于零，即可得到最优解：

(eq: 3)

于是  就是协方差矩阵的最大似然估计。到这里问题就出来了，再回到 Gaussian 分布的概率密度函数 (eq: 1)，如果把 带进去的话，可以看到是需要对它进行求逆的，这也是我被问得最多的问题：由于  singular 了，导致求逆失败，于是代码也就运行失败了。

从 (eq: 3) 可以大致看到什么时候  会 singular：由于该式的形式，可以知道  的 rank 最大不会超过 ，所以如果 ，也就是数据的维度大于数据点的个数的话，得到的估计值  肯定是不满秩的，于是就 singular 了。当然，即使不是这么极端的情况，如果  比较小或者  比较大的时候，仍然会有很大的几率出现不满秩的。

这个问题还可以从更抽象一点的角度来看。我们之前计算过估计 Covariance 矩阵实际上是要估计  这么多个参数，如果  比较大的话，这个数目也会变得非常多，也就意味着模型的复杂度很大，对于很复杂的模型，如果训练数据的量不够多，则无法确定一个唯一的最优模型出来，这是机器学习中非常常见的一个问题，通常称为overfitting。解决 overfitting 问题的方法有不少，最直接的就是用更多的训练数据（也就是增大 ），当然，有时候训练数据只有那么多，于是就只好从反面入手，降低模型复杂度。

在降低模型复杂度方面最常用的方法大概是正则化了，而正则化最简单的例子也许是 Ridge Regression 了，通过在目标函数后面添加一项关于参数的 -norm 的项来对模型参数进行限制；此外也有诸如 LASSO 之类的使用 -norm 来实现稀疏性的正则化，这个我们在之前也曾经简单介绍过。

在介绍估计 Covariance 的正则化方法之前，我们首先来看一种更直接的方法，直接限制模型的复杂度：特别地，对于我们这里 Gaussian 的例子，比如，我们可以要求协方差矩阵是一个对角阵，这样一来，对于协方差矩阵，我们需要估计的参数个数就变成了  个对角元素，而不是原来的  个。大大减少参数个数的情况下，overfitting 的可能性也极大地降低了。

注意 Covariance Matrix 的非对角线上的元素是对应了某两个维度之间的相关性 (Correlation) 的，为零就表示不相关。而对于 Gaussian 分布来说，不相关和独立 (Independence) 是等价的，所以说在我们的假设下，向量  的各个维度之间是相互独立的，换句话说，他们的联合概率密度可以分解为每个维度各自的概率密度（仍然是 Gaussian 分布）的乘积，这个特性使得我们可以很方便地对协方差矩阵进行估计：只有每个维度上当做一维的 Gaussian 分布各自独立地进行参数估计就可以了。注意这个时候除非某个维度上的坐标全部是相等的，否则 Covariance 矩阵对角线上的元素都能保证大于零，也就不会出现 singular 的情况了。

那么直观上将协方差矩阵限制为对角阵会产生什么样的模型呢？我们不妨看一个二维的简单例子。下面的代码利用了scikit-learn 里的 GMM 模块来估计一个单个 component 的 Gaussian，它在选项里已经提供了 Covariance 的限制，这里我们除了原始的 full 和对角阵的 diag 之外，还给出一个 spherical，它假定协方差矩阵是  这样的形式，也就是说，不仅是对角阵，而且所有对角元素都相等。下面的代码把三种情况下 fit 出来的模型画出来，在图 (fig: 1) 中可以看到。

  1. import numpy as np
  2. import pylab as pl
  3. from sklearn import mixture
  5. n_samples = 300
  6. c_types = ['full', 'diag', 'spherical']
  8. np.random.seed(0)
  9. C = np.array([[0., -0.7], [3.5, 1.7]])
  10. X_train = np.dot(np.random.randn(n_samples, 2), C)
  12. pl.figure(dpi=100,figsize=(3,3))
  13. pl.scatter(X_train[:, 0], X_train[:, 1], .8)
  14. pl.axis('tight')
  15. pl.savefig('GaussianFit-data.svg')
  16. pl.close()
  18. for c_type in c_types:
  19. clf = mixture.GMM(n_components=1, covariance_type=c_type)
  20. clf.fit(X_train)
  22. x = np.linspace(-15.0, 20.0, num=200)
  23. y = np.linspace(-10.0, 10.0, num=200)
  24. X, Y = np.meshgrid(x, y)
  25. XX = np.c_[X.ravel(), Y.ravel()]
  26. Z = np.log(-clf.eval(XX)[0])
  27. Z = Z.reshape(X.shape)
  29. pl.figure(dpi=100,figsize=(3,3))
  30. CS = pl.contour(X, Y, Z)
  31. pl.scatter(X_train[:, 0], X_train[:, 1], .8)
  33. pl.axis('tight')
  34. pl.savefig('GaussianFit-%s.svg' % c_type)
  35. pl.close()
图 1
Training samples.
Spherical Covariance.
Diagonal Covariance.
Full Covariance.

可以看到，从直观上来讲：spherical Covariance 的时候就是各个维度上的方差都是一样的，所以形成的等高线在二维上是正圆。而 diagonal Covariance 则可以 capture 更多的特征，因为它允许形成椭圆的等高线，但是它的限制是这个椭圆是不能旋转的，永远是正放的。Full Covariance 是 fit 最好的，但是前提是数据量足够来估计这么多的参数，虽然对于 2 维这样的低维度通常并不是问题，但是维度高了之后就经常会变得很困难。

一般来说，使用对角线的协方差矩阵是一种不错的降低模型复杂度的方式，在 scikit-learn 的 GMM 模块中这甚至是默认的。不过如果有时候你觉得这样的限制实在是和你的真实数据想去甚远的话，也有另外的处理方式。最简单的处理办法是在求出了 Covariance 矩阵的估计值之后直接加上一个 ，也就是在对角线上统一加上一个很小的正数 ，通常称作正则化系数。例如，Matlab 自带的统计工具箱里的 gmdistribution.fit 函数就有一个可选参数叫做Regularize，就是我们这里说的 。为什么加上  之后就不会 singular 了呢？首先  本身至少肯定是半正定的，所以 ：

所以这样之后肯定就是正定的了。此外在 scikit-learn 的 GMM 学习的函数中也有一个类似的参数 min_covar，它的文档里说的是“Floor on the diagonal of the covariance matrix to prevent overfitting. Defaults to 1e-3.”感觉好像是如果原来的 Covariance 矩阵对角线元素小于该参数的时候就将它设置为该参数，不过这样的做法是不是一定产生正定的协方差矩阵似乎就不像上面那种那样直观可以得出结论了。

不过这样的做法虽然能够解决 singular 的问题，但是总让人有些难以信服：你莫名其妙地在我的协方差矩阵上加上了一些东西，这到底是什么意思呢？最简单的解释可以从式 (eq: 1) 那里协方差矩阵你的地方来看，对于矩阵求逆或者解方程的时候出现 singular 的情况，加上一个  也算是数值上的一种标准处理方式，叫做 Tikhonov Regularization，Ridge Regression 本身也可以通过这个角度来解释。从数值计算的角度来说，这样处理能够增加数值稳定性，直观地来讲，稳定性是指  元素的微小数值变化，不会使得求逆（或者解方程之类的）之后的解产生巨大的数值变化，如果 是 singular 的话，通常就不具有这样的稳定性了。此外，随着 regularization coefficient  逐渐趋向于零，对应的解也会逐渐趋向于真实的解。

如果这个解释还不能令出信服的话，我们还可以从 Bayes 推断的角度来看这个问题。不过这本身是一个很大的话题，我在这里只能简略地讲一个思路，想了解更多的同学可以参考《Pattern Recognition and Machine Learning》第二章或者其他相关的书籍。

总的来说，我们这里要通过 Maximum a posteriori (MAP) Estimation 来对协方差矩阵进行估计。做法是首先为  定义一个先验分布 (prior) ，然后对于数据 ，我们根据 Bayes 公式，有

(eq: 4)

其中等式左边称作后验 (posterior)，右边红色部分是似然 (likelihood)，蓝色部分是用于将概率分布 normalize 到全概率为 1 的，通常并不需要直接通过  去计算，而只要对求出的概率分布进行归一化就可以了。MAP 的思路是将 posterior 最大的那个  作为最佳估计值，直观上很好理解，就是给定了数据点之后“最可能的”那个 。而计算上我们是通过 (eq: 4) 来实现的，因为分母是与  无关的，所以在计算的时候（因为只要比较相对大小嘛）忽略掉，这样一来，就可以看到，其实 MAP 方法和刚才我们用的最大似然 (MLE) 方法唯一的区别就在于前面乘了一个先验分布。

虽然从最终的计算公式上来看差别很细微，但是 MAP 却有很多好处。和本文最相关的当然是先验的引入，这在很大程度上和我们平时用正则化的目的是一样的：加入我们自己对模型的先验知识或者假设或者要求，这样可以避免在数据不够的时候产生 overfitting 。此外还有另一个值得一提的好玩特性就是 MAP 可以在线计算：注意到后验分布本身也是关于  的一个合法分布，所以在我们计算出后验之后，如果又得到了新的训练数据，那么我们可以将之前算出来的后验又作为先验，代入新的一轮的计算，这样可以一直不断地重复下去。:D

不过，虽然 MAP 框架看上去很美，但是在 general 的情况下计算却可能会变得很复杂，因此，虽然理论上来说，我们可以把任何先验知识通过先验分布的方式加入到模型中来，但是在实际中先验的选择往往是根据“哪个更好计算”而不是根据“哪个更合理”的准则来做出的。基于这个准则，有了 Conjugate Prior 的概念，简单地来说，如果一个 prior 和 likelihood 相乘再归一化之后得到的 posterior 形式上是和 prior 一样的话，那么它就是该 likelihood 的一个 conjugate prior。

选择 conjugate prior 的好处是显而易见的，由于先验和后验的形式是一样的，都是某一类型的分布，只是参数的值发生了变化，所以说可以看成是在观察到数据之后对模型参数进行修正（注意这里的模型是关于  的，而  本身又是关于数据的模型——一个 Gaussian 分布——的参数，不要搞混淆了），并且这种修正有时候可以得到非常直观的解释，不过我在这里就不细说了。此外，由于形式没有发生变化，所以在使用在线计算的时候，每一次观察到新的数据所使用的计算公式都还是一样的。

回到我们的问题，关于 Gaussian 分布的协方差矩阵的 Conjugate Prior 是一个叫做 Inverse Wishart Distribution 的家伙，它的概率密度函数是这个样子的：

其中  和参数  都是  的正定矩阵，而  则是 Multivariate Gamma Function。是不是看上去有点恐怖？^_^bbbb我也觉得有点恐怖……接下来让我们来看 MAP，首先将这个 prior 乘以我们的 likelihood，和 MLE 一样的，我们在  空间中计算，所以对整个式子取对数，然后丢掉和  无关的常数项（因为是关于  最大化嘛），得到下面的式子：

可以看到现在已经变得不是那么复杂了，形式上和我们之前的 (eq: 2) 是一样的，只是多了红色的一项，接下来对  进行求导并令导数等于零，类似地得到如下方程：

于是得到最优解为

这里  是之前的最大似然估计。接下来如果我们将我们 prior 中的参数  选为 ，就可以得到 MAP 估计为：

就得到了我们刚才的正则化的形式，虽然有一点细微的区别就是这里前面多了一个缩放系数。这样一来，我们就得到了这种正则化方法的一个看起来比较靠谱的解释：实际上就是我们先验地假设协方差矩阵满足一个均值为 spherical 形状（所有对角元素都相等的一个对角阵）的 Inverse Wishart 分布（关于 Inverse Wishart 分布的均值的公式可以参见wikipedia），然后通过数据来对这个先验进行修正得到后验概率，并从后验中取了概率（密度）最大的那个  作为最终的估计。

这种方式比暴力地直接强制要求协方差矩阵具有 diagonal 或者 spherical 形式要温柔得多，并且随着训练数据的增多，先验的影响也会逐渐减小，从而趋向于得到真实的参数，而强制结构限制的结果则是不论你最后通过什么手段搞到了多少数据，如果真实模型本身不是 diagonal 或者 spherical 的话，那么你的参数估计的 bias 将永远都无法被消除，这个我们已经在图 (fig: 1) 中直观地看到了。当然，在训练数据本身就不足的情况下，两种选择谁好谁好就没有定论了，不要忘了，强制结构限制的计算复杂度要小很多。:)

最后，对于不想看长篇大论的人，我做一下简单的总结：如果碰到估计 Gaussian 的协方差矩阵时结果 singular 的情况，解决方法一般有两种：

  1. 直接限制协方差矩阵的结构，比如要求它是一个对角阵，这种情况下数据的各个维度将会是独立的，所以可以把每个维度看成一个一维的 Gaussian 分布来单独进行估计。
  2. 在估计出来的协方差矩阵的对角线上统一加上一个很小的正数，使得它变成正定的。这种方法的合理性可以通过正则化或者 MAP 估计来进行解释。

然后再补充几点不知道该放在哪里的注意事项：

  • prior、posterior 和 likelihood 比较容易搞混淆，记住 prior 和 posterior 都是关于参数（比如我们这里就是 ）的分布，prior 是我们所假设的参数本来的分布，而 posterior 则是在观察到训练数据之后得到的条件分布。而 likelihood 则完全不一样，一方面它是关于数据  的，另一方面它没有被归一化，所以也并不是一个合法的概率分布。
  • 标量值函数关于矩阵变量的求导原理上其实就是和多元标量值函数求导一样的，不过求起来非常繁琐，一般不会自己算，网上可以找到一个叫做《Matrix Cookbook》的小册子，里面有搜集了各种常用的矩阵求导的公式，基本上直接从那里查询就可以解决大部分问题了。

GMM的EM算法实现

http://blog.csdn.net/abcjennifer/article/details/8198352

在 聚类算法K-Means, K-Medoids, GMM, Spectral clustering，Ncut一文中我们给出了GMM算法的基本模型与似然函数，在EM算法原理中对EM算法的实现与收敛性证明进行了详细说明。本文主要针对如何用EM算法在混合高斯模型下进行聚类进行代码上的分析说明。

1. GMM模型：

每个 GMM 由 K 个 Gaussian 分布组成，每个 Gaussian 称为一个“Component”，这些 Component 线性加成在一起就组成了 GMM 的概率密度函数：

根据上面的式子，如果我们要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K个Gaussian Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数 pi(k) ，选中了 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

那么如何用 GMM 来做 clustering 呢？其实很简单，现在我们有了数据，假定它们是由 GMM 生成出来的，那么我们只要根据数据推出 GMM 的概率分布来就可以了，然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算概率密度通常被称作 density estimation ，特别地，当我们在已知（或假定）了概率密度函数的形式，而要估计其中的参数的过程被称作“参数估计”。

2. 参数与似然函数：

现在假设我们有 N 个数据点，并假设它们服从某个分布（记作 p(x) ），现在要确定里面的一些参数的值，例如，在 GMM 中，我们就需要确定 影响因子pi(k)、各类均值pMiu(k) 和 各类协方差pSigma(k) 这些参数。 我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于  ，我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和 sum_{i=1}^N log p(x_i)，得到 log-likelihood function 。接下来我们只要将这个函数最大化（通常的做法是求导并令导数等于零，然后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，我们就认为这是最合适的参数，这样就完成了参数估计的过程。

下面让我们来看一看 GMM 的 log-likelihood function ：

由于在对数函数里面又有加和，我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们采取之前从 GMM 中随机选点的办法：分成两步，实际上也就类似于K-means 的两步。

3. 算法流程：

1.  估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率）：对于每个数据 x_i 来说，它由第 k 个 Component 生成的概率为

其中N（xi | μk,Σk）就是后验概率。

2. 通过极大似然估计可以通过求到令参数=0得到参数pMiu，pSigma的值。具体请见这篇文章第三部分。

其中 N_k = sum_{i=1}^N gamma(i, k) ，并且 pi_k 也顺理成章地可以估计为 N_k/N 。

3. 重复迭代前面两步，直到似然函数的值收敛为止。

4. matlab实现GMM聚类代码与解释：

 

说明：fea为训练样本数据，gnd为样本标号。算法中的思想和上面写的一模一样，在最后的判断accuracy方面，由于聚类和分类不同，只是得到一些 cluster ，而并不知道这些 cluster 应该被打上什么标签，或者说。由于我们的目的是衡量聚类算法的 performance ，因此直接假定这一步能实现最优的对应关系，将每个 cluster 对应到一类上去。一种办法是枚举所有可能的情况并选出最优解，另外，对于这样的问题，我们还可以用 Hungarian algorithm 来求解。具体的Hungarian代码我放在了资源里，调用方法已经写在下面函数中了。

注意：资源里我放的是Kmeans的代码，大家下载的时候只要用bestMap.m等几个文件就好~

 

1. gmm.m，最核心的函数，进行模型与参数确定。

[cpp] view plain copy
 
  1. function varargout = gmm(X, K_or_centroids)  
  2. % ============================================================  
  3. % Expectation-Maximization iteration implementation of  
  4. % Gaussian Mixture Model.  
  5. %  
  6. % PX = GMM(X, K_OR_CENTROIDS)  
  7. % [PX MODEL] = GMM(X, K_OR_CENTROIDS)  
  8. %  
  9. %  - X: N-by-D data matrix.  
  10. %  - K_OR_CENTROIDS: either K indicating the number of  
  11. %       components or a K-by-D matrix indicating the  
  12. %       choosing of the initial K centroids.  
  13. %  
  14. %  - PX: N-by-K matrix indicating the probability of each  
  15. %       component generating each point.  
  16. %  - MODEL: a structure containing the parameters for a GMM:  
  17. %       MODEL.Miu: a K-by-D matrix.  
  18. %       MODEL.Sigma: a D-by-D-by-K matrix.  
  19. %       MODEL.Pi: a 1-by-K vector.  
  20. % ============================================================  
  21. % @SourceCode Author: Pluskid (http://blog.pluskid.org)  
  22. % @Appended by : Sophia_qing (http://blog.csdn.net/abcjennifer)  
  23.       
  24.   
  25. %% Generate Initial Centroids  
  26.     threshold = 1e-15;  
  27.     [N, D] = size(X);  
  28.    
  29.     if isscalar(K_or_centroids) %if K_or_centroid is a 1*1 number  
  30.         K = K_or_centroids;  
  31.         Rn_index = randperm(N); %random index N samples  
  32.         centroids = X(Rn_index(1:K), :); %generate K random centroid  
  33.     else % K_or_centroid is a initial K centroid  
  34.         K = size(K_or_centroids, 1);   
  35.         centroids = K_or_centroids;  
  36.     end  
  37.    
  38.     %% initial values  
  39.     [pMiu pPi pSigma] = init_params();  
  40.    
  41.     Lprev = -inf; %上一次聚类的误差  
  42.       
  43.     %% EM Algorithm  
  44.     while true  
  45.         %% Estimation Step  
  46.         Px = calc_prob();  
  47.    
  48.         % new value for pGamma(N*k), pGamma(i,k) = Xi由第k个Gaussian生成的概率  
  49.         % 或者说xi中有pGamma(i,k)是由第k个Gaussian生成的  
  50.         pGamma = Px .* repmat(pPi, N, 1); %分子 = pi(k) * N(xi | pMiu(k), pSigma(k))  
  51.         pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K); %分母 = pi(j) * N(xi | pMiu(j), pSigma(j))对所有j求和  
  52.    
  53.         %% Maximization Step - through Maximize likelihood Estimation  
  54.           
  55.         Nk = sum(pGamma, 1); %Nk(1*k) = 第k个高斯生成每个样本的概率的和，所有Nk的总和为N。  
  56.           
  57.         % update pMiu  
  58.         pMiu = diag(1./Nk) * pGamma' * X; %update pMiu through MLE(通过令导数 = 0得到)  
  59.         pPi = Nk/N;  
  60.           
  61.         % update k个 pSigma  
  62.         for kk = 1:K   
  63.             Xshift = X-repmat(pMiu(kk, :), N, 1);  
  64.             pSigma(:, :, kk) = (Xshift' * ...  
  65.                 (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);  
  66.         end  
  67.    
  68.         % check for convergence  
  69.         L = sum(log(Px*pPi'));  
  70.         if L-Lprev < threshold  
  71.             break;  
  72.         end  
  73.         Lprev = L;  
  74.     end  
  75.    
  76.     if nargout == 1  
  77.         varargout = {Px};  
  78.     else  
  79.         model = [];  
  80.         model.Miu = pMiu;  
  81.         model.Sigma = pSigma;  
  82.         model.Pi = pPi;  
  83.         varargout = {Px, model};  
  84.     end  
  85.    
  86.     %% Function Definition  
  87.       
  88.     function [pMiu pPi pSigma] = init_params()  
  89.         pMiu = centroids; %k*D, 即k类的中心点  
  90.         pPi = zeros(1, K); %k类GMM所占权重（influence factor）  
  91.         pSigma = zeros(D, D, K); %k类GMM的协方差矩阵，每个是D*D的  
  92.    
  93.         % 距离矩阵，计算N*K的矩阵（x-pMiu）^2 = x^2+pMiu^2-2*x*Miu  
  94.         distmat = repmat(sum(X.*X, 2), 1, K) + ... %x^2, N*1的矩阵replicateK列  
  95.             repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...%pMiu^2，1*K的矩阵replicateN行  
  96.             2*X*pMiu';  
  97.         [~, labels] = min(distmat, [], 2);%Return the minimum from each row  
  98.    
  99.         for k=1:K  
  100.             Xk = X(labels == k, :);  
  101.             pPi(k) = size(Xk, 1)/N;  
  102.             pSigma(:, :, k) = cov(Xk);  
  103.         end  
  104.     end  
  105.    
  106.     function Px = calc_prob()   
  107.         %Gaussian posterior probability   
  108.         %N(x|pMiu,pSigma) = 1/((2pi)^(D/2))*(1/(abs(sigma))^0.5)*exp(-1/2*(x-pMiu)'pSigma^(-1)*(x-pMiu))  
  109.         Px = zeros(N, K);  
  110.         for k = 1:K  
  111.             Xshift = X-repmat(pMiu(k, :), N, 1); %X-pMiu  
  112.             inv_pSigma = inv(pSigma(:, :, k));  
  113.             tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);  
  114.             coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));  
  115.             Px(:, k) = coef * exp(-0.5*tmp);  
  116.         end  
  117.     end  
  118. end  

2. gmm_accuracy.m调用gmm.m，计算准确率：

[cpp] view plain copy
 
  1. function [ Accuracy ] = gmm_accuracy( Data_fea, gnd_label, K )  
  2. %Calculate the accuracy Clustered by GMM model  
  3.   
  4. px = gmm(Data_fea,K);  
  5. [~, cls_ind] = max(px,[],1); %cls_ind = cluster label   
  6. Accuracy = cal_accuracy(cls_ind, gnd_label);  
  7.   
  8.     function [acc] = cal_accuracy(gnd,estimate_label)  
  9.         res = bestMap(gnd,estimate_label);  
  10.         acc = length(find(gnd == res))/length(gnd);  
  11.     end  
  12.   
  13. end  

3. 主函数调用

gmm_acc = gmm_accuracy(fea,gnd,N_classes);

写了本文进行总结后自己很受益，也希望大家可以好好YM下上面pluskid的gmm.m，不光是算法，其中的矩阵处理代码也写的很简洁，很值得学习。

另外看了两份东西非常受益，一个是pluskid大牛的《漫谈 Clustering (3): Gaussian Mixture Model》，一个是JerryLead的EM算法详解，大家有兴趣也可以看一下，写的很好。

关于Machine Learning更多的学习资料与相关讨论将继续更新，敬请关注本博客和新浪微博Sophia_qing。

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/xfzhang/p/3789559.html