高斯混合模型

http://blog.csdn.net/crzy_sparrow/article/details/7413019

好吧，我承认这个题目有点噱头，其实本文要讲的一般的高斯混合模型，基于matlab实现，没有涉及到opencv。之所以作为opencv的学习笔记之一是因为之后打算讲下基于高斯混合模型的背景建模（实现目标跟踪），所以把这个也放上来了。

    混合高斯模型的原理说白了就一句话：任意形状的概率分布都可以用多个高斯分布函数去近似。换句话说，高斯混合模型（GMM）可以平滑任意形状的概率分布。其参数求解方法一般使用极大似然估计求解，但使用极大似然估计法往往不能获得完整数据（比如样本已知，但样本类别（属于哪个高斯分布）未知），于是出现了EM（最大期望值）求解方法。

    虽然上面说的简单，但是混合高斯模型和EM求解的理论还是比较复杂的，我把我所找到的我认为能够快速掌握高斯混合模型的资料打包到了附件中，大家可以去下载，了解混合高斯模型以及EM的完整推导过程。

附件下载地址：

http://download.csdn.net/detail/crzy_sparrow/4187994

大致抽取下高斯混合模型的重要概念：

1）任意数据分布可用高斯混合模型（M个单高斯）表示（（1）式）

（1）

其中：

    （2）

   （3）   表示第j个高斯混合模型

2）N个样本集X的log似然函数如下：

     （4）

其中：

参数：，   （5）

下面具体讲讲基于EM迭代的混合高斯模型参数求解算法流程：

1）初始参数由k-means（其实也是一种特殊的高斯混合模型）决定：

[plain]  view plain copy

1. function [Priors, Mu, Sigma] = EM_init_kmeans(Data, nbStates)  
2. % Inputs -----------------------------------------------------------------  
3. %   o Data:     D x N array representing N datapoints of D dimensions.  
4. %   o nbStates: Number K of GMM components.  
5. % Outputs ----------------------------------------------------------------  
6. %   o Priors:   1 x K array representing the prior probabilities of the  
7. %               K GMM components.  
8. %   o Mu:       D x K array representing the centers of the K GMM components.  
9. %   o Sigma:    D x D x K array representing the covariance matrices of the   
10. %               K GMM components.  
11. % Comments ---------------------------------------------------------------  
12. %   o This function uses the 'kmeans' function from the MATLAB Statistics   
13. %     toolbox. If you are using a version of the 'netlab' toolbox that also  
14. %     uses a function named 'kmeans', please rename the netlab function to  
15. %     'kmeans_netlab.m' to avoid conflicts.   
16.   
17. [nbVar, nbData] = size(Data);  
18.   
19. %Use of the 'kmeans' function from the MATLAB Statistics toolbox  
20. [Data_id, Centers] = kmeans(Data', nbStates);   
21. Mu = Centers';  
22. for i=1:nbStates  
23.   idtmp = find(Data_id==i);  
24.   Priors(i) = length(idtmp);  
25.   Sigma(:,:,i) = cov([Data(:,idtmp) Data(:,idtmp)]');  
26.   %Add a tiny variance to avoid numerical instability  
27.   Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(nbVar,1));  
28. end  
29. Priors = Priors ./ sum(Priors);  

2）E步（求期望）

求取：

（7）

其实，这里最贴切的式子应该是log似然函数的期望表达式

事实上，3）中参数的求取也是用log似然函数的期望表达式求偏导等于0解得的简化后的式子，而这些式子至于（7）式有关，于是E步可以只求（7）式，以此简化计算，不需要每次都求偏导了。

[plain]  view plain copy

1. %% E-step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
2. for i=1:nbStates  
3.   %Compute probability p(x|i)  
4.   Pxi(:,i) = gaussPDF(Data, Mu(:,i), Sigma(:,:,i));  
5. end  
6. %Compute posterior probability p(i|x)  
7. Pix_tmp = repmat(Priors,[nbData 1]).*Pxi;  
8. Pix = Pix_tmp ./ repmat(sum(Pix_tmp,2),[1 nbStates]);  
9. %Compute cumulated posterior probability  
10. E = sum(Pix);  

[plain]  view plain copy

1. function prob = gaussPDF(Data, Mu, Sigma)  
2. % Inputs -----------------------------------------------------------------  
3. %   o Data:  D x N array representing N datapoints of D dimensions.  
4. %   o Mu:    D x K array representing the centers of the K GMM components.  
5. %   o Sigma: D x D x K array representing the covariance matrices of the   
6. %            K GMM components.  
7. % Outputs ----------------------------------------------------------------  
8. %   o prob:  1 x N array representing the probabilities for the   
9. %            N datapoints.       
10.   
11. [nbVar,nbData] = size(Data);  
12.   
13. Data = Data' - repmat(Mu',nbData,1);  
14. prob = sum((Data*inv(Sigma)).*Data, 2);  
15. prob = exp(-0.5*prob) / sqrt((2*pi)^nbVar * (abs(det(Sigma))+realmin));  

3）M步（最大化步骤）

更新权值：

更新均值：

更新方差矩阵：

[cpp]  view plain copy

1. %% M-step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
2.  for i=1:nbStates  
3.    %Update the priors  
4.    Priors(i) = E(i) / nbData;  
5.    %Update the centers  
6.    Mu(:,i) = Data*Pix(:,i) / E(i);  
7.    %Update the covariance matrices  
8.    Data_tmp1 = Data - repmat(Mu(:,i),1,nbData);  
9.    Sigma(:,:,i) = (repmat(Pix(:,i)',nbVar, 1) .* Data_tmp1*Data_tmp1') / E(i);  
10.    %% Add a tiny variance to avoid numerical instability  
11.    Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(nbVar,1));  
12.  end  

4）截止条件

不断迭代EM步骤，更新参数，直到似然函数前后差值小于一个阈值，或者参数前后之间的差（一般选择欧式距离）小于某个阈值，终止迭代，这里选择前者。附上结合EM步骤的代码：

[cpp]  view plain copy

1. while 1  
2.   %% E-step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
3.   for i=1:nbStates  
4.     %Compute probability p(x|i)  
5.     Pxi(:,i) = gaussPDF(Data, Mu(:,i), Sigma(:,:,i));  
6.   end  
7.   %Compute posterior probability p(i|x)  
8.   Pix_tmp = repmat(Priors,[nbData 1]).*Pxi;  
9.   Pix = Pix_tmp ./ repmat(sum(Pix_tmp,2),[1 nbStates]);  
10.   %Compute cumulated posterior probability  
11.   E = sum(Pix);  
12.   %% M-step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
13.   for i=1:nbStates  
14.     %Update the priors  
15.     Priors(i) = E(i) / nbData;  
16.     %Update the centers  
17.     Mu(:,i) = Data*Pix(:,i) / E(i);  
18.     %Update the covariance matrices  
19.     Data_tmp1 = Data - repmat(Mu(:,i),1,nbData);  
20.     Sigma(:,:,i) = (repmat(Pix(:,i)',nbVar, 1) .* Data_tmp1*Data_tmp1') / E(i);  
21.     %% Add a tiny variance to avoid numerical instability  
22.     Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(nbVar,1));  
23.   end  
24.   %% Stopping criterion %%%%%%%%%%%%%%%%%%%%  
25.   for i=1:nbStates  
26.     %Compute the new probability p(x|i)  
27.     Pxi(:,i) = gaussPDF(Data, Mu(:,i), Sigma(:,:,i));  
28.   end  
29.   %Compute the log likelihood  
30.   F = Pxi*Priors';  
31.   F(find(F<realmin)) = realmin;  
32.   loglik = mean(log(F));  
33.   %Stop the process depending on the increase of the log likelihood   
34.   if abs((loglik/loglik_old)-1) < loglik_threshold  
35.     break;  
36.   end  
37.   loglik_old = loglik;  
38.   nbStep = nbStep+1;  
39. end  

结果测试（第一幅为样本点集，第二幅表示求取的高斯混合模型，第三幅为回归模型）：

忘了说了，传统的GMM和k-means一样，需要给定K值（即:要有几个高斯函数）。