高斯混合模型（GMM）及其求解（期望最大化（EM）算法）

转载请注明出处：http://blog.****.net/u014540876/article/details/79115805

1、高斯混合模型的公式表达

高斯混合模型是指随机变量x具有如下形式的分布（概率密度函数）：
(公式1)
其中，参数θ$\theta$代表所有混合成分的参数（均值向量μ与协方差矩阵Σ）的集合：
(公式2)
每个混合成分的概率密度函数为：
(公式3)
k$k$表示该高斯混合分布由k$k$个高斯混合成分组成。αi${\alpha }_i$为相应高斯混合成分的混合系数，由于高斯混合分布也是一个分布，所以有
(公式4)
即：∑ki=1αi=1(公式5)$\sum _{i=1}^k{\alpha }_i=1(公式5)$

2、高斯混合模型的样本生成过程

高斯混合模型是聚类算法的一种。k−means$k-means$算法是确切给出每个样本被分配到某一个簇，称为硬分配；而高斯混合模型则是给出每个样本被分配到每个簇的概率，最后从中选取一个最大的概率对应的簇作为该样本被分配到的簇，称为软分配。

高斯混合模型是假定按照以下生成过程生成所有样本的（以下过程都是在把αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$当做已知的常数之后来思考）：首先，根据α1,α2,⋅⋅⋅,αk${\alpha }_1,{\alpha }_2,\cdot \cdot \cdot ,{\alpha }_k$定义的先验概率（所谓“先验概率”可以理解为上帝创造的概率，不以其他力量为转移，我们只能通过已经发生的事情结合条件概率公式来推导之，这里的“先验”指的就是这个概率是先于人类能够观察和感受到的经验就已经存在的。）选择高斯混合成分，也就是说αi${\alpha }_i$就是第i个混合成分所占据整个高斯混合成分的比重，这个比重是适应于所有样本的；然后，根据被选择的混合成分的概率密度函数进行采样(样本的生成是指样本从无到有的过程，这里不要认为是先有样本再有高斯混合分布。而我们求解高斯混合分布的过程是由已知的样本倒推高斯混合分布，也就是求解所有的αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$)。

3、高斯混合模型的求解过程

求解高斯混合分布的过程，也就是求解所有的模型参数αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$，此时已知样本集D$D$和预先设置的混合成分的个数k$k$（这种需要预先设置好的参数称为超参数）。

对模型参数进行参数估计，我们都会先想到极大似然估计（高斯（Gauss$Gauss$）最早提出最大似然估计，并将其用于研究测量误差的分布。后来该方法被Fisher$Fisher$总结应用于求参数估计，并命名为极大似然估计），于是，首先写出极大似然函数的对数似然函数：
(公式6)

若训练集D={x1,x2,...,xm}$D=\{x_1,x_2,...,x_m\}$由上述过程生成，令随机变量zj∈{1,2,⋅⋅⋅,k}$z_j\in \{1,2,\cdot \cdot \cdot ,k\}$表示生成样本xj$x_j$的高斯混合成分的序号，其取值未知，显然，zj$z_j$的先验概率P(zj=i)$P(z_j=i)$对应于αi(i=1,2,⋅⋅⋅,k)${\alpha }_i(i=1,2,\cdot \cdot \cdot ,k)$.根据贝叶斯定理，zj$z_j$的后验概率对应于

(公式7)
在此解释一下这个公式。有些人可能在看西瓜书时疑惑于书中这个公式中一会儿用pm$p_m$表示概率，一会儿用p$p$表示概率，其实用p$p$或者pm$p_m$表示都可以，因为它们都是表示的随机变量或者条件随机变量的概率，是等价的。

EM算法的E步

在理想情况下，每个样本应该是只由一个混合成分生成，这个混合成分对应的也就是该样本被分配到的簇，即样本xj只由第i个混合成分生成，亦即p(zj=i∣xj)=1$p(z_j=i\mid x_j)=1$，并且p(zn=i∣xn)=0(当n≠j时)$p(z_n=i\mid x_n)=0(当n\ne j时)$。但是，由于我们的无知，我们预先并不知道这样的理想的高斯混合分布是怎么样的，所以我们只能根据已经观察到的数据集D$D$来获得每个样本由每个混合成分生成的概率，这个概率也就是这个公式所表达的值，所以这也是这个公式为什么是一个后验概率的表达式（已知样本xj$x_j$的情况下，样本xj由第i个混合成分生成的概率），这也就相当于尽我们所能。如果用数学语言表达就是：记随机变量 hji$h_{ji}$表示样本 xj$x_j$是否由第i个混合成分生成（hji$h_{ji}$也就是EM算法中的隐变量），如果是，则hji$h_{ji}$记为1，否则hji$h_{ji}$记为0，根据这个定义我们知道，hj1,hj2,hj3,⋅⋅⋅,hjk$h_{j1},h_{j2},h_{j3},\cdot \cdot \cdot ,h_{jk}$这k$k$个数中，只由一个数为1（说明xj$x_j$只由唯一的一个混合成分生成），其余所有数都为0，那么随机变量hji$h_{ji}$的分布列是

hji（i=1,2,⋅⋅⋅,k）$h_{ji}（i=1,2,\cdot \cdot \cdot ,k）$	0	1
概率	p(hji=0∣xj)$p(h_{ji}=0\mid x_j)$	p(hji=1∣xj)$p(h_{ji}=1\mid x_j)$

获得了隐变量的分布列，接下来求解隐变量的期望：
(公式8)

通过以上过程，我们就完成了EM算法中的E步（即明确隐变量和求解隐变量的期望），也得到了极大似然函数的对数似然函数.接下来执行M步.

EM算法的M步

所谓M步，也就是最大化对数似然函数，因为我们要求解的是αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$作为自变量变化的过程中，对数似然函数的极大值，所以需要将对数似然函数对αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$分别求偏导，接下来以对μi${\mu }_i$求偏导详细说明公式推导过程.

先对μi${\mu }_i$求偏导，这里会涉及到函数链式求导法则和列向量的实数函数对该列向量求导，后者会在后面加以说明。

其中
这一步用到了以下一个列向量的实数函数对该列向量求导的公式：
若A$A$是n$n$阶方阵，x$x$是n$n$维列向量，则

∂(xTAx)∂x=(A+AT)x$\frac{\mathrm{\partial }(x^{T}Ax)}{\mathrm{\partial }x}=(A+A^T)x$，

特殊地，当A$A$是n$n$阶对称方阵时，A=AT$A=A^T$，即

∂(xTAx)∂x=2Ax$\frac{\mathrm{\partial }(x^{T}Ax)}{\mathrm{\partial }x}=2Ax$

而这里的Σ−1i${\mathrm{\Sigma }}_i^{-1}$恰好是m阶对称方阵。

接下来，因为Σ−1i${\mathrm{\Sigma }}_i^{-1}$是第i$i$个混合成分的协方差矩阵的逆，所以它一定是非奇异矩阵，对于非奇异矩阵，有以下定理：
如果n$n$阶方阵A$A$是可逆矩阵，x$x$是n$n$维列向量，那么方程Ax=0$Ax=0$有且只有平凡解（零解），即x=0$x=0$。
令，
则 (公式9)
记，将其带入公式9得：

（公式10）
类似地，可以得到：
（公式11）
对于求解混合系数αi,由于在极大化对数似然函数的情况下，还需要保证
，
所以这是一个含等式和不等式约束的优化问题（其实这里的对数似然函数LL(D)$LL(D)$是相对于αi${\alpha }_i$的凹函数，对它求极大就相当于对LL(D)$LL(D)$求极小，也就是变成了一个标准的凸优化问题），为了简单起见，这里只需要先考虑等式约束，如果在等式约束条件下就能够取得LL(D)$LL(D)$的极小值并且同时能满足不等式约束，那么等式约束条件下的最优解就是这个约束条件下的最优解。所以接下来构造拉格朗日公式：
（公式12）
对其求αi${\alpha }_i$的偏导并令其偏导为0：
（公式13）
对（公式13）左右两边同时乘以αi${\alpha }_i$,得到：
（公式14）
对 （公式14）左右两边同时对i=1,2,3,⋅⋅⋅,k$i=1,2,3,\cdot \cdot \cdot ,k$求和，得到：
（公式15）
其中∑ki=1∑mj=1γji$\sum _{i=1}^k\sum _{j=1}^m{\gamma }_{ji}$代表的是所有的样本由所有的样本生成的后验概率之和，这个值一定是样本的个数，也就是m$m$。
由（公式15）可以得到λ=−m$\lambda =-m$，再结合（公式14），可以得到：

（公式16）说明每个高斯成分的混合系数等于所有属于该混合成分的样本的平均后验概率。
综上，可以获得高斯混合模型的EM算法求解：
先给定混合成分个数k$k$，再初始化高斯混合模型的模型参数(αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$,其中i=1,2,3,⋅⋅⋅,k$i=1,2,3,\cdot \cdot \cdot ,k$)。再执行E步：计算γji${\gamma }_{ji}$，接下来根据（公式10）、（公式11）和（公式16）来更新计算αi,μi,Σi${\alpha }_i,{\mu }_i,{\mathrm{\Sigma }}_i$，循环执行E步和M步，直至算法满足停止条件（例如已达到最大迭代轮数，或者对数似然函数增长很少或者不再增长等等）
求得所有参数之后，再根据

求得xj$x_j$所分配到的簇，输出即可。

后记，第一次写博客，有些地方写的很啰嗦，推导过程很详细，因此也适合初学者查看。错漏之处热烈欢迎各位批评指正，如需详细讨论可发送邮件至我的邮箱:[email protected]