高斯混合模型GMM及EM算法的求解

高斯混合模型GMM与EM算法的求解

一、基础知识

1.1 高斯分布与概率密度函数

二、高斯混合模型GMM介绍

2.1 示例

2.2 高斯混合模型

三、EM算法

3.1 EM算法估计参数

3.2 EM算法流程

一、基础知识

1.1 高斯分布与概率密度函数

高斯分布也就是我们所熟知的正态分布，是研究表明，在数学、物理科学和经济学等学科中，大量数据的分布通常是服从高斯分布， 应用范围非常广泛。高斯分布的定义为：假设随机变量$X$服从高斯分布，即：
$$\begin{array}{c}\mathcal{X}\sim \mathcal{N}(\mu ,{\sigma }^2)\end{array}$$
高斯分布的概率密度函数为：
$$\begin{array}{c}f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}\end{array}$$
其中$\mu$为数学期望，${\sigma }^2$为方差。
如图1所示，均值与方差决定分布的位置与幅度。
图1 高斯分布示例图

对于高维高斯模型，与一维类似，自变量变成了多维，多维变量$\begin{array}{c}X=（x_1,x_2,\dots x_n）\end{array}$,联合概率密度函数为：
$$\begin{array}{l}f(x)=\frac{1}{(2\pi {)}^{\frac{D}{2}}\sqrt{|\Sigma |}}\mid \exp \left(-\frac{(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}{2}\right)\end{array}$$
其中$D$表示数据维度，$\Sigma$表示协方差矩阵，描述各维向量之间的相关度。

二、高斯混合模型GMM介绍

2.1 示例

举一个简单的例子，如图2所示， 对高斯分布做参数估计，得到一个高斯分布模型，很明显这样的分布不太合理，数据点明显分成两个聚类，按照这个模型，椭圆中心的样本数量极少，一般来讲越靠近椭圆中心样本出现的概率越大，这是由概率密度函数所决定的，所以，样本数据服从单高斯分布不合理。

图2 单高斯模型分析样本

这时求解两个高斯模型，通过一定的权重将两个高斯模型融合成一个模型，即最终的混合高斯模型，图3是混合高斯模型对样本进行描述。

图3 混合高斯模型分析样本

2.2 高斯混合模型

高斯混合模型GMM（Gaussian Mixture Model）是一类聚类算法。它是多个高斯分布函数的线性组合，通常用于解决统一集合下的数据包含多种不同分布的情况。
设随机变量$X$，则高斯混合模型的可以描述为：
$$\begin{array}{c}P(x)={\sum }_{k=1}^K{\pi }_{k}N(x|{\mu }_k,{\Sigma }_k)\end{array}$$
其中，$K$为混合模型中子高斯模型的数量，${\pi }_k$ 是混合系数，表示观测数据属于第$k$个子模型的概率，且满足$\begin{array}{c}{\pi }_k\ge 0,{\sum }_{k=1}^K{\pi }_k=1\end{array}，$也可以理解为每个成分的权重， $\begin{array}{c}\mathcal{N}（x|{\mu }_k,{\Sigma }_k）\end{array}$是第$k$个子模型的高斯分布密度函数，$\mu$ 为均值， $\Sigma$是协方差矩阵。

引入一个隐变量$y_j$，它是一个$K$维向量，每一次采样，选择第$k$个高斯模型的概率，它只在某个特定元素取值为1，其它元素取值为0。将高斯混合模型改写为：$$\begin{array}{c}P(x)={\sum }_{k=1}^{K}p(y_j)p(x|y_j)\end{array}$$
对比可得：${\pi }_k=p(y_j)$，$p(x|y_j)=N(x|{\mu }_k,{\Sigma }_k)$,利用这些公式根据贝叶斯理论求解后验概率$p(y_j|x)$。
记$\begin{array}{c}{{\rm Y}}_{jk}=p(y_j|x)\end{array}$，且满足$\begin{array}{c}{\sum }_{k=1}^K{{\rm Y}}_{jk}=1\end{array}$。

三、EM算法

3.1 EM算法估计参数

解GMM模型，实际上就是确定GMM模型的参数$(\mu ,\pi ,\Sigma )$， 由这些参数确定GMM模型最有可能产生的样本。求解参数的流程为：首先写出似然函数，对似然函数取对数，接着对其求导，并令导数为0，得出模型参数。
最大似然函数为：
$$\begin{array}{c}\begin{array}{l}LL={\sum }_{j=1}^M\log \left({\sum }_{k=1}^K{\pi }_{k}N\left(x\mid {\mu }_k,{\Sigma }_k\right)\right)\\ ={\sum }_{j=1}^M\log \left({\sum }_{k=1}^K\frac{{{\rm Y}}_{jk}{\pi }_{k}N\left(x\mid {\mu }_k,{\Sigma }_k\right)}{{{\rm Y}}_{jk}}\right)\\ \ge {\sum }_{j=1}^M{\sum }_{k=1}^K{{\rm Y}}_{jk}\log \left(\frac{{\pi }_{k}N\left(x\mid {\mu }_k,{\Sigma }_k\right)}{{{\rm Y}}_{jk}}\right)\end{array}\end{array}$$

记$\begin{array}{c}H={\sum }_{j=1}^M{\sum }_{k=1}^K{{\rm Y}}_{jk}log(\frac{{\pi }_{k}N(x_n|{\mu }_k,{\Sigma }_k)}{{{\rm Y}}_{jk}})\end{array}$，分别对$\begin{array}{c}{\mu }_k、{\Sigma }_k和{\pi }_k\end{array}$求导，并令结果等于0，得出结果。
$$\begin{array}{}\end{array}\begin{array}{l}H={\sum }_{j=1}^M{\sum }_{k=1}^K{{\rm Y}}_{jk}\log \left(\frac{{\pi }_{k}N\left(x_n\mid {\mu }_k,{\Sigma }_k\right)}{{{\rm Y}}_{jk}}\right)\\ ={\sum }_{j=1}^M{\sum }_{k=1}^K{{\rm Y}}_{jk}\log [\frac{{\pi }_k}{{{\rm Y}}_{jk}}\cdot \frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}\cdot \exp \left(-\frac{1}{2}{\left(x_n-{\mu }_k\right)}^T{\Sigma }_k^{-1}\left(x_n-{\mu }_k\right)\right]\end{array}$$
对${\Sigma }_k$求导
$$\begin{array}{l}\frac{\partial H}{\partial {\Sigma }_k}={\sum }_{j=1}^M{{\rm Y}}_{jk}\left[\frac{{\pi }_k}{{{\rm Y}}_{jk}}\cdot \frac{1}{-{|\sum |}^{\frac{1}{2}}}\cdot \frac{1}{2}{\Sigma }^{-\frac{1}{2}}\left|{\sum }_k\right|{\left({\Sigma }^{-1}\right)}^T+\frac{1}{2}{\left({\Sigma }_k^{-1}\right)}^T{\left(x_n-{\mu }_k\right)}^T\left(x_n-{\mu }_k\right){\left({\Sigma }_k^{-1}\right)}^T\right]\\ ={\sum }_{j=1}^M{{\rm Y}}_{jk}\left[-\frac{1}{2}{\left({\Sigma }_k^{-1}\right)}^T+\frac{1}{2}{\left({\Sigma }_k^{-1}\right)}^T\left(x_n-{\mu }_k\right){\left(x_n-{\mu }_k\right)}^T{\left({\Sigma }_k^{-1}\right)}^T\right]\end{array}$$
令导数为0，即：
$$\begin{array}{l}\frac{\partial H}{\partial {\Sigma }_k}={\sum }_{j=1}^M{{\rm Y}}_{jk}\left[-\frac{1}{2}{\left({\Sigma }_k^{-1}\right)}^T+\frac{1}{2}{\left({\Sigma }_k^{-1}\right)}^T\left(x_n-{\mu }_k\right){\left(x_n-{\mu }_k\right)}^T{\left({\Sigma }_k^{-1}\right)}^T\right]\end{array}=0$$
先左乘${\Sigma }_k^T$，再右乘${\Sigma }_k^T$，注意$\begin{array}{l}{\Sigma }_k={\Sigma }_k^T\end{array}$。

$$\begin{array}{l}{\Sigma }_{j=1}^M{{\rm Y}}_{jk}[-{\Sigma }_k+(x_n-{\mu }_k)(x_n-{\mu }_k{)}^T]\end{array}=0$$
$$\begin{array}{l}-{\Sigma }_{j=1}^M{{\rm Y}}_{jk}{\Sigma }_k+{\Sigma }_{j=1}^M{{\rm Y}}_{jk}(x_n-{\mu }_k)(x_n-{\mu }_k{)}^T\end{array}=0$$
得出：
$${\Sigma }_k=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}\left(x_n-{\mu }_k\right){\left(x_n-{\mu }_k\right)}^T}{{\sum }_{j=1}^M{{\rm Y}}_{jk}}$$
同理，求${\mu }_k$
$$\begin{array}{l}\frac{\partial H}{\partial {\mu }_k}={\sum }_{j=1}^M{{\rm Y}}_{jk}\left[\frac{1}{2}\left({\left({\Sigma }_k^{-1}\right)}^T+{\Sigma }_k^{-1}\right)\left(x_n-{\mu }_k\right)\right]\end{array}=0$$
两边同时左乘${\Sigma }_k$，得到：
$${\mu }_k=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}{x}_n}{{\sum }_{j=1}^M{{\rm Y}}_{jk}}$$
最后，求解${\pi }_k$，${\pi }_k$有限制条件：$\begin{array}{l}{\sum }_{k=1}^K{\pi }_{jk}=1\end{array}$，需要加入拉格朗日算子进行计算：
$$\begin{array}{l}F+\lambda ({\sum }_{k=1}^K{\pi }_k-1)=0\end{array}$$
对${\pi }_k$求导：
$$\begin{array}{l}\frac{\partial F}{\partial {\pi }_k}+\lambda ({\sum }_{k=1}^K{\pi }_k-1)=0\end{array}$$
解出：
$$\begin{array}{l}{\pi }_k=-\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}}{\lambda }\end{array}$$
对和$\lambda$求导：
$$\begin{array}{l}\frac{\partial F}{\partial \lambda }={\sum }_{k=1}^K{\pi }_k-1=0\end{array}$$
将${\pi }_k$值代入上式，其中${\sum }_{k=1}^K{{\rm Y}}_{jk}=1$，${\sum }_{k=1}^K{\sum }_{j=1}^M{{\rm Y}}_{jk}={\sum }_{j=1}^M{\sum }_{k=1}^K{{\rm Y}}_{jk}$，得出：
$$\lambda =-M$$
将$\lambda$代回${\pi }_k$中，得出${\pi }_k$：
$$\begin{array}{l}{\pi }_k=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}}{M}\end{array}$$

3.2 EM算法流程

EM算法是一种迭代算法，具体计算步骤为：
（1）首先初始化参数
（2）E-step，求期望，依据当前的参数，计算每个数据$j$来自模型$k$的可能性，
$$\begin{array}{l}{{\rm Y}}_{jk}=\frac{{\pi }_{k}N\left(x\mid {\mu }_k,{\Sigma }_k\right)}{{\sum }_l{\pi }_{l}N\left(x\mid {\mu }_l,{\Sigma }_l\right)}\end{array}$$
（3）M-step：求极大，计算新一轮迭代的模型参数
$${\mu }_k^{new}=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}{x}_n}{{\sum }_{j=1}^M{{\rm Y}}_{jk}}$$
$${\sum }_k^{new}=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}\left(x_n-{\mu }_k^{new}\right){\left(x_n-{\mu }_k^{new}\right)}^T}{{\sum }_{j=1}^M{{\rm Y}}_{jk}}$$
$$\begin{array}{l}{\pi }_k^{\text{new }}=\frac{{\sum }_{j=1}^M{{\rm Y}}_{jk}}{M}\end{array}$$
（4）重复计算E-step和M-step：直至收敛。

参考:
[1]高斯混合模型
[2]EM算法
[3]GMM 模型与EM算法求解详细推导
[4]高斯混合模型（GMM）及其EM算法的理解