高斯混合模型(GaussianMixture Model)

0 前言

  高斯混合模型（Gaussian Mixture Model）通常简称GMM，是一种广泛使用的聚类算法，该方法使用了高斯分布作为参数模型，并使用了期望最大（Expectation Maximization，简称EM）算法进行训练。

1 单高斯模型

  首先，当随机变量X属于一维情况下的高斯概率密度函数：
$$P(x;\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp(-\frac{(x-\mu {)}^2}{2{\sigma }^2})\text{\hspace{1em}(1.1)}$$
  其中$\mu$为均值，${\sigma }^2$表示方差，其概率密度函数图像如下：

  当随机变量X是多维特征数据时，高斯分布服从下方概率密度函数：

$$P(x,\mathbf{\mu },\mathbf{\Sigma })=\frac{1}{(2\pi {)}^{\frac{D}{2}}{\Sigma }^{\frac{1}{2}}}exp(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu ))\text{\hspace{1em}(1.2)}$$
  其中，$\mu$均值，$\Sigma$为协方差矩阵，D为样本X的维数。其概率密度函数图像如下：

  对于单高斯模型，如果明确训练样本是属于某个高斯模型，我们可以通过极大似然估计求得该高斯模型的参数。假设有样本集 D = { x 1 , x 2 , . . . , x m } D=\{x_1,x_2,...,x_m\} D={x1​,x2​,...,xm​}，其中$x_i\in R^d$，每个样本点都是独立的。通过概率密度函数，可以得到样本集的似然函数：
$$L(\mu ,\Sigma )=\prod _{i=1}^{m}P(x_i)=\prod _{i=1}^m\frac{1}{(2\pi {)}^{\frac{D}{2}}{\Sigma }^{\frac{1}{2}}}exp(-\frac{1}{2}(x_i-\mu {)}^T{\Sigma }^{-1}(x_i-\mu ))\text{\hspace{1em}(1.3)}$$
由于每个点发生的概率都很小，乘积会变得极其小，不利于计算和观察，因此通常我们用 Maximum Log-Likelihood 来计算:
$$l(\mu ,\Sigma )=logL((\mu ,\Sigma )=\sum _{i=1}^m-\frac{D}{2}log2\pi -\frac{1}{2}log|\Sigma |-\frac{1}{2}(x_i-\mu {)}^T{\Sigma }^{-1}(x_i-\mu )\text{\hspace{1em}(1.4)}$$
对最大对数似然函数关于各个参数求偏导，并令结果等于0，可解得各个参数。
$$\begin{gathered} \frac{\partial l(\mu ,\Sigma )}{\partial \mu }=\sum _{i=1}^m-\frac{1}{2}\ast 2{\Sigma }^{-1}(x_i-\mu )=0(\Sigma 为正定矩阵) \\ m\mu =\sum _{i=1}^m{x}_i\Rightarrow \mu =\frac{1}{m}\sum _{i=1}^m{x}_i\text{\hspace{1em}(1.5)} \end{gathered}$$
对$\Sigma$求偏导，我们先了解下面的矩阵求导知识：

• $tr[ABC]=tr[CAB]=tr[BCA]$
• $x^{T}Ax=tr(x^{T}Ax)=tr(x{x}^{T}A)$
• $\frac{\partial tr(AB)}{\partial A}=B^T$
• $\frac{\partial log|A|}{\partial A}=A^{-T}$

$$\begin{array}{rl}l(\mu ,\Sigma )& =\sum _{i=1}^m-\frac{D}{2}log2\pi -\frac{1}{2}log|\Sigma |-\frac{1}{2}(x_i-\mu {)}^T{\Sigma }^{-1}(x_i-\mu )\\ & =\sum _{i=1}^m-\frac{D}{2}log2\pi +\frac{1}{2}log|{\Sigma }^{-1}|-\frac{1}{2}tr[(x_i-\mu )(x_i-\mu {)}^T{\Sigma }^{-1}]\end{array}\text{\hspace{1em}(1.6)}$$
$$\begin{gathered} \frac{\partial l(\mu ,\Sigma )}{\partial {\Sigma }^{-1}}=\frac{m}{2}\Sigma -\frac{1}{2}\sum _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T=0 \\ \Sigma =\frac{1}{m}\sum _{i=1}^m(x_i-\mu )(x_i-\mu {)}^T\text{\hspace{1em}(1.7)} \end{gathered}$$

2 混合高斯模型

  假设设有随机变量X，混合高斯模型由K个高斯模型组成（即数据包含K个类），则GMM的概率密度函数如下：
$$P(x)=\sum _{i=1}^{K}p(k)p(x|k)=\sum _{i=1}^K{\pi }_{i}N(x;{\mu }_i,{\Sigma }_i)\text{\hspace{1em}(2.1)}$$
  其中，$p(x|k)=N(x;{\mu }_k,{\Sigma }_k)$是第k个高斯模型的概率密度函数，${\pi }_k$第k个高斯模型的权重，称作选择第k个模型的先验概率，且满足${\sum }_{i=1}^K{\pi }_k=1$。
m个样本点的联合概率为：
$$L(\mu ,\Sigma ,\pi )=\prod _{i=1}^{m}p(x_i)=\prod _{i=1}^m\left[\sum _{j=1}^K{\pi }_{j}N(x;{\mu }_j,{\Sigma }_j)\right]\text{\hspace{1em}(2.2)}$$
对数似然函数表示为：
$$l(\theta )=logL(\mu ,\Sigma ,\pi )=\sum _{i=1}^{m}log\sum _{j=1}^K{\pi }_{j}N(x;{\mu }_j,{\Sigma }_j)\text{\hspace{1em}(2.3)}$$
对似然函数求导，令导数为0，发现不能求得参数的解析解，下面我们将使用EM算法求解混合高斯模型的参数。
在介绍推导过程之前，先引入几个公式：

离散随机变量X，p(x)表示X的概率密度分布函数，g(x)表示X的某一函数。
离散变量数学期望：$E[g(X)]={\sum }_{i}g(x_i)p(x_i)$
如果f是凸函数，X是随机变量。由Jensen不等式：$E[f(X)]\ge f(EX)$

• E-step
  引入隐变量$z_{ij}$，将一个样本点表示为：
  $$({\mathbf{x}}_{\mathbf{i}},z_1,z_2,...,z_K),i=1,2,...,m$$
  其中 z j ∈ { 0 , 1 } z_j \in \{0,1\} zj​∈{0,1}，$z_j=1$ ,表示样本$x_i$来自第j个高斯模型，且${\sum }_{j=1}^K{z}_j=1$。
  样本集的似然函数为：
  $$L(\theta )=\prod _{i=1}^{m}p(x_i;\theta )=\prod _{i=1}^m\sum _{j=1}^{K}p(x_i,z_j;\theta )\text{\hspace{1em}(2.4)}$$
  对数似然函数表示为：
  $$\begin{array}{rl}l(\theta )& =logL(\theta )=\sum _{i=1}^{m}log\left(\sum _{j=1}^{K}p(x_i,z_j;\theta )\right)\\ & =\sum _{i=1}^{m}log\sum _{j=1}^K{Q}_i(z_j)\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\end{array}\text{\hspace{1em}(2.5)}$$
  其中，${\sum }_{j=1}^K{Q}_i(z_j)=1,Q_i(z)\ge 0$ 。令：
  $$p(z_{ij})=Q_i(z_j)f(z_{ij})=\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}$$
  由期望公式：$E[f(X)]={\sum }_{i}f(x_i)p(x_i)$，可知
  $$\sum _{j=1}^K{Q}_i(z_j)\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}是\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}的数学期望\text{\hspace{1em}(2.6)}$$
  由于$f(x)=logx$是凹函数，根据jensen不等式：
  $$l(\theta )=\sum _{i=1}^{m}log\left[E\left(\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\right)\right]\ge \sum _{i=1}^{m}E\left[log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\right]\text{\hspace{1em}(2.7)}$$
  根据期望的定义，可将不等式右边展开：
  $$\sum _{i=1}^{m}E\left[log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\right]=\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\text{\hspace{1em}(2.8)}$$
  可以得到似然函数关于参数$\theta$的下界：
  $$l(\theta )\ge \sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}=J(\theta )\text{\hspace{1em}(2.9)}$$
  假设$\theta$已经给定，那么$l(\theta )$的值就取决于$Q_i(z_j)$和$p(x_i,z_j;\theta )$了。我们可以通过调整这两个概率使下界$J(\theta )$不断上升，以逼近$l(\theta )$的真实值，当下界$J(\theta )$与似然函数$l(\theta )$的值相等，然后固定$Q_i(z_j)$，利用极大似然估计调整$\theta$，使下界$J(\theta )$值达到最大值，得到的$\theta$为新的${\theta }^{t+1}$；再固定$\theta$，调整$Q_i(z_j)$…直到收敛到似然函数$l(\theta )$的最大值${\theta }^{\ast }$处。过程示意图如下：

接下来我们需要寻找使不等式等号成立的条件，容易知道，不等式等号成立，就是Jensen不等式等号成立：
$$\sum _{i=1}^{m}log\left[E\left(\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\right)\right]\ge \sum _{i=1}^{m}E\left[log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\right]\text{\hspace{1em}(2.10)}$$
要使上面不等式等号成立，当且仅当：
$$\begin{array}{rl}& \frac{p(x_i,z_j;\theta )}{Q_i(z_j)}=C\\ & \Rightarrow p(x_i,z_j;\theta )=C{Q}_i(z_j)\\ & \Rightarrow \sum _{j=1}^{K}p(x_i,z_j;\theta )=C\sum _{j=1}^K{Q}_i(z_j)\because \sum _{j=1}^K{Q}_i(z_j)=1\\ & \Rightarrow \sum _{j=1}^{K}p(x_i,z_j;\theta )=p(x_i;\theta )=C\\ & \therefore Q_i(z_j)=\frac{p(x_i,z_j;\theta )}{{\sum }_{j=1}^{K}p(x_i,z_j;\theta )}=\frac{p(x_i,z_j;\theta )}{p(x_i;\theta )}=p(z_j|x_i;\theta )\end{array}\text{\hspace{1em}(2.11)}$$

1. M-step
   通过E步，我们得到了每个样本i关于$z_j$的分布$Q_i(z_j)$，即已知$Q_i(z_j)$，此时固定$Q_i(z_j)$，求下界$J(\theta )$最大值：
   $$J(\theta )=\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\text{\hspace{1em}(2.12)}$$
   对下界函数$J(\theta )$求关于$\theta$的偏导数：
   $$\begin{array}{rl}J(\theta )& =\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log\frac{p(x_i,z_j;\theta )}{Q_i(z_j)}\\ & =\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)logp(x_i,z_j;\theta )-Q_i(z_j)log{Q}_i(z_j)\\ & =\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log({\pi }_{j}N({\mu }_j,{\Sigma }_j))-Q_i(z_j)log{Q}_i(z_j)\\ & =\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log\left[{\pi }_j\frac{1}{(2\pi {)}^{\frac{D}{2}}|{\Sigma }_j{|}^{\frac{1}{2}}}exp(-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j))\right]-Q_i(z_j)log{Q}_i(z_j)\\ & =\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)\left[log{\pi }_j-\frac{D}{2}log2\pi -\frac{1}{2}log|{\Sigma }_j|-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j)\right]-Q_i(z_j)log{Q}_i(z_j)\end{array}\text{\hspace{1em}(2.13)}$$
   对${\pi }_p$求偏导，需要使用拉格朗日乘子法：
   $$\begin{array}{rl}\frac{\partial J(\theta )}{\partial {\pi }_p}& =\partial \left[\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)log{\pi }_j+\lambda (1-\sum _{j=1}^K{\pi }_j)\right]/\partial {\pi }_p\\ & =\frac{{\sum }_{i=1}^m{Q}_i(z_p)}{{\pi }_p}-\lambda =0\\ & \Rightarrow \sum _{i=1}^m{Q}_i(z_p)=\lambda {\pi }_p\\ & \Rightarrow \sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)=\lambda \sum _{j=1}^K{\pi }_j\\ & \because \sum _{j=1}^K{Q}_i(z_j)=1,\sum _{j=1}^K{\pi }_j=1\\ & \therefore \lambda =m\Rightarrow {\pi }_p=\frac{{\sum }_{i=1}^m{Q}_i(z_p)}{m}\end{array}\text{\hspace{1em}(2.14)}$$
   对${\mu }_p$求偏导，只看和${\mu }_p$有关的部分：
   $$\begin{array}{rl}\frac{\partial J(\theta )}{\partial {\mu }_p}& =\partial \left[\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)(-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j))\right]/\partial {\mu }_p\\ & =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)\frac{\partial (x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p)}{\partial {\mu }_p}\because \frac{\partial \mathbf{u}\mathbf{v}}{\partial \mathbf{x}}=\frac{\partial {\mathbf{u}}^{\mathbf{T}}\mathbf{v}}{\partial \mathbf{x}}=\frac{\partial \mathbf{u}}{\partial \mathbf{x}}\mathbf{v}+\frac{\partial \mathbf{v}}{\partial \mathbf{x}}\mathbf{u}\\ & =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)\left[\frac{\partial (x_i-{\mu }_p)}{\partial {\mu }_p}{\Sigma }_p^{-1}(x_i-{\mu }_p)+\frac{\partial ({\Sigma }_p^{-1}{x}_i-{\Sigma }_p^{-1}{\mu }_p)}{\partial {\mu }_p}(x_i-{\mu }_p)\right]\\ & =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)\left[\frac{\partial (-{\mu }_p)}{\partial {\mu }_p}{\Sigma }_p^{-1}(x_i-{\mu }_p)+\frac{\partial (-{\Sigma }_p^{-1}{\mu }_p)}{\partial {\mu }_p}(x_i-{\mu }_p)\right]\\ & =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)\left[-{\Sigma }_p^{-1}(x_i-{\mu }_p)-({\Sigma }_p^{-1}{)}^T(x_i-{\mu }_p)\right]\\ & =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)\left[-{\Sigma }_p^{-1}(x_i-{\mu }_p)-{\Sigma }_p^{-1}(x_i-{\mu }_p)\right]\\ & =\sum _{i=1}^m{Q}_i(z_p)[{\Sigma }_p^{-1}(x_i-{\mu }_p)]=0\\ & \Rightarrow {\mu }_p=\frac{{\sum }_{i=1}^m{Q}_i(z_p)x_i}{{\sum }_{i=1}^m{Q}_i(z_p)}\end{array}\text{\hspace{1em}(2.15)}$$
   对${\Sigma }_p$求偏导，只看和${\Sigma }_p$有关的部分：
   $$\begin{array}{rl}\frac{\partial J(\theta )}{\partial {\Sigma }_p}& =\partial \left[\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)\left(-\frac{1}{2}log|{\Sigma }_j|-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j)\right)\right]/\partial {\Sigma }_p\\ & =\partial \left[\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)(log|{\Sigma }_p|+(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p))\right]/\partial {\Sigma }_p\\ & \because \frac{\partial log|{\Sigma }_p|}{\partial {\Sigma }_p}=\frac{1}{|{\Sigma }_p|}\frac{|{\Sigma }_p|}{\partial {\Sigma }_p}=\frac{1}{|{\Sigma }_p|}|{\Sigma }_p|({\Sigma }_p^{-1}{)}^T={\Sigma }_p^{-1}\\ & \because \frac{\partial [(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p)]}{\partial {\Sigma }_p}=\frac{\partial tr[(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p)]}{\partial {\Sigma }_p}=\frac{\partial tr[{\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T]}{\partial {\Sigma }_p}\\ & \frac{\partial tr[{\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T]}{\partial {\Sigma }_p}=-({\Sigma }_p^{-1}{)}^T[(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{]}^T({\Sigma }_p^{-1}{)}^T=-{\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}\\ \frac{\partial J(\theta )}{\partial {\Sigma }_p}& =\sum _{i=1}^m-\frac{1}{2}{Q}_i(z_p)({\Sigma }_p^{-1}-{\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1})=0\\ & \Rightarrow \sum _{i=1}^m{Q}_i(z_p){\Sigma }_p^{-1}=\sum _{i=1}^m{Q}_i(z_p){\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(右边同时乘以{\Sigma }_p)\\ & \Rightarrow \sum _{i=1}^m{Q}_i(z_p)I=\sum _{i=1}^m{Q}_i(z_p){\Sigma }_p^{-1}(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T(左边同时乘以{\Sigma }_p)\\ & \Rightarrow \sum _{i=1}^m{Q}_i(z_p){\Sigma }_p=\sum _{i=1}^m{Q}_i(z_p)(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T\\ & \Rightarrow {\Sigma }_p=\frac{{\sum }_{i=1}^m{Q}_i(z_p)(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T}{{\sum }_{i=1}^m{Q}_i(z_p)}\end{array}\text{\hspace{1em}(2.16)}$$

对${\Sigma }_p$求偏导的另一种方式，对${\Sigma }_p^{-1}求偏导$：
$$\begin{array}{rl}\frac{\partial J(\theta )}{\partial {\Sigma }_p^{-1}}& =\partial \left[\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)\left(-\frac{1}{2}log|{\Sigma }_j|-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j)\right)\right]/\partial {\Sigma }_p^{-1}\\ & =\partial \left[\sum _{i=1}^m\sum _{j=1}^K{Q}_i(z_j)\left(\frac{1}{2}log|{\Sigma }_j^{-1}|-\frac{1}{2}(x_i-{\mu }_j{)}^T{\Sigma }_j^{-1}(x_i-{\mu }_j)\right)\right]/\partial {\Sigma }_p^{-1}\\ & =\partial \left[\sum _{i=1}^m\frac{1}{2}{Q}_i(z_p)(log|{\Sigma }_p^{-1}|-(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p))\right]/\partial {\Sigma }_p^{-1}\\ & =\sum _{i=1}^m\frac{1}{2}{Q}_i(z_p)\left[(({\Sigma }_p^{-1}{)}^{-1}{)}^T-\frac{\partial tr[(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}(x_i-{\mu }_p)]}{\partial {\Sigma }_p^{-1}}\right]\\ & =\sum _{i=1}^m\frac{1}{2}{Q}_i(z_p)\left[{\Sigma }_p-\frac{\partial tr[(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{\Sigma }_p^{-1}]}{\partial {\Sigma }_p^{-1}}\right]\\ & =\sum _{i=1}^m\frac{1}{2}{Q}_i(z_p)\left[{\Sigma }_p-[(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T{]}^T\right]\\ & =\sum _{i=1}^m\frac{1}{2}{Q}_i(z_p)\left[{\Sigma }_p-(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T\right]=0\\ & \Rightarrow {\Sigma }_p=\frac{{\sum }_{i=1}^m{Q}_i(z_p)(x_i-{\mu }_p)(x_i-{\mu }_p{)}^T}{{\sum }_{i=1}^m{Q}_i(z_p)}\end{array}\text{\hspace{1em}(2.17)}$$

3 EM算法

通过上面对EM算法的理论推导，现在我们对EM算法的实现步骤进行一个总结：

1. 初始化参数 K , θ = { π j , μ j , Σ j } j = 1 , 2 , . . , K K,\theta=\{ \pi_j,\mu_j,\Sigma_j \} \quad j=1,2,..,K K,θ={πj​,μj​,Σj​}j=1,2,..,K
2. E-step:根据当前参数${\pi }_j,{\mu }_j,{\Sigma }_j$计算后验概率$Q_i(z_j)$：
   $$Q_i(z_j)=\frac{{\pi }_{j}N(x_i;{\mu }_j,{\Sigma }_j)}{{\sum }_{p=1}^K{\pi }_{p}N(x_i;{\mu }_p,{\Sigma }_p)}\text{\hspace{1em}(3.1)}$$
3. M-step:根据E步中计算的$Q_i(z_j)$，重新计算参数 θ = { π j , μ j , Σ j } j = 1 , 2 , . . , K \theta=\{ \pi_j,\mu_j,\Sigma_j \} \quad j=1,2,..,K θ={πj​,μj​,Σj​}j=1,2,..,K
4. 计算对数似然函数，观察对数似然函数值的变化是否收敛，收敛则结束，否则转到步骤2，直到收敛。

注：通过EM求解混合高斯模型，得到的解不一定是全局最优解，可能收敛到局部最优。

4 代码实现

import numpy as np
import matplotlib.pyplot as plt

#高斯概率密度函数
def gaussian(x, Mean, Cov_matrix):
    dim = np.shape(cov)[0]  # 维度
    # 之所以加入单位矩阵是为了防止行列式为0的情况
    covdet = np.linalg.det(Cov_matrix + np.eye(dim) * 0.01)  # 协方差矩阵的行列式
    covinv = np.linalg.inv(Cov_matrix + np.eye(dim) * 0.01)  # 协方差矩阵的逆
    xdiff = x - Mean
    # 概率密度
    prob = 1.0 / np.power(2 * np.pi, 1.0 * dim / 2) / np.sqrt(np.abs(covdet)) * np.exp(
        -0.5 * xdiff.T @ covinv @ xdiff)
    return prob

def GMM(X, K, iter_num=10):
    global  imgnum
    m, n = X.shape
    # 初始化参数
    oldpi = pi = np.full(K, 1.0 / K)
    oldmu = mu = [X[i] for i in np.roll(np.arange(K), np.random.choice(m))]
    olsigma = sigma = [np.eye(n) for i in range(K)]
    Q = np.zeros((m, K))
    errorlost = 0
    for itn in range(iter_num):
        plt.close()
        plt.plot(X[:, 0], X[:, 1], '.')
        # E_step
        for i in range(m):
            pxz = [pi[j] * gaussian(X[i], mu[j], sigma[j]) for j in range(K)]
            pxzsum = np.sum(pxz)
            Q[i] = pxz / pxzsum
        Qcol = np.sum(Q, axis=0)
        pi = Qcol / m
        mu = [np.sum(Q[:, j:j + 1] * X, axis=0) / Qcol[j] for j in range(K)]
        for j in range(K):
            xdif = X - mu[j]
            sigma[j] = 1.0 / Qcol[j] * np.sum(
                [Q[i][j] * (xdif[i:i + 1].T @ xdif[i:i + 1]) for i in range(m)], axis=0)
        # 画出每个高斯分布的均值
        for j in range(K):
            plt.plot(mu[j][0], mu[j][1], marker='o', c='r')
        plt.savefig(str(imgnum)+".png")
        imgnum += 1
        plt.show()
        curerror = sum(np.power((oldpi - pi).tolist(), 2)) + sum(np.power(np.ravel(oldmu) - np.ravel(mu), 2)) + sum(
            np.power(np.ravel(olsigma) - np.ravel(sigma), 2))
        if abs(curerror - errorlost) < 0.0001:
            break
        errorlost = curerror
        print("第{:d}次迭代参数平方损失{:6f}：".format(itn, errorlost))
    return pi, mu, sigma

if __name__ == '__main__':
    # 初始化参数
    mean = [2, 2]
    cov = [[1, 0], [0, 1]]
    s1 = np.random.multivariate_normal(mean, cov, 80, "raise")
    mean = [10, 8]
    cov = [[1, 0], [0, 4]]
    s2 = np.random.multivariate_normal(mean, cov, 100, "raise")
    mean = [4, 13]
    cov = [[2, 1], [1, 3]]
    s3 = np.random.multivariate_normal(mean, cov, 120, "raise")

    sall = np.vstack((s1, s2, s3))
    np.random.shuffle(sall)
    # 高斯分量的个数
    K = 3
    # 最大迭代次数
    iter_numbers = 60
    pi, mu, sigma = GMM(sall, K, iter_numbers)
    pointnum = 100
    xrange = np.linspace(sall[:, 0].min(), sall[:, 0].max(), pointnum)
    yrange = np.linspace(sall[:, 1].min(), sall[:, 1].max(), pointnum)
    XX, YY = np.meshgrid(xrange, yrange)
    plt.close()
    plt.scatter(s1[:, 0], s1[:, 1], marker='.', c='r')
    plt.scatter(s2[:, 0], s2[:, 1], marker='o', c='b')
    plt.scatter(s3[:, 0], s3[:, 1], marker='x', c='g')
    plt_data = np.dstack((XX, YY))
    for j in range(K):
        Z = [[gaussian(plt_data[q][p], mu[j], sigma[j]) for p in range(pointnum)] for q in range(pointnum)]
        cs = plt.contour(XX, YY, Z)
        plt.clabel(cs)
    plt.savefig(str(imgnum) + ".png")
    plt.show()

5 参考

高斯混合模型的终极理解