混合高斯模型

一、标准差椭圆

1.1 标准差椭圆概念

1. 这算法最早是由美国南加州大学（UniversityofSouthern California）社会学教授
   韦尔蒂.利菲弗（D. Welty Lefever）提出。
2. 算法又称为Lefever‘s “StandardDeviationalEllipse”（利菲弗方向性分布）
3. 算法最大的特点，是用来度量一组数据的方向和分布的，生成的结果会输出一个椭圆 。

示例图：
（红色的点是伤寒发病的案例，蓝色的区域，发病的数据方向与长江的流向方向基本一致，而范围较大。）

1.2 标准差椭圆（椭圆心 与 椭圆方向）计算

其实算法很简单，要画出一个椭圆：

1. 确定 圆心。
2. 确定 旋转角度。
3. 确定 XY 轴的长度。

SDEx和SDEy为椭圆的圆心
$$SD{E}_x=\sqrt{\frac{{\sum }_{i=1}^n(x_i-\bar{X}{)}^2}{n}}SD{E}_y=\sqrt{\frac{{\sum }_{i=1}^n(y_i-\bar{Y}{)}^2}{n}}$$
其中，$X_i$和$Y_i$是每个要素的空间位置坐标，$X$和$Y$是算数平均中心

椭圆的方向

• 以X轴为准，正北方（12点方向）为0度，顺时针旋转，计算公式如下：
  $$\begin{array}{rl}\tan \theta & =\frac{A+B}{C}\\ \text{A}& =\sum _{i=1}^n{\tilde{x}}_i^2-\sum _{i=1}^n{\tilde{y}}_i^2\\ \text{B}& =\sqrt{{\left(\sum _{i=1}^n{\tilde{x}}_i^2-\sum _{i=1}^n{\tilde{y}}_i^2\right)}^2+4{\left(\sum _{i=1}^n{\tilde{x}}_i{\tilde{y}}_i\right)}^2}\\ \text{C}& =2\sum _{i=1}^n{\tilde{x}}_i{\tilde{y}}_i\end{array}$$
  其中：${\tilde{x}}_i$ 和 ${\tilde{y}}_i$ 是平均中心和 $xy$坐标的差

确定 XY 轴的长度
$$\begin{array}{rl}{\sigma }_x& =\sqrt{2}\sqrt{\frac{{\sum }_{i=1}^n({\tilde{x}}_i\cos \theta -{\tilde{y}}_i\sin \theta {)}^2}{n}}\\ {\sigma }_y& =\sqrt{2}\sqrt{\frac{{\sum }_{i=1}^n({\tilde{x}}_i\sin \theta +{\tilde{y}}_i\cos \theta {)}^2}{n}}\end{array}$$

使用方向：

1. 可用来在地图上标示一组犯罪行为的分布趋势，并且能够确定该行为与特定要素 （一系列酒吧或餐馆、某条特定街道等）的关系。

2. 在地图上标示地下水井样本的特定污染，可以指示毒素的扩散方式，这在部署应急防灾策略时非常有用。

3. 对各个物种所在区域的椭圆的大小、形状和重叠部分进行比较可以分析与物种入侵或者隔离相关的深入信息。

4. 绘制一段时间内疾病爆发情况的椭圆可用于建立疾病传播的模型。

二、高斯方程

高斯分布即正态分布，其曲线图如下：

二维高斯分布，其曲线图如下：

二维高斯函数的一阶偏导数：
$$\frac{\partial G}{\partial x}=\left(-\frac{1}{2\pi {\sigma }^4}\right)x{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

$$\frac{\partial G}{\partial y}=\left(-\frac{1}{2\pi {\sigma }^4}\right)y{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

二维高斯函数的二阶偏导数:
$$\frac{{\partial }^{2}G}{\partial x^2}=\left(-\frac{1}{2\pi {\sigma }^4}\right)\left(1-\frac{x^2}{{\sigma }^2}\right)e^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

$$\frac{{\partial }^{2}G}{\partial y^2}=\left(-\frac{1}{2\pi {\sigma }^4}\right)\left(1-\frac{y^2}{{\sigma }^2}\right)e^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

$$\frac{{\partial }^{2}G}{\partial x\partial y}=\left(\frac{xy}{2\pi {\sigma }^6}\right)e^{-\frac{x^2+y^2}{2{\sigma }^2}}$$

三、高斯混合模型（Gaussian Mixed Model）

$$\exp (-(x^2+y^2))$$
高斯混合模型指的是多个高斯分布函数的线性组合，理论上GMM可以拟合出任意类型的分布:

• 通常用于解决同一集合下的数据包含多个不同的分布的情况（或者是同一类分布但参数不一样， 或者是不同类型的分布，比如正态分布和伯努利分布）。

• 图1中的点在我们看来明显分成两个聚类。这两个聚类中的点分别通过两个不同的正态分布随机生成而来。

• 但是如果没有GMM，那么只能用一个的二维高斯分布来描述图1中的数据。图中的椭圆即为二倍标准差的正态分布椭圆。这显然不太合理，毕竟肉眼一看就觉得应该把它们分成两类。

使用GMM。如图2，数据在平面上的空间分布和图1一样，这时使用两个二维高斯分布来描述图2中的数据， 分别记为$N({\mu }_1,{\Sigma }_1)和N({\mu }_2,{\Sigma }_2)$. 图中的两个椭圆分别是这两个高斯分布的二倍标准差椭圆。

可以看到使用两个二维高斯分布来描述图中的数据显然更合理。 实际上图中的两个聚类的中的点是通过两个不同的正态分布随机生成而来。如果将两个二维高斯分布$N({\mu }_1,{\Sigma }_1)和N({\mu }_2,{\Sigma }_2)$合成一个二维的分布，那么，就可以用合成后的分布来描述图2中的所有点。
最直观的方法就是对这两个二维高斯分布做线性组合，用线性组合后的分布来描述整个集合中的数据。这就是高斯混合模型（GMM）。

3.1 高斯混合模型（GMM） 概念

GMM 由K个gaussian 组成，每个高斯gaussian由Component 组成 设有随机变量$X$，则混合高斯模型可以用表示（GMM 图 1）：

其中$N(x|{\mu }_k,{\Sigma }_k)$称为混合模型中的第$k$个分量（component）。如前面图2中的例子，有两个聚类，可以用两个二维高斯分布来表示，那么分量数$K$=2. ${\pi }_k$是混合系数（mixture coefficient），且满足（GMM 图 2）。

3.2 GMM 用于聚类分析

如果要从 GMM 的分布中随机地取一个点的话，分为两步：

1. 首先随机地在这 K 个 Component 之中选一个，（每个 Component 被选中的概率实际上就是它的系数${\pi }_k$ ）

2. 再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为已知的问题。

   将 GMM 用于聚类时，假设数据服从混合高斯分布， 那么只要根据数据推出 GMM 的概率分布来就可以了；
   然后 GMM 的 K 个 Component 实际上对应K个 cluster 。
   根据数据来推算概率密度通常被称作 density estimation 。
   特别地，当我已知（或假定）概率密度函数的形式，而要估计其中的参数的过程被称作『参数估计』。

如图2的例子，很明显有两个聚类，可以定义K=2. 那么对应的GMM形式如下:
$$p(x)={\pi }_{1}N(x|{\mu }_1,{\Sigma }_1)+{\pi }_{2}N(x|{\mu }_2,{\Sigma }_2)$$

上式中未知的参数有六个： $({\pi }_1,{\mu }_1,{\Sigma }_1;{\pi }_2,{\mu }_2,{\Sigma }_2)$。

之前提到GMM聚类时分为两步:

第一步是随机地在这K个分量中选一个，每个分量被选中的概率即为混合系数 ${\pi }_k$. 可以设定 ${\pi }_1={\pi }_2=0.5$，
表示每个分量被选中的概率是0.5，即从中抽出一个点，这个点属于第一类的概率和第二类的概率各占一半。
但实际应用中事先指定πk的值是很笨的做法，当问题一般化后，会出现一个问题：
当从图2中的集合随机选取一个点，怎么知道这个点是来自 $N(x|{\mu }_1,{\Sigma }_1)$还是 $N(x|{\mu }_2,{\Sigma }_2)$呢？
换言之怎么根据数据自动确定 ${\pi }_1$和 ${\pi }_2$的值？这就是GMM参数估计的问题。要解决这个问题，可以使用EM算法。
通过EM算法，我们可以迭代计算出GMM中的参数： $({\pi }_k,x_k,{\Sigma }_k)$.

3.3 GMM 似然函数 参数估计

假设我们有$N$个数据点，并服从GMM 分布记作（$p(x)$）需要确定${\pi }_k,{\mu }_k,{\Sigma }_k$这些参数。

找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大。
而这个概率实际上就等于 ${\prod }_{i=1}^{N}p(x_i)$这个乘积称作：似然函数。

通常取对数，把乘积变成和 ${\sum }_{i=1}^{N}logp(x_i)$,求函数最大化，即似然函数最大化（这时的参数最优）。

GMM 似然函数（log-likelihood function）推导 :

$$\sum _{i=1}^{N}log\left\{\sum _{k=1}^K{\pi }_{k}N(x_i|{\mu }_k,{\Sigma }_k)\right\}$$
GMM 中随机选点分成两步：（类似k_means 的两步）

1. 每个 Component 生成的概率（不是每个 Component 被选中的概率）；
   估算数据来自哪个Component（分量/组份）；
   对于每个数据xi来说，它由第k个Component 生成的概率（由EM算法中的概率迭代示例可知）：

   下式中的μ和Σ也是待估计的值，因此采样迭代法：
   在计算γ(zk)时假定μ和Σ已知；需要先验给定μ和Σ。
   

2. EM算法估计GMM参数
   估计每个Component 参数的粗略值 (估计每个组份的参数)
   上面得到的 $\gamma (z_k)$ 由Component k 生成的概率 （或者说该Component在生成这个数据上所做的贡献），
   
   反复迭代上述两步，直到似然函数收敛为止。

GMM 总结

分析中我们可以看到 GMM 和 K-means 的迭代求解法其实非常相似（都可以追溯到 EM 算法），
因此也有和 K-means 同样的问题──并不能保证总是能取到全局最优，如果运气比较差，
取到不好的初始值，就有可能得到很差的结果。
对于 K-means 的情况，我们通常是重复一定次数然后取最好的结果，不过 GMM每一次迭代的计算量比 K-means 要大许多。

更流行的做法是先用 K-means（已经重复并取最优值了）得到一个粗略的结果，然后将其作为初值
（只要将 K- means 所得的 centroids 传入 gmm 函数即可），再用 GMM 进行细致迭代。

3.4 EM算法

1. 定义分量数目K，对每个分量k设置${\pi }_k，{\mu }_k和{\Sigma }_k$的初始值,
2. 然后计算对数似然函数:$$p(x|\pi ,\mu ,\Sigma )=\sum _{k=1}^K{\pi }_{k}N(x|{\mu }_k,{\Sigma }_k)$$
3. E step
         ~~~~~      根据当前的${\pi }_k，{\mu }_k,{\Sigma }_k$计算后验概率$\gamma (z_{nk})$
   $$\gamma (z_{nk})=\frac{{\pi }_{k}N(x_n|{\mu }_n,{\Sigma }_n)}{{\sum }_{j=1}^K{\pi }_{j}N(x_n|{\mu }_j,{\Sigma }_j)}$$
4. M step
         ~~~~~      根据E step中计算的$\gamma (z_{nk})$再计算新的${\pi }_k，{\mu }_k,{\Sigma }_k$
   $$\begin{array}{rl}{\mathbf{\mu }}_k^{new}& =\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk}){\mathbf{x}}_n\\ {\mathbf{\Sigma }}_k^{new}& =\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})({\mathbf{x}}_n-{\mathbf{\mu }}_k^{new})({\mathbf{x}}_n-{\mathbf{\mu }}_k^{new}{)}^T\\ {\pi }_k^{new}& =\frac{N_k}{N}\end{array}$$
   其中：
   $$Nk=\sum _{n=1}^N\gamma (z_{nk})$$
5. 计算2中的对数似然函数，得到：
   $$lnp(x|\pi ,\mu ,\Sigma )=\sum _{n=1}^{N}ln\left\{\sum _{k=1}^K{\pi }_{k}N(x_k|{\mu }_k,{\Sigma }_k)\right\}$$
6. 检查参数是否收敛或对数似然函数是否收敛，若不收敛，则返回第2步。