请和我一起学习机器学习算法（高斯混合聚类）

1. 高斯混合聚类的原理分析

高斯混合的含义

高斯混合聚类有一个重要的前提：
假定样本集是从若干个（需要聚类的个数）满足高斯分布的集合中按照一定的比例随机抽取生成的。而高斯混合聚类的任务就是上面假定的逆过程，从这些杂乱的高斯分布中抽取出来的样本集，重新归类到集合中。

定类别

既然高斯混合聚类是从若干个满足高斯分布的集合中抽取而成，那么高斯分布是怎么样的呢？根据定义多元的高斯分布的概率密度函数长下面这个样子。
$$p(\mathbf{x})=\frac{1}{(2\pi {)}^{n/2}|\Sigma {|}^{-1}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T{\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu })}$$
这个样子简直太麻烦了，所以我们换成下面这种简单的：
$$p(\mathbf{x}|\mathbf{\mu },\Sigma )$$
表示概率密度值由后面两个变量决定，一个是均值，一个是协方差矩阵。均值大家肯定都理解，协方差矩阵不好理解，先放着，有时间在看。
那接下来就是说混合高斯应该是什么样子的呢？不就是若干个（假设为k个）简单加起来吗，一定概率，则是说每个高斯类别具有不同的权重，简单来说就是多了个系数$a_i$呗。
$$p_M(\mathbf{x})=\sum _{i=1}^k{a}_{i}p(\mathbf{x}|{\mathbf{\mu }}_i,{\Sigma }_i)$$
既然$a_i$表示为概率，则总和为1，表示样本只能从这几个高斯分布中抽取。官方说法叫做 混合系数
$$\sum _{i=i}^k{a}_i=1$$
现在，假定我从样本集中随机选出一个样本${\mathbf{x}}_j$，则这个样本是属于这k个高斯集合中的哪个的呢？很显然，这就是我们聚类要干的事情。虽然不知道是哪个，但是属于任何一个的概率我们是知道的。设这个样本属于$z_j$
$$p(z_j=i)=a_i对于任意的i=1,2,3,...,k$$
为什么？因为样本集我们就是按照$a_i$权重从各个集合中取出的呀。

接下来就是最为关键的推到了，对于一个已知的一个样本${\mathbf{x}}_j$，$p(z_j=i)$又是多少呢？（前面的是随机抽取一个样本，会得到什么结果我们未知.） 这明显是一个条件概率的问题。也就是在条件${\mathbf{x}}_j$下，求取$p(z_j=i)$的值。根据贝叶斯公式$p(x|y)\cdot p(y)=p(y|x)\cdot p(x)$

$$\begin{gathered} p_M(z_j=i|{\mathbf{x}}_j)=\frac{p(z_j=i)p_M({\mathbf{x}}_j|z_j=i)}{p_M({\mathbf{x}}_j)} \\ {p}_M(z_j=i|{\mathbf{x}}_j)=\frac{a_i\cdot p({\mathbf{x}}_j|{\mathbf{x}}_i,{\mathbf{\Sigma }}_i)}{{\sum }_{l=1}^k{a}_{l}p({\mathbf{x}}_j|{\mathbf{\mu }}_l,{\Sigma }_l)} \end{gathered}$$
转化很明显，$p_M({\mathbf{x}}_j|z_j=i)$ 后面的条件表面只属于第i个集合，则概率记为第i个集合的概率，分母则是混合高斯的组成公式的展开。
简单记作：
$${\gamma }_{ji}=p_M(z_j=i|{\mathbf{x}}_j)$$
表示在$x_j$的条件下，$x_j$属于第i个集合的概率，这时一个很关键的概率。
如果混合高斯分布已知，则我们只需要将所有i（1~k）的概率都计算出来，概率最大的即可以认为$x_j$ 属于这个概率最大的高斯分布（假设为第${\lambda }_j个$），则：
λ j = a r g max ⁡ i ∈ { 1 , 2 , . . k } γ j i \lambda_j= arg \max_{i \in \{1,2,..k\}} \gamma_{ji} λj​=argi∈{1,2,..k}max​γji​
讲解到这里，问题已经很明显了：
我们的目的是求出使得参数${\gamma }_{ji}$最大的变量i，作为${\mathbf{x}}_j$的类别标记。
但是要求解参数i， 则需要获得混合高斯分布$p_M(\mathbf{x})$的参数
问题转化为：如何求解混合高斯分布$p_M(\mathbf{x})$的各个参数，从而求解每一个样本对应的$\lambda$（也就是我们所说的标记）

估参数

我们从未知的分布$p_M(\mathbf{x})$(也就是混合高斯分布)中采取了样本集，现在要计算分布中的各个参数，这不就是传说中的参数估计问题吗？由于我们要用这些参数去计算概率，所以必然是点估计，点估计中与条件有关的，最好使的不就是极大似然估计吗？
要使用极大似然估计，必然要先搞清楚，我们需要估计的参数是什么对吧。是时候把我们需要估计的模型（也就是刚刚的概率分布公式啦）拿出来了。
$$p_M(\mathbf{x})=\sum _{i=1}^k{a}_{i}p(\mathbf{x}|{\mathbf{\mu }}_i,{\Sigma }_i)$$
明显的，我们需要估计{$a_i,{\mu }_i,{\Sigma }_i$|i=1,2,…k} 也就是k组这三个量。
接下来计算极大似然函数（取对数形式），对于含由m个样本的集合D，
$$\begin{array}{rl}LL(D)& =\ln (\prod _{j=1}^m{p}_M({\mathbf{x}}_j))\\ & =\sum _{j=1}^m\ln (\sum _{i=1}^k{a}_{i}p({\mathbf{x}}_j|{\mathbf{\mu }}_i,{\Sigma }_i))\end{array}$$
也就说在已知{$m,k,\mathbf{x}$}的情况下，计算使得$LL(D)$最大的参数{$a_i,{\mu }_i,{\Sigma }_i|1\le i\le k$} 。
假设我们有一组参数{$a_i,{\mu }_i,{\Sigma }_i|1\le i\le k$} 使得值最大，则：
$$\begin{array}{rl}\frac{\partial LL(D)}{\partial {\mathbf{\mu }}_i}& =0\\ & =\sum _{j=1}^m\frac{a_i\cdot p({\mathbf{x}}_j|{\mathbf{\mu }}_i,{\Sigma }_i)}{{\sum }_{i=1}^k{a}_{i}p({\mathbf{x}}_j|{\mathbf{\mu }}_i,{\Sigma }_i)}({\mathbf{x}}_j-{\mathbf{\mu }}_i)\\ & =\sum _{j=1}^m{\gamma }_{ij}({\mathbf{x}}_j-{\mathbf{\mu }}_i)=0\end{array}$$
所以有
$${\mathbf{\mu }}_i=\frac{{\sum }_{j=1}^m{\gamma }_{ji}{\mathbf{x}}_j}{{\sum }_{j=1}^m{\gamma }_{ji}}$$
类似的，根据
$$\begin{array}{rl}\frac{\partial LL(D)}{\partial {\mathbf{\Sigma }}_i}& =0\end{array}$$
可以得
【实际上，这些公式的推导都十分麻烦，包含了大量的矩阵内容，我还没搞懂，后面搞懂了补上】
$${\mathbf{\Sigma }}_i=\frac{{\sum }_{j=1}^m{\gamma }_{ji}({\mathbf{x}}_j-{\mathbf{\mu }}_i)({\mathbf{x}}_j-{\mathbf{\mu }}_i{)}^T}{{\sum }_{j=1}^m{\gamma }_{ji}}$$
对于参数a，有条件$a_i\ge 0$，并且${\sum }_{i=1}^m{a}_i=1,$
$$\begin{array}{rl}\frac{\partial (LL(D)+({\sum }_{i=1}^m{a}_i-1))}{\partial a_i}& =0\\ & =\sum _{j=1}^m\frac{p({\mathbf{x}}_j|{\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)}{{\sum }_{l=1}^k{a}_l\ast p({\mathbf{x}}_j|{\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)}+\lambda =0\end{array}$$
两边同时乘以$a_i$, 对所有样本求和可知，
$$\begin{array}{rl}\lambda & =-m\\ a_i& =\frac{1}{m}\sum _{j=1}^m{\gamma }_{ji}\end{array}$$

分类估参联合求解(EM)

接下来就是EM算法了，

1. 先初始化当前参数。
2. 根据当前参数计算${\lambda }_{ji}$,【E步】
   $${\gamma }_{ji}=p_M(z_j=i|{\mathbf{x}}_j)$$
3. 根据计算出来的${\lambda }_{ji}$，反推参数。【M步】
   $$\begin{array}{rl}{\mathbf{\mu }}_i& =\frac{{\sum }_{j=1}^m{\gamma }_{ji}{\mathbf{x}}_j}{{\sum }_{j=1}^m{\gamma }_{ji}}\\ {\mathbf{\Sigma }}_i& =\frac{{\sum }_{j=1}^m{\gamma }_{ji}({\mathbf{x}}_j-{\mathbf{\mu }}_i)({\mathbf{x}}_j-{\mathbf{\mu }}_i{)}^T}{{\sum }_{j=1}^m{\gamma }_{ji}}\\ a_i& =\frac{1}{m}\sum _{j=1}^m{\gamma }_{ji}\end{array}$$
4. 重复2-3步，达到EM终止条件（比如：达到最大迭代次数，似然函数$LL(D)$增长很少，等）。
5. 算法结束

2. 高斯混合聚类的代码实现

python sklearn 库中有现成的算法：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.mixture import GaussianMixture as GMM
from sklearn.datasets.samples_generator import make_blobs

if __name__ == "__main__":
    #Generate some data for train
    X, y_true = make_blobs(n_samples=400, centers=4,
                       cluster_std=0.60, random_state=0)
    X = X[:, ::-1] # flip axes for better plotting
    #plt.figure(1)
    ax=plt.subplot(1, 2, 1) 
    ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)
  
    gmm = GMM(n_components=4, random_state=42)
    labels = gmm.fit(X).predict(X)

    bx=plt.subplot(1, 2, 2)
    bx.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
    plt.show()

代码引用https://jakevdp.github.io/PythonDataScienceHandbook/05.12-gaussian-mixtures.html
原理部分主要参考西瓜书
实际上是有一点粗糙，后续代码实战过程中手撕它。