呵呵，第一次写技术类博文，请多指教。因为做GMM时遇到些小麻烦，所以干脆把问题解决的过程都写下来，帮助学习GMM的童鞋省些工夫。闲话少说，开始。

   使用的数据库是http://www.zjucadcg.cn/dengcai/Data/FaceData.html所提供的。这个数据是YaleB中的10个subject加上Extended YaleB的28个subject，从原来的9个pose里面取出frontal pose，64种光照变化保留，所以一共有38*64 = 2432个。有些图象在获取时损坏被去掉，所以还剩2414个。YaleB_32x32.mat里面包含了所有的脸部图象（即fea，2414x1024的那个，因为每个图象被chop为32x32，所以reshape为1x1024）。YaleB_32x32_corrected.mat里面是经过光照修正的脸部图象。此网站还提供了一些split，就是随机序号，用来生成训练集和测试集，20Train就是每个subject的64个图象中的20个用于训练，依此类推。

   read_data是读入数据并归一化（除以255）。注意归一化很重要，如果不做结果会很差。

   做GMM之前先使用PCA把1024维的特征降到150维，否则难以收敛。coefs是主成分，scores每个sample是投射到低维空间后的坐标。

   显示GMM center时要把mu乘以coefs还原为1024维。

   我的GMM是采用了general的假设，即协方差未知、均值未知，但是假设所有的协方差矩阵均为对角阵（即除了对角线元素都为0）。也可以弱化为假设所有类协方差均相等，但是我没有试。关于GMM的原理和计算，请参考CMU教授Andrew Moore的讲义：http://www.autonlab.org/tutorials/gmm.html。我觉得既有趣又好懂。GMM基本上是和EM（Expectation Maximazation）绑在一起的，所以如果不知道什么是EM也请看此文。

   几个关键之处：

   1. 使用kmeans初始化，即GMM的各类中心设为kmeans给出的c个中心；

   2. 各协方差矩阵初始化为单位阵，不能设为kmeans给出的协方差，否则会不收敛（在这个上面忙了很久）；

   3. sigma = sigma + 1e-4是防止sigma中出现0，只要加上一个小正数即可，不一定为1e-4，这个数是经验；

   4. thresh，即判定收敛与否的阈值，也是试出来的……但是MATLAB的gmdistribution.fit没有让用户设定阈值，它会自己判断，我想可能是判断error占mu的norm加sigma的norm的百分比吧，不知道它具体怎么实现的。

  总而言之，GMM的收敛问题很难办，一个细节没搞好就不收敛了……最后效果还不错，有几张脸和subject很像，但是有些脸基本上是全黑的，无法辨认；这个问题大概是因为GMM把所有光照不足黑乎乎的脸都认为是一类了。无监督学习容易出现这种问题。

  下面附GMM代码：

  主程序：

  main_gmm.m

===============================================

[fea_Train fea_Test gnd_Train gnd_Test] = read_data(50, 1);
fea_All = cat(1, fea_Train, fea_Test);

[coefs scores] = princomp(fea_All);　　% princomp是MATLAB的PCA函数
coefs_gmm = coefs(:,1:150);　　　　　% 降到150维
scores_gmm = scores(:,1:150);
[mu sigma] = my_gmm(scores_gmm, 38);

gmm_face = mu*coefs_gmm';

for i = 1:38
    subplot(4,10,i);
    Y = reshape(gmm_face(i,:),[32 32]);
    imagesc(Y);axis off;colormap(gray)　　% colormap会自动做直方图均衡化
    title(['Subject ', num2str(i)]);
    hold on
end

=============================================

  函数：

  my_gmm.m

function [mu sigma] = my_gmm(X,c)
% Fit Gaussian Mixture Model for input data X, c classes. For simplicity
% and stability covariance are set to be diagonal.
% Input:
% X:[nxd]. n observations by d variables.
% c: no. of classes.
% Output:
% mu: [cxd]. Each row is the mean of a class.
% sigma:[dxc].

% Initialization
% vl_kmeans can be replaced by MATLAB's kmeans, but this one is
% significantly faster. The toolbox can be downloaded at vlfeat.org.
[mu idx] = vl_kmeans(X',c);
clear idx
mu = mu';
mu_prev = mu;

% A priori probability of each class.
[n d] = size(X);
prior = ones(c,1)/c;
prior_prev = prior;

% A posterior probability of each sample belonging to each class.
post = ones(n,c)/c;
% set sigma as an array representing diagonal matrix
sigma = ones(1,d,c);sigma_prev = sigma;
likeli = zeros(n,c);

maxIter = 100;
thresh = 1;
error = 0;

for t=1:maxIter
    % E-step
    for i = 1:c
        post(:,i)= mvnpdf(X,mu_prev(i,:),sigma_prev(:,:,i))*prior_prev(i);
        likeli(:,i) = mvnpdf(X,mu_prev(i,:),sigma_prev(:,:,i));
    end
    for j = 1:n
        post(j,:) = post(j,:)/sum(post(j,:));
    end   
    % M-step   
    mu = post'*X./repmat(sum(post)',1,d);
    for i = 1:c
       sigma(:,:,i) = diag((X'-repmat(mu(i,:)',1,n))*diag(post(:,i))*(X'-repmat(mu(i,:)',1,n))'/sum(post(:,i)));  % 这个地方当时卡了很久T_T乘后验概率矩阵的时候一定要将其对角化！否则sigma不正定！diag(post(:,i))就是取每个数据属于第i类的后验概率向量（2414x1）将其变为2414x2414维的对角阵，这样X‘PX才与公式中的相等
    end
   
    % Regularize
    sigma = sigma + 1e-4;
    for i = 1:c
        prior(i) = sum(post(:,i))/n;
    end
   
    for i = 1:c
        error = error + norm(mu(i,:)-mu_prev(i,:)) + norm(sigma(:,:,i)-sigma_prev(:,:,i));
    end
   
    if error < thresh
        break
    end
    mu_prev = mu;
    sigma_prev = sigma;
    prior_prev = prior;
    error = 0;
end

fprintf('Iteration times: %d\n',t);

============================================================

  利用gmdistribution.fit函数进行参数估计的结果：

[fea_Train fea_Test gnd_Train gnd_Test] = read_data(50, 1);
fea_All = cat(1, fea_Train, fea_Test);
options = statset('Display', 'final');
obj = gmdistribution.fit(fea_All, 38, 'Options', options, 'CovType', 'diagonal', 'Regularize', 1e-4);
for i = 1:38
    subplot(4,10,i);
     Y = reshape(obj.mu(i,:),[32 32]);
     imagesc(Y);axis off;colormap(gray);
     title(['Subject ', num2str(i)]);
     hold on

end

=============================================================

  如果此文对您的学习有任何帮助，欢迎在下面回复告诉我！谢谢！