最近在看高斯混合模型（Gaussian Mixture Model, GMM），涉及到高斯分布的参数。为此特意回顾了概率论的二维高斯分布的相关概念，并分析了参数对二维高斯分布曲面的影响。
 
1、多维高斯分布的概率密度函数
 
　　　 
 　　多维变量$X = ({x_1},{x_2},...{x_n})$的联合概率密度函数为： 
 　　　　　　　$f(X) = \frac{1}{{{{(2\pi )}^{d/2}}{{\left| \Sigma  \right|}^{1/2}}}}\exp [ - \frac{1}{2}{(X - u)^T}{\Sigma ^{ - 1}}(X - u)],X = ({x_1},{x_2}...{x_n})$ 
 　　其中： 
 　　d：变量维度。对于二维高斯分布，有d=2; 
 　　$u = \left( \begin{array}{l}  
{u_1}\  
{u_2}\  
…\  
{u_n}  
\end{array} \right)$：各位变量的均值； 
 　　$\Sigma$：协方差矩阵，描述各维变量之间的相关度。对于二维高斯分布，有：

 
   
    Σ=(δ11δ21δ12δ22)
   
    
     
      Σ
     
     
      =
     
     
      
       (
      
      
       
        
         
          
           
            
             
              δ
             
             
              
               11
              
             
            
           
          
         
         
          
           
            
             
              δ
             
             
              
               12
              
             
            
           
          
         
        
        
         
          
           
            
             
              δ
             
             
              
               21
              
             
            
           
          
         
         
          
           
            
             
              δ
             
             
              
               22
              
             
            
           
          
         
        
       
      
      
       )
      
     
    

 $$\Sigma  = \left( {\begin{array}{*{20}{c}}
{{\delta _{11}}}&{{\delta _{12}}}\\
{{\delta _{21}}}&{{\delta _{22}}}
\end{array}} \right)$$  
 　　后文主要分析均值和协方差矩阵对二维高斯分布的影响。 
 
 
2、均值和协方差矩阵对二维高斯分布的影响
 
2.1 $u = \left( \begin{array}{l}  
0\  
0  
\end{array} \right),\Sigma  = \left( {\begin{array}{*{20}{c}}
3&0\\
0&3
\end{array}} \right)$ 
  
 2.2 $u = \left( \begin{array}{l}  
4\  
4  
\end{array} \right),\Sigma  = \left( {\begin{array}{*{20}{c}}
3&0\\
0&3
\end{array}} \right)$ 
  
 2.3 $u = \left( \begin{array}{l}  
0\  
0  
\end{array} \right),\Sigma  = \left( {\begin{array}{*{20}{c}}
3&0\\
0&{10}
\end{array}} \right)$ 
  
 2.4 $u = \left( \begin{array}{l}  
0\  
0  
\end{array} \right),\Sigma  = \left( {\begin{array}{*{20}{c}}
1&{0.8}\\
{0.8}&1
\end{array}} \right)$ 
  
 2.5 $u = \left( \begin{array}{l}  
0\  
0  
\end{array} \right),\Sigma  = \left( {\begin{array}{*{20}{c}}
1&{ - 0.8}\\
{ - 0.8}&1
\end{array}} \right)$ 
 
 
3、总结
 
①均值表征的是各维变量的中心，其对二维高斯曲面的影响较好理解，它使得整个二维高斯曲面在xoy平面上移动； 
 ②对于协方差矩阵，对角线上的两个元素，即${\delta _{11}}$和${\delta _{22}}$表征的是x维和y维变量的方差，决定了整个高斯曲面在某一维度上的“跨度”，方差越大，“跨度”越大； 
 ③协方差矩阵的斜对角线上面的两个元素，即${\delta _{12}}$和${\delta _{21}}$（${\delta _{12}}$=${\delta _{21}}$）表征的是各维变量之间的相关性：${\delta _{12}}$>0说明x与y呈正相关（x越大，y越大），其值越大，正相关程度越大；${\delta _{12}}$<0呈负相关；否则不相关。

   参考文章
   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
   《统计学习方法》 李航
   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。
   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义
高斯混合函数
        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）
        二维高斯混合函数，如(2)式所示。
   （2）

       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？
协方差矩阵
       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。
        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。
 
         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。


   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。
    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。
EM算法
   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。
   我主要写一下自己的理解，也为自己以后复习用。
   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。
   基本思想：
   E步，假设参数已知，去估计隐藏变量的期望；
   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。
  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。
学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：
1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，

%产生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,'*');






  2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑
   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。
 程序：GMM.m

%程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.
%Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。
%Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据
%产生2个二维正态数据
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
%生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
%显示
scatter(X(:,1),X(:,2),10,'.');
%====================
K=2;
[N,D]=size(X);
Gamma=zeros(N,K);
Psi=zeros(N,K);
Mu=zeros(K,D);
LM=zeros(K,D);
Sigma =zeros(D, D, K); 
Pi=zeros(1,K);
%选择随机的两个样本点作为期望迭代初值
Mu(1,:)=X(randint(1,1,[1 N/2]),:);
Mu(2,:)=X(randint(1,1,[N/2 N]),:);
%所有数据的协方差作为协方差初值
for k=1:K
  Pi(k)=1/K;
  Sigma(:, :, k)=cov(X);
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
while true
%Estimation Step  
for k = 1:K
  Y = X - repmat(Mu(k,:),N,1);
  Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值
end
Gamma_SUM=zeros(D,D);
for j = 1:N
  for k=1:K
  Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值
  end
end
%Maximization Step
for k = 1:K
%update Mu
  Mu_SUM= zeros(1,D);
  for j=1:N
     Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   
  end
  Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
%update Sigma
 Sigma_SUM= zeros(D,D);
 for j = 1:N
   Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
 end
 Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
 %update Pi
 Pi_SUM=0;
 for j=1:N
   Pi_SUM=Pi_SUM+Gamma(j,k);
 end
 Pi(1,k)=Pi_SUM/N;
end

R_Mu=sum(sum(abs(LMu- Mu)));
R_Sigma=sum(sum(abs(LSigma- Sigma)));
R_Pi=sum(sum(abs(LPi- Pi)));
R=R_Mu+R_Sigma+R_Pi;
if ( R<1e-10)
    disp('期望');
    disp(Mu);
    disp('协方差矩阵');
    disp(Sigma);
    disp('权值系数');
    disp(Pi);
    obj=gmdistribution(Mu,Sigma);
    figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
    break;
end
 LMu=Mu;
 LSigma=Sigma;
 LPi=Pi;
end
%=====================


%%%%%%%%运行结果%%%%%%%
%
% 期望
%    1.0334    2.0479
%    -0.9703   -0.9879
% 
% 协方差矩阵
% 
% (:,:,1) =
% 
%     0.9604   -0.0177
%    -0.0177    0.4463
% 
% 
% (:,:,2) =
% 
%     0.9976    0.0667
%     0.0667    1.0249
% 
% 权值系数
%     0.4935    0.5065
%%%%%%%%%%%%%%%%%%%%%%%%%
 
     

参考文章
 
   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
 
   《统计学习方法》 李航
 
   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。
 
   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义
 
高斯混合函数
 
        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）
 
 
二维高斯混合函数，如(2)式所示。
 
 
   （2）
 
       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？
 
协方差矩阵
 
       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。
 
        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。
 
 
 
         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。
 
   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。
 
    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。
 
EM算法
 
   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。
 
   我主要写一下自己的理解，也为自己以后复习用。
 
   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。
 
   基本思想：
 
   E步，假设参数已知，去估计隐藏变量的期望；
 
   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。
 
  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。
 
学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：
 
1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，
 

 %产生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,'*');
 

   2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑
 
   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。
 
 程序：GMM.m
 

 %程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.
 %Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。
 %Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据
 %产生2个二维正态数据
 MU1 = [1 2];
 SIGMA1 = [1 0; 0 0.5];
 MU2 = [-1 -1];
 SIGMA2 = [1 0; 0 1];
 %生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数
 X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
 %显示
 scatter(X(:,1),X(:,2),10,'.');
 %====================
 K=2;
 [N,D]=size(X);
 Gamma=zeros(N,K);
 Psi=zeros(N,K);
 Mu=zeros(K,D);
 LM=zeros(K,D);
 Sigma =zeros(D, D, K); 
 Pi=zeros(1,K);
 %选择随机的两个样本点作为期望迭代初值
 Mu(1,:)=X(randint(1,1,[1 N/2]),:);
 Mu(2,:)=X(randint(1,1,[N/2 N]),:);
 %所有数据的协方差作为协方差初值
 for k=1:K
   Pi(k)=1/K;
   Sigma(:, :, k)=cov(X);
 end
 LMu=Mu;
 LSigma=Sigma;
 LPi=Pi;
 while true
 %Estimation Step  
 for k = 1:K
   Y = X - repmat(Mu(k,:),N,1);
   Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值
 end
 Gamma_SUM=zeros(D,D);
 for j = 1:N
   for k=1:K
   Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值
   end
 end
 %Maximization Step
 for k = 1:K
 %update Mu
   Mu_SUM= zeros(1,D);
   for j=1:N
      Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   
   end
   Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
 %update Sigma
  Sigma_SUM= zeros(D,D);
  for j = 1:N
    Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
  end
  Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
  %update Pi
  Pi_SUM=0;
  for j=1:N
    Pi_SUM=Pi_SUM+Gamma(j,k);
  end
  Pi(1,k)=Pi_SUM/N;
 end
 
R_Mu=sum(sum(abs(LMu- Mu)));
 R_Sigma=sum(sum(abs(LSigma- Sigma)));
 R_Pi=sum(sum(abs(LPi- Pi)));
 R=R_Mu+R_Sigma+R_Pi;
 if ( R<1e-10)
     disp('期望');
     disp(Mu);
     disp('协方差矩阵');
     disp(Sigma);
     disp('权值系数');
     disp(Pi);
     obj=gmdistribution(Mu,Sigma);
     figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
     break;
 end
  LMu=Mu;
  LSigma=Sigma;
  LPi=Pi;
 end
 %=====================
 

 %%%%%%%%运行结果%%%%%%%
 %
 % 期望
 %    1.0334    2.0479
 %    -0.9703   -0.9879
 % 
 % 协方差矩阵
 % 
 % (:,:,1) =
 % 
 %     0.9604   -0.0177
 %    -0.0177    0.4463
 % 
 % 
 % (:,:,2) =
 % 
 %     0.9976    0.0667
 %     0.0667    1.0249
 % 
 % 权值系数
 %     0.4935    0.5065
 %%%%%%%%%%%%%%%%%%%%%%%%%
  
 
参考文章
 
   http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
 
   《统计学习方法》 李航
 
   EM算法是一种迭代算法，1977年由Dempster等人总结出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望(expectation);M步，求极大(maximization).所以这一算法称为期望极大算法(expectation maximization algorithm),简称EM算法。
 
   首先准备一些预备知识，如：二维高斯混合函数的形式，协方差的定义
 
高斯混合函数
 
        对于一维混合高斯函数，如(1)式所示。(所谓的混合高斯，就是将一些高斯函数利用一些权值系数组合起来)  （1）
 
        二维高斯混合函数，如(2)式所示。
 
   （2）
 
       可以看出二维高斯混合模型，其中的d为数据的维度，在这里等于2；∑  代表协方差矩阵，两行两列；这里的x以及μ都是二维的数据，用矩阵表示就是一行两列。这里提到了协方差矩阵，那么什么是协方差矩阵，与方差有什么区别？
 
协方差矩阵
 
       理解协方差矩阵的关键在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间。拿到一个样本矩阵，我们最先明确的是一行是一个样本还是一个维度（重要）。
 
        方差表示的是每个样本值与全体样本值的平均数之差的平方值的平均数，反映了一组数据的离散程度。
 
 
 
         既然有了方差还要协方差干啥？显然，方差只能衡量一维数据的离散程度，对于二维数据，就用到了协方差，如下式所示。其中X，Y代表一组数据的两个维度。COV(X,Y)的值大于零，表示两个维度的数据成正相关，小于零则表示成负相关，等于零表示两个维度的数据相互独立。显然当X=Y时，COV(X,X)计算的就是方差。由协方差的定义可以得到COV(X,Y)=COV(Y,X)。
 
   二维数据的协方差矩阵表示为 ，显然协方差矩阵为对角矩阵。
 
    三维数据的协方差矩阵表示为，同理知道N维数据的协方差矩阵。
 
EM算法
 
   假设一组数据x1,x2,...,xn是来自上述二维高斯混合模型，这里我们利用EM算法确定各个高斯部件的参数，充分拟合给定的数据。算法目的具体的理论知识与算法步骤参考《统计学习方法》李航一书中对于EM算法的讲解(李航的这一本书，绝对值得学习)。
 
   我主要写一下自己的理解，也为自己以后复习用。
 
   EM算法其实就是概率模型参数的极大似然估计法，只不过存在一些隐变量，从而不能像一般的概率模型参数估计一样直接求似然函数的最大值，因为这里的似然函数存在隐变量，因而不存在解析解，就不能直接对 似然函数求偏导来求极大值，只能迭代求解。
 
   基本思想：
 
   E步，假设参数已知，去估计隐藏变量的期望；
 
   M步，利用E步求得的隐藏变量的期望，根据最大似然估计求得参数，然后新得到的参数值重新被用于E步。。。直至收敛到局部最优解。
 
  需要注意的一点是：理论中，E步最重要的是求Q函数（李航 9.29式），即完全数据的对数似然函数的期望，M步要做的是用Q函数对各参数求偏导，从而确定满足Q函数最大的各参数的值。而在程序中，因为M步Q函数对参数求偏导后，各参数的值由Q函数中的Gamma(j,k)决定，所以E步只需计算Gamma(j,k)即可，Gamma(j,k)代表第j个样本属于第k个高斯模型的概率，这样M步各参数的值也就确定了，所以切记实例理论与程序实现的区别，好像是程序只是利用了理论的精简结果。
 
学习完理论之后，最好写个程序练习一下，加深自己的理解。下面的程序(matlab编写)，是这样一个过程：
 
1.首先产生两组数据点。每一组数据点都是符合给定均值Mu和协方差矩阵∑的。程序和图如下所示，
 

 %产生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,'*');
 

   2.对于上述产生的数据，利用EM算法预测这组数据对应的均值μ和协方差矩阵∑
 
   对于各参数初值的设置：两个均值初始值各设置为两组数据点中的随机一个样本点，权值系数初始值设置为1/k,k是高斯混合模型中模型的个数，协方差矩阵初值设置为所有数据的协方差。
 
 程序：GMM.m
 

 %程序中Psi对应二维高斯函数，其中的(x-μ)与其转置的顺序与上述提到的高斯函数不同，这样是为了保证矩阵可乘性，不影响结果.
 %Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。
 %Mu为期望，Sigma为协方差矩阵%Pi为各模型的权值系数%产生2个二维正态数据
 %产生2个二维正态数据
 MU1 = [1 2];
 SIGMA1 = [1 0; 0 0.5];
 MU2 = [-1 -1];
 SIGMA2 = [1 0; 0 1];
 %生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数
 X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
 %显示
 scatter(X(:,1),X(:,2),10,'.');
 %====================
 K=2;
 [N,D]=size(X);
 Gamma=zeros(N,K);
 Psi=zeros(N,K);
 Mu=zeros(K,D);
 LM=zeros(K,D);
 Sigma =zeros(D, D, K); 
 Pi=zeros(1,K);
 %选择随机的两个样本点作为期望迭代初值
 Mu(1,:)=X(randint(1,1,[1 N/2]),:);
 Mu(2,:)=X(randint(1,1,[N/2 N]),:);
 %所有数据的协方差作为协方差初值
 for k=1:K
   Pi(k)=1/K;
   Sigma(:, :, k)=cov(X);
 end
 LMu=Mu;
 LSigma=Sigma;
 LPi=Pi;
 while true
 %Estimation Step  
 for k = 1:K
   Y = X - repmat(Mu(k,:),N,1);
   Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y')));      %Psi第一列代表第一个高斯分布对所有数据的取值
 end
 Gamma_SUM=zeros(D,D);
 for j = 1:N
   for k=1:K
   Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');                                               %Psi的第一行分别代表两个高斯分布对第一个数据的取值
   end
 end
 %Maximization Step
 for k = 1:K
 %update Mu
   Mu_SUM= zeros(1,D);
   for j=1:N
      Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);   
   end
   Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
 %update Sigma
  Sigma_SUM= zeros(D,D);
  for j = 1:N
    Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
  end
  Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
  %update Pi
  Pi_SUM=0;
  for j=1:N
    Pi_SUM=Pi_SUM+Gamma(j,k);
  end
  Pi(1,k)=Pi_SUM/N;
 end
 
R_Mu=sum(sum(abs(LMu- Mu)));
 R_Sigma=sum(sum(abs(LSigma- Sigma)));
 R_Pi=sum(sum(abs(LPi- Pi)));
 R=R_Mu+R_Sigma+R_Pi;
 if ( R<1e-10)
     disp('期望');
     disp(Mu);
     disp('协方差矩阵');
     disp(Sigma);
     disp('权值系数');
     disp(Pi);
     obj=gmdistribution(Mu,Sigma);
     figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
     break;
 end
  LMu=Mu;
  LSigma=Sigma;
  LPi=Pi;
 end
 %=====================
 

 %%%%%%%%运行结果%%%%%%%
 %
 % 期望
 %    1.0334    2.0479
 %    -0.9703   -0.9879
 % 
 % 协方差矩阵
 % 
 % (:,:,1) =
 % 
 %     0.9604   -0.0177
 %    -0.0177    0.4463
 % 
 % 
 % (:,:,2) =
 % 
 %     0.9976    0.0667
 %     0.0667    1.0249
 % 
 % 权值系数
 %     0.4935    0.5065
 %%%%%%%%%%%%%%%%%%%%%%%%%
  
      
 
 
     
 
--------------------- 
 作者：lpkinging 
 原文：https://blog.csdn.net/u011582757/article/details/70003351 

https://blog.csdn.net/lin_limin/article/details/81048411
 
https://blog.csdn.net/farmwang/article/details/78699926
 
http://blog.sina.com.cn/s/blog_a36a563e0102y2ec.html
 
http://blog.sina.com.cn/s/blog_a36a563e0102y2ec.html
 
 
 
 
为X,Y的相关系数！为0；σ1=σ2=σ
 
二维高斯曲面的公式（x,y代表像素的模板坐标，模板中心位置为原点）
 
 
 
 
根据这个公式，我们可以计算得到不同σ的高斯模板。下面是C语言程序实现：
 
当σ即半径为0.7时：
 
 
 
 
 
 
#include "stdafx.h"
 
 
#include <fstream>
 
 
#include <math.h>
 
 
#include <iomanip>
 
 
using namespace std;
 
 
#define N 4
 
 
#define PI 3.141592653
 
 
int _tmain(int argc, _TCHAR* argv[])
 
 
{
 
 
char* path = "C:\\path1.txt";
 
 
ofstream fout(path);
 
 
double a[2 * N + 1][2 * N + 1]; // 高斯模板;
 
 
double r = 1; // 高斯半径;
 
 
double A = 1 / (2 * PI * r * r);
 
 
if (fout)
 
 
{
 
 
for (int i = -1 * N; i <= N; i++)
 
 
{
 
 
for (int j = -1 * N; j <= N; j++)
 
 
{
 
 
a[i + N][j + N] = A*exp((-1)*(i*i + j*j) / (2 * r*r));
 
 
fout << setiosflags(ios::fixed) << setprecision(6) << a[i + N][j + N] << " ";
 
 
}
 
 
fout << endl;
 
 
}
 
 
}
 
 
fout.close();
 
 
return 0;
 
 
}
 
最后得到的结果如下图所示
 
 
那么，我们还有个问题，如何确定模板的大小与标准差之间的关系。经过我们的不断验证，即改变上述的σ和模板大小的值，可以得知。当σ越大时，要求的模板也就是越大。即当σ为1时，我们可以得到如下的高斯模板：
 
 
 
 
上图的意思是，与水平面平行的x，y平面，也就是上面高斯模板中的坐标(x,y)，z轴表示的灰度值，也就是上面高斯模板中的灰度值。
 
通过上图我们可以得知，7*7的模板已经不能满足我们，我们需要9*9的模板，这时如果用5*5的模板，则会取中间部分的5*5模板。模板和σ的选择取决于我们图像的大小。如果图像是几万的分辨率，则需要高斯模板和σ值都要很大，如果几十的分辨率，用很小的如3*3的模板或者5*5的模板就可以了，这时σ一般可以取0.5左右（也就是matlab里fspecial函数里默认的σ值）。
 
得到的绘出的高斯曲面：
 
 
 
 
现在众所周知用的3*3或者5*5的模板，都是对高斯曲面的一个整数除法形式的近似，也就是对高斯半径约为0.849时的一个近似。
 
1 2 1 
 2 4 2 /16 
 1 2 1
 
16代表的是权重，即模板内所有的数之和。5*5的模板如下，也是对高斯曲面的一个整数除法形式的近似。
 
 
 
 
最近在看晓川老（shi）师(shu)的博士论文，接触了混合高斯模型（Gaussian mixture model, GMM）和EM(Expectation Maximization)算法，不禁被论文中庞大的数学公式所吓退。本文通过查阅相关资料，在复杂巧妙的推理公式中融入了自己的理解，详细梳理了混合高斯模型和EM算法。
 1 单高斯模型(Gaussian single model, GSM）
 
  简单回顾一下概率论讲过的高斯模型。
 高斯模型是一种常用的变量分布模型，在数理统计领域有着广泛的应用（……好吧读了这么多年书没白费，教科书般的话语已植入骨髓）。一维高斯分布的概率密度函数如下：
 f(x)=12π√σexp(−(x−μ)22σ2)(1)
 f(x)=2π
 ​σ1​exp(−2σ2(x−μ)2​)(1)
 μμ和σ2σ2分别是高斯分布的均值和方差。
 譬如将男生身高视为变量X, 假设男生的身高服从高斯分布，则X∼N(μ,σ2)X∼N(μ,σ2)，女生亦如此。只是男女生身高分布可能具有不同的均值和方差。图1是从谷歌图库中搜索到的男女生身高分布图，来源不清，个人觉得男生的均值身高虚高……四个记号分别表示0~3σ
 
σ准则。
 图1 男女生身高分布差异
 图1 男女生身高分布差异
 
  多维变量X=(x1,x2,...xn)
 X=(x1​,x2​,...xn​)的联合概率密度函数为：
 f(X)=1(2π)d/2∣Σ∣1/2exp[−12(X−u)TΣ−1(X−u)],X=(x1,x2...xn)(2)f(X)=(2π)d/2∣Σ∣1/21​exp[−21​(X−u)TΣ−1(X−u)],X=(x1​,x2​...xn​)(2)
 其中：
 d：变量维度。对于二维高斯分布，有d=2;
 u=⎛⎝⎜⎜⎜u1u2...un⎞⎠⎟⎟⎟u=⎝⎜⎜⎛​u1​u2​...un​​⎠⎟⎟⎞​：各维变量的均值；
 ΣΣ：协方差矩阵，描述各维变量之间的相关度。对于二维高斯分布，有：
 Σ=[δ11δ21δ12δ22](3)
 
Σ=[δ11​δ21​​δ12​δ22​​](3)
 
这里写图片描述
 图2 二维高斯数据分布
 
  图2是二维高斯分布产生的数据示例，参数设定为：u=(00),Σ=[10.80.85]
 
u=(00​),Σ=[10.8​0.85​]。关于二维高斯分布的参数设定对为高斯曲面的影响，这篇文章二维高斯分布（Two-dimensional Gaussian distribution）的参数分析有提及，主要是为了下文理解混合高斯分布做铺垫。服从二维高斯分布的数据主要集中在一个椭圆内部，服从三维的数据集中在一个椭球内部。
 2 混合高斯模型（Gaussian mixture model, GMM）
 2.1 为什么要有混合高斯模型
 
先来看一组数据。
 

 这里写图片描述
 图3 混合高斯分布所产生数据
 
  如果我们假设这组数据是由某个高斯分布产生的，利用极大似然估计（后文还会提及）对这个高斯分布做参数估计，得到一个最佳的高斯分布模型如下。
 
这里写图片描述
 图4 用单高斯模型对样本作分析不合理示意
 
  有什么问题吗？一般来讲越靠近椭圆的中心样本出现的概率越大，这是由概率密度函数决定的，但是这个高斯分布的椭圆中心的样本量却极少。显然样本服从单高斯分布的假设并不合理。单高斯模型无法产生这样的样本。
 实际上，这是用两个不同的高斯分布模型产生的数据。
 
这里写图片描述
 图5 混合高斯模型对样本作分析示意
 
  正当单高斯模型抓耳挠腮的时候，混合高斯模型就大摇大摆地进场了。它通过求解两个高斯模型，并通过一定的权重将两个高斯模型融合成一个模型，即最终的混合高斯模型。这个混合高斯模型可以产生这样的样本。
 更一般化的描述为：假设混合高斯模型由K个高斯模型组成（即数据包含K个类），则GMM的概率密度函数如下：
 p(x)=∑Kk=1p(k)p(x∣k)=∑Kk=1πkN(x∣uk,Σk)(3)
 p(x)=k=1∑K​p(k)p(x∣k)=k=1∑K​πk​N(x∣uk​,Σk​)(3)
 其中，p(x∣k)=N(x∣uk,Σk)p(x∣k)=N(x∣uk​,Σk​)是第k个高斯模型的概率密度函数，可以看成选定第k个模型后，该模型产生x的概率；p(k)=πkp(k)=πk​是第k个高斯模型的权重，称作选择第k个模型的先验概率，且满足∑Kk=1πk=1
 
k=1∑K​πk​=1。
 所以，混合高斯模型并不是什么新奇的东西，它的本质就是融合几个单高斯模型，来使得模型更加复杂，从而产生更复杂的样本。理论上，如果某个混合高斯模型融合的高斯模型个数足够多，它们之间的权重设定得足够合理，这个混合模型可以拟合任意分布的样本。
 2.2 直观上理解混合高斯模型
 
  下面通过几张图片来帮助理解混合高斯模型。
 首先从简单的一维混合高斯模型说起。
 
这里写图片描述
 图6 一维混合高斯模型
 
  在图6中，y1,y2和y3分别表示三个一维高斯模型，他们的参数设定如图所示。y4表示将三个模型的概率密度函数直接相加，注意的是这并不是一个混合高斯模型，因为不满足∑Kk=1πk=1
 
k=1∑K​πk​=1的条件。而y5和y6分别是由三个相同的高斯模型融合生成的不同混合模型。由此可见，调整权重将极大影响混合模型的概率密度函数曲线。另一方面也可以直观地理解混合高斯模型可以更好地拟合样本的原因：它有更复杂更多变的概率密度函数曲线。理论上，混合高斯模型的概率密度函数曲线可以是任意形状的非线性函数。
 下面再给出一个二维空间3个高斯模型混合的例子。
 
这里写图片描述
 (a) 3个类别高斯分布截面轮廓线
 
这里写图片描述
 (b) 混合高斯分布截面轮廓线
 
这里写图片描述
 © 二维混合高斯分布概率密度函数图
 
图7 二维混合高斯模型
 
  (a) 图表示的是3个高斯模型的截面轮廓图，3个模型的权重系数已在图中注明，由截面轮廓图可知3个模型之间存在很多重叠区域。其实这也正是混合高斯模型所希望的。因为如果它们之间的重叠区域较少，那么生成的混合高斯模型一般较为简单，难以生成较为复杂的样本。
 设定好了3个高斯模型和它们之间的权重系数之后，就可以确定二维混合高斯分布的概率密度函数曲面，如图©所示。图(b)是对于图©概率密度曲面的截面轮廓线。从图7也可以看出，整个混合高斯分布曲面相对比于单高斯分布曲面已经异常复杂。实际上，通过调整混合高斯分布的系数(π,μ,Σ)
 
(π,μ,Σ)，可以使得图©的概率密度曲面去拟合任意的三维曲面，从而采样生成所需要的数据样本。
 3 极大似然估计(Maximum Likehood Estimate, MLE)（最大化对数似然函数）
 ● 最大化对数似然函数（log-likelihood function）的意义
 
  首先直观化地解释一下最大化对数似然函数要解决的是什么问题。
 假设我们采样得到一组样本yt
 yt​,而且我们知道变量Y服从高斯分布（本文只提及高斯分布，其他变量分布模型类似），数学形式表示为Y∼N(μ,Σ)Y∼N(μ,Σ)。采样的样本如图8所示，我们的目的就是找到一个合适的高斯分布（也就是确定高斯分布的参数μ,Σ
 
μ,Σ），使得这个高斯分布能产生这组样本的可能性尽可能大。
 
这里写图片描述
 图8 最大化似然函数的意义
 
  那怎么找到这个合适的高斯分布呢（在图8中的表示就是1~4哪个分布较为合适）？这时候似然函数就闪亮登场了。
 似然函数数学化：设有样本集Y=y1,y2...yN
 Y=y1​,y2​...yN​。p(yn∣μ,Σ)p(yn​∣μ,Σ)是高斯分布的概率分布函数，表示变量Y=ynY=yn​的概率。假设样本的抽样是独立的，那么我们同时抽到这N个样本的概率是抽到每个样本概率的乘积，也就是样本集Y的联合概率。此联合概率即为似然函数：
 L(μ,Σ)=L(y1,y2...yN;μ,Σ)=∏Nn=1p(yn;μ,Σ)(4)
 
L(μ,Σ)=L(y1​,y2​...yN​;μ,Σ)=n=1∏N​p(yn​;μ,Σ)(4)
 对式子(4)进行求导并令导数为0（即最大化似然函数，一般还会先转化为对数似然函数再最大化），所求出的参数就是最佳的高斯分布对应的参数。
 所以最大化似然函数的意义就是：通过使得样本集的联合概率最大来对参数进行估计，从而选择最佳的分布模型。
 对于图8产生的样本用最大化似然函数的方法，最终可以得到序号1对应的高斯分布模型是最佳的模型。
 4 EM算法（最大化Q函数）
 4.1 为什么要有EM算法（EM算法与极大似然估计分别适用于什么问题）
 ● 尝试用极大似然估计的方法来解GMM模型
 
  解GMM模型，实际上就是确定GMM模型的参数(μ,Σ,π)
 (μ,Σ,π)，使得由这组参数确定的GMM模型最有可能产生采样的样本。
 先试试看用极大似然估计的方法来解GMM模型会出现什么样的问题。
 如第3小节所述，要利用极大似然估计求解模型最重要的一步就是求出似然函数，即样本集出现的联合概率。而对于混合高斯模型，如何求解某个样本ytyt​的概率？显然我们得先知道这个样本来源于哪一类高斯模型，然后求这个高斯模型生成这个样本的概率p(yt)p(yt​)。
 但是问题来了：我们只有样本。不知道样本到底来源于哪一类的高斯模型。那么如何求解样本的生成概率p(yt)p(yt​)？
 先引入一个隐变量γγ。它是一个K维二值随机变量，在它的K维取值中只有某个特定的元素γkγk​的取值为1，其它元素的取值为0。实际上，隐变量描述的就是：每一次采样，选择第k个高斯模型的概率，故有：
 p(γk=1)=πk(5)p(γk​=1)=πk​(5)
 当给定了γγ的一个特定的值之后（也就是知道了这个样本从哪一个高斯模型进行采样），可以得到样本y的条件分布是一个高斯分布，满足：
 p(y∣γk=1)=N(y∣μk,Σk)(6)p(y∣γk​=1)=N(y∣μk​,Σk​)(6)
 而实际上，每个样本到底是从这K个高斯模型中哪个模型进行采样的，是都有可能的。故样本y的概率为：
 p(y)=∑γp(γ)p(y∣γ)=∑Kk=1πkN(y∣μk,Σk)(7)p(y)=∑γ​p(γ)p(y∣γ)=k=1∑K​πk​N(y∣μk​,Σk​)(7)
 样本集Y(n个样本点)的联合概率为：
 L(μ,Σ,π)=L(y1,y2...yN;μ,Σ,π)=∏Nn=1p(yn;μ,Σ,π)=∏Nn=1∑Kk=1πkN(yn∣μk,Σk)(8)L(μ,Σ,π)=L(y1​,y2​...yN​;μ,Σ,π)=n=1∏N​p(yn​;μ,Σ,π)=n=1∏N​k=1∑K​πk​N(yn​∣μk​,Σk​)(8)
 对数似然函数表示为：
 lnL(μ,Σ,π)=∑Nn=1ln∑Kk=1πkN(yn∣μk,Σk)(9)lnL(μ,Σ,π)=n=1∑N​lnk=1∑K​πk​N(yn​∣μk​,Σk​)(9)
 好了，然后求导，令导数为0，得到模型参数(μ,Σ,π)
 
(μ,Σ,π)。
 貌似问题已经解决了，喜大普奔。
 然而仔细观察可以发现，对数似然函数里面，对数里面还有求和。实际上没有办法通过求导的方法来求这个对数似然函数的最大值。
 MLE（极大似然估计）略显沮丧。这时候EM算法走过来，安慰着说：兄弟别哭，老哥帮你。
 ● 极大似然估计与EM算法适用问题分析
 
  下面先阐述一下极大似然估计与EM算法分别适用于解决什么样的问题。
 
这里写图片描述
 图9 极大似然估计适用问题
 
这里写图片描述
 图10 EM算法适用问题
 
  如果我们已经清楚了某个变量服从的高斯分布，而且通过采样得到了这个变量的样本数据，想求高斯分布的参数，这时候极大似然估计可以胜任这个任务；而如果我们要求解的是一个混合模型，只知道混合模型中各个类的分布模型（譬如都是高斯分布）和对应的采样数据，而不知道这些采样数据分别来源于哪一类（隐变量），那这时候就可以借鉴EM算法。EM算法可以用于解决数据缺失的参数估计问题（隐变量的存在实际上就是数据缺失问题，缺失了各个样本来源于哪一类的记录）。
 下面将介绍EM算法的两个步骤：E-step(expectation-step，期望步)和M-step(Maximization-step,最大化步);
 4.2 E-step
 
我们现有样本集Y=(y1,y2...yT)
 Y=(y1​,y2​...yT​)，通过隐变量γt,kγt,k​（表示ytyt​这个样本来源于第k个模型）的引入，可以将数据展开成完全数据：
 (yt,γt,1,γt,2...γt,K),t=1,2...T
 
(yt​,γt​,1​,γt,2​...γt,K​),t=1,2...T
 
所谓的完全数据，就是不缺失的数据。只有样本集Y=(y1,y2...yT)
 Y=(y1​,y2​...yT​)的数据是不完整的，存在信息缺失的。若ytyt​由第1类采样而来，则有γt,1=1,γt,2=0...γt,K=0γt​,1​=1,γt,2​=0...γt,K​=0，表示为(yt,1,0,...0)(yt​,1,0,...0)。
 所以要求能采到这组数据的可能性，需要分两步走：①第t个样本由哪一类产生？②如果第t个样本由第k类产生，那么第k类产生第t个样本的概率为多少？
 综合考虑上面两步，有了完全数据的似然函数：
 p(y,γ∣μ,Σ,π)=∏Tt=1p(yt,γt,1,γt,2...γt,K∣μ,Σ,π) =∏Tt=1∏Kk=1(πkN(yt;μk,Σk))γt,k =∏Kk=1π∑Tt=1γt,kk∏Tt=1(N(yt;μk,Σk))γt,k(10)p(y,γ∣μ,Σ,π)=t=1∏T​p(yt​,γt​,1​,γt,2​...γt,K​∣μ,Σ,π) =t=1∏T​k=1∏K​(πk​N(yt​;μk​,Σk​))γt,k​ =k=1∏K​πk∑t=1T​γt,k​​t=1∏T​(N(yt​;μk​,Σk​))γt,k​​(10)
 第1个等号到第2个等号的理解：若ytyt​由第1类采样而来，则有γt,1=1,γt,2=0...γt,K=0γt​,1​=1,γt,2​=0...γt,K​=0，
 p(yt,γt,1,γt,2...γt,K∣μ,Σ,π)=∏Kk=1(πkN(yt;μk,Σk))γt,k =(π1N(yt;μ1,Σ1))γt,1(π2N(yt;μ2,Σ2))γt,2...(πKN(yt;μK,ΣK))γt,K =(π1N(yt;μ1,Σ1))1(π2N(yt;μ2,Σ2))0...(πKN(yt;μK,ΣK))0 =π1N(yt;μ1,Σ1)(11)p(yt​,γt​,1​,γt,2​...γt,K​∣μ,Σ,π)=k=1∏K​(πk​N(yt​;μk​,Σk​))γt,k​ =(π1​N(yt​;μ1​,Σ1​))γt,1​(π2​N(yt​;μ2​,Σ2​))γt,2​...(πK​N(yt​;μK​,ΣK​))γt,K​ =(π1​N(yt​;μ1​,Σ1​))1(π2​N(yt​;μ2​,Σ2​))0...(πK​N(yt​;μK​,ΣK​))0 =π1​N(yt​;μ1​,Σ1​)​(11)
 注意式子(11)与式子(7)的差别：如果求p(yt)p(yt​)则需要考虑ytyt​有可能来源于k个类；而如果求的是p(yt,γt,1,γt,2...γt,K)p(yt​,γt​,1​,γt,2​...γt,K​)则已经限定了ytyt​只会来源于某个类。
 第2个等式到第3个等式的理解：先交换累乘符号。由于πkπk​与t无关，故而可以从内部的累乘符号中提取出来。
 实际上完全数据的似然函数描述的就是采集到这些样本的可能性。
 完全数据的对数似然函数为：
 lnp(y,γ∣μ,Σ,π)=∑Kk=1(∑Tt=1γt,k)lnπk+∑Tt=1γt,k(−ln(2π)−12ln∣Σk∣−12(yt−μt)T(Σk)−1(yt−μt))(12)lnp(y,γ∣μ,Σ,π)=k=1∑K​(t=1∑T​γt,k​)lnπk​+t=1∑T​γt,k​(−ln(2π)−21​ln∣Σk​∣−21​(yt​−μt​)T(Σk​)−1(yt​−μt​))(12)
 这一步应该没啥问题吧。。。注意的是，此处考虑的是二维高斯分布的情况，对应于式子(2)中的d=2。
 我们的目标就是找出一组参数(μ∗,Σ∗,π∗)(μ∗,Σ∗,π∗)使得lnp(y,γ∣μ,Σ,π)lnp(y,γ∣μ,Σ,π)最大。
 那么问题来了：lnp(y,γ∣μ,Σ,π)lnp(y,γ∣μ,Σ,π)中含有隐变量γγ，γγ的存在使得我们没法最大化lnp(y,γ∣μ,Σ,π)lnp(y,γ∣μ,Σ,π) 。如果我们知道了γγ，那么最大化lnp(y,γ∣μ,Σ,π)lnp(y,γ∣μ,Σ,π)就显得水到渠成。
 但是坑爹的就是：我们只有采样数据ytyt​，γγ未知。
 那么怎么办呢？对γγ来一个估计。
 猜想我们给了一组起始参数(μ0,Σ0,π0)(μ0,Σ0,π0)或者优化过的第i次迭代的参数(μi,Σi,πi)(μi,Σi,πi)，也就是说每一个高斯分布的参数我们都有了，γγ做的事不就是决定每个样本由哪一个高斯分布产生的嘛，有了每个高斯分布的参数那我们就可以猜想每个样本最有可能来源于哪个高斯分布没错吧！Done!
 为此我们不最大化lnp(y,γ∣μ,Σ,π)lnp(y,γ∣μ,Σ,π)（也无法最大化它），而是最大化Q函数。Q函数如下：
 Q(μ,Σ,π,μi,Σi,πi)=Eγ[lnp(y,γ∣μ,Σ,π)∣Y,μi,Σi,πi]=Eγ[∑Kk=1(∑Tt=1γt,k∣yt,μi,Σi,πi)lnπk+∑Tt=1(γt,k∣yt,μi,Σi,πi)(−ln(2π)−12ln∣Σk∣−12(yt−μt)T(Σk)−1(yt−μt))]=∑Kk=1(∑Tt=1E(γt,k∣yt,μi,Σi,πi)lnπk+∑Tt=1E(γt,k∣yt,μi,Σi,πi)(−ln(2π)−12ln∣Σk∣−12(yt−μt)T(Σk)−1(yt−μt)))Q(μ,Σ,π,μi,Σi,πi)=Eγ​[lnp(y,γ∣μ,Σ,π)∣Y,μi,Σi,πi]=Eγ​[k=1∑K​(t=1∑T​γt,k​∣yt​,μi,Σi,πi)lnπk​+t=1∑T​(γt,k​∣yt​,μi,Σi,πi)(−ln(2π)−21​ln∣Σk​∣−21​(yt​−μt​)T(Σk​)−1(yt​−μt​))]=k=1∑K​(t=1∑T​E(γt,k​∣yt​,μi,Σi,πi)lnπk​+t=1∑T​E(γt,k​∣yt​,μi,Σi,πi)(−ln(2π)−21​ln∣Σk​∣−21​(yt​−μt​)T(Σk​)−1(yt​−μt​)))​
 其中，E(γt,k∣yt,μi,Σi,πi)E(γt,k​∣yt​,μi,Σi,πi)就是对γγ的估计：
 E(γt,k∣yt,μi,Σi,πi)=p(γt,k=1∣yt,μi,Σi,πi) =p(γt,k=1,yt∣μi,Σi,πi)p(yt) =p(γt,k=1,yt∣μi,Σi,πi)∑Kk=1p(γt,k=1,yt∣μi,Σi,πi) =p(yt∣γt,k=1,μi,Σi,πi)p(γt,k=1∣μi,Σi,πi)∑Kk=1p(yt∣γt,k=1,μi,Σi,πi)p(γt,k=1∣μi,Σi,πi) =πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)(14)E(γt,k​∣yt​,μi,Σi,πi)=p(γt,k​=1∣yt​,μi,Σi,πi) =p(yt​)p(γt,k​=1,yt​∣μi,Σi,πi)​ =∑k=1K​p(γt,k​=1,yt​∣μi,Σi,πi)p(γt,k​=1,yt​∣μi,Σi,πi)​ =∑k=1K​p(yt​∣γt,k​=1,μi,Σi,πi)p(γt,k​=1∣μi,Σi,πi)p(yt​∣γt,k​=1,μi,Σi,πi)p(γt,k​=1∣μi,Σi,πi)​ =∑k=1K​πki​N(yt​;μki​,Σki​)πki​N(yt​;μki​,Σki​)​​(14)
 这公式是不是很可怕？？别急，带上几点声明，再去看公式就很好理解了！
 ① Q函数描述的其实就是在给定(μi,Σi,πi)(μi,Σi,πi)参数下，先对样本Y做一个最有可能的划分（每个样本来源于各个类的可能性，即对γγ的估计E(γt,k∣yt,μi,Σi,πi)E(γt,k​∣yt​,μi,Σi,πi)），再描述能够产生这组样本的可能性（Q函数）；
 ② 有了对于γγ的估计之后，Q函数只和样本有关（传统意义上的似然函数亦如此，完全数据的似然函数还与γ
 
γ有关），而不再含有隐变量，从而使得最大化Q函数成为可能；
 ③ 最大化Q函数的过程实则就是使得能够产生这组样本的可能性最大，与最大化似然函数的思路如出一辙。
 4.3 M-step
 
有个Q函数，就可以对Q函数进行最大化，得到下一次迭代的模型参数了，即：
 μi+1,Σi+1,πi+1=argmaxQ(μ,Σ,π,μi,Σi,πi)(15)
 μi+1,Σi+1,πi+1=argmaxQ(μ,Σ,π,μi,Σi,πi)(15)
 对Q函数进行求导，并另其导数为0，可得：
 μi+1k=∑Tt=1πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)ytE(γt,k∣yt,μi,Σi,πi),k=1,2...K(16)μki+1​=E(γt,k​∣yt​,μi,Σi,πi)∑t=1T​∑k=1K​πki​N(yt​;μki​,Σki​)πki​N(yt​;μki​,Σki​)​yt​​,k=1,2...K(16)
 Σi+1k=∑Tt=1πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)(yt−μik)2E(γt,k∣yt,μi,Σi,πi),k=1,2...K(17)Σki+1​=E(γt,k​∣yt​,μi,Σi,πi)∑t=1T​∑k=1K​πki​N(yt​;μki​,Σki​)πki​N(yt​;μki​,Σki​)​(yt​−μki​)2​,k=1,2...K(17)
 πi+1k=E(γt,k∣yt,μi,Σi,πi)T,k=1,2...K(18)πki+1​=TE(γt,k​∣yt​,μi,Σi,πi)​,k=1,2...K(18)
 其中μi+1k,Σi+1k,πi+1k
 
μki+1​,Σki+1​,πki+1​分别表示第（i+1）次迭代，第k个类的均值，协方差矩阵和所占的权重。
 4.4 一个例子梳理EM算法的整个过程
 
  EM算法的核心思想是：通过迭代的过程来找到一组最优的参数(μ∗,Σ∗,π∗)
 
(μ∗,Σ∗,π∗)，使得这组参数表示的模型最有可能产生现有的采样数据。每次迭代的过程就是参数矫正的过程。
 
这里写图片描述
 图11 EM算法参数优化过程
 
  现假设初始化一组参数(μ0,Σ0,π0)
 (μ0,Σ0,π0)。在这组参数下，2类二维高斯分布如图11绿色椭圆所示。然后利用现有的参数，E-step开始对样本数据进行划分（对γγ进行估计）。蓝色的样本大多都被划分给第1类模型，橘黄色的样本大多都被划分给第2类模型。但是第1类模型还有优化空间：第1类模型还不能使得蓝色样本出现的联合概率达到最大。第2类模型也是如此。M-step便优化了2类模型的参数，得到新的参数(μ1,Σ1,π1)(μ1,Σ1,π1)，使得优化后2类高斯分布如图11红色椭圆所示。其中，第1类模型主要优化的是模型均值（即椭圆的中心），第二类模型主要优化的是模型协方差矩阵（即椭圆的长轴、短轴和长短轴的方向）。然后重复进行E-step和M-step，直到参数(μ,Σ,π)(μ,Σ,π)收敛。
 最后谈谈混合高斯模型的参数ππ。
 混合高斯模型的参数μ,Σμ,Σ比较好理解，用于描述各个高斯分布的形状，对于它们的调整也比较直观：使得本高斯分布能够更好地接纳被划分到这类分布的样本。而为什么要有参数ππ？它描述的是各个高斯分布所占的比重，如果不加“歧视”的话（样本来源于各个高斯分布的可能性一致），则有πk=1/Kπk​=1/K；而如果对于某一类高斯分布（即为i）有侧重的话，则相应的πiπi​较大，体现在图11中就是被分配给各个类的样本数占样本总数的比例。如果一轮优化后，某一类高斯分布又接纳了更多样本，则其πiπi​变大，反之变小（所以图11从绿色椭圆调整为红色椭圆实际上两个类所对应的权重也被优化了）。
 而从本质上来看参数ππ，则是为了混合高斯模型能有更好的曲面拟合能力。当参数π
 
π退化为某一类高斯分布的权重远远大于其他类高斯分布的时候，混合高斯模型就退化成了单高斯模型！
 5 总结
 
  图12和图13梳理了高斯分布和混合高斯分布参数估计的逻辑流程。
 
这里写图片描述
 图12 高斯分布参数估计逻辑流程
 
这里写图片描述
 图13 混合高斯分布参数估计逻辑流程
 
  相对比于高斯分布的参数估计，混合高斯分布的参数估计更加复杂。主要原因在于隐变量的存在。而为什么混合高斯分布的参数估计需要多次迭代循环进行？是因为EM算法中对于γ
 γ的估计利用的是初始化或者第i步迭代的参数(μi,Σi,πi)(μi,Σi,πi)，这对于样本的分类划分是有误差的。所以它只能通过多次迭代优化寻找更佳的参数来抵消这一误差。
 终于把这篇文章梳理完了。世界杯要结束了，伪球迷也想见证一下冠军诞生。至此，本文结束。
 ---------------------
 作者：林立民爱洗澡
 来源：CSDN
 原文：https://blog.csdn.net/lin_limin/article/details/81048411
 版权声明：本文为博主原创文章，转载请附上博文链接！