本文翻译自 https://jakevdp.github.io/PythonDataScienceHandbook/05.12-gaussian-mixtures.html

上一节中探讨的k-means聚类模型简单易懂，但其简单性导致其应用中存在实际挑战。具体而言，k-means的非概率特性及简单地计算点与类蔟中心的欧式距离来判定归属，会导致其在许多真实的场景中性能较差。本节，我们将探讨高斯混合模型(GMMs)，其可以看成k-means的延伸，更可以看成一个强有力的估计工具，而不仅仅是聚类。

我们将以一个标准的import开始

我们看下k-means的缺陷，思考下如何提高聚类模型。正如上一节所示，给定简单，易于分类的数据，k-means能找到合适的聚类结果。

举例而言，假设我们有些简单的数据点，k-means算法能以某种方式很快地将它们聚类，跟我们肉眼分辨的结果很接近：

从直观的角度来看，我可能期望聚类分配时，某些点比其他的更确定：举例而言，中间两个聚类之间似乎存在非常轻微的重叠，这样我们可能对这些数据点的分配没有完全的信心。不幸的是，k-means模型没有聚类分配的概率或不确定性的内在度量（尽管可能使用bootstrap 的方式来估计这种不确定性）。为此，我们必须考虑泛化这种模型。

k-means模型的一种理解思路是，它在每个类蔟的中心放置了一个圈（或者，更高维度超球面)，其半径由聚类中最远的点确定。该半径充当训练集中聚类分配的一个硬截断：任何圈外的数据点不被视为该类的成员。我们可以使用以下函数可视化这个聚类模型：

观察k-means的一个重要发现，这些聚类模式必须是圆形的。k-means没有内置的方法来计算椭圆形或椭圆形的簇。因此，举例而言，假设我们将相同的数据点作变换，这种聚类分配方式最终变得混乱：

高斯混合模型（GMM）试图找到一个多维高斯概率分布的混合，以模拟任何输入数据集。在最简单的情况下，GMM可用于以与k-means相同的方式聚类。

但因为GMM包含概率模型，因此可以找到聚类分配的概率方式 - 在Scikit-Learn中，通过调用predict_proba方法实现。它将返回一个大小为[n_samples, n_clusters]的矩阵，用于衡量每个点属于给定类别的概率：

我们可以可视化这种不确定性，比如每个点的大小与预测的确定性成比例；如下图，我们可以看到正是群集之间边界处的点反映了群集分配的不确定性：

本质上说，高斯混合模型与k-means非常相似：它使用期望-最大化的方式，定性地执行以下操作：

有了这个，我们可以看看四成分的GMM为我们的初始数据提供了什么：

同样，我们可以使用GMM方法来拟合我们的拉伸数据集；允许full的协方差，该模型甚至可以适应非常椭圆形，伸展的聚类模式：

这清楚地表明GMM解决了以前遇到的k-means的两个主要实际问题。

如果看了之前拟合的细节，你将看到covariance_type选项在每个中都设置不同。该超参数控制每个类簇的形状的自由度；对于任意给定的问题，必须仔细设置。默认值为covariance_type =“diag”，这意味着可以独立设置沿每个维度的类蔟大小，并将得到的椭圆约束为与轴对齐。一个稍微简单和快速的模型是covariance_type =“spherical”，它约束了类簇的形状，使得所有维度都相等。尽管它并不完全等效，其产生的聚类将具有与k均值相似的特征。更复杂且计算量更大的模型（特别是随着维数的增长）是使用covariance_type =“full”，这允许将每个簇建模为具有任意方向的椭圆。

对于一个类蔟，下图我们可以看到这三个选项的可视化表示：

尽管GMM通常被归类为聚类算法，但从根本上说它是一种密度估算算法。也就是说，GMM适合某些数据的结果在技术上不是聚类模型，而是描述数据分布的生成概率模型。

例如，考虑一下Scikit-Learn的make_moons函数生成的一些数据：

如果我们尝试用视为聚类模型的双成分的GMM模拟数据，则结果不是特别有用：

但是如果我们使用更多成分的GMM模型，并忽视聚类的类别，我们会发现更接近输入数据的拟合：

这里，16个高斯分布的混合不是为了找到分离的数据簇，而是为了对输入数据的整体分布进行建模。这是分布的一个生成模型，这意味着GMM为我们提供了生成与我们的输入类似分布的新随机数据的方法。例如，以下是从这个16分量GMM拟合到我们原始数据的400个新点：

GMM非常方便，可以灵活地建模任意多维数据分布。

GMM是一种生成模型这一事实为我们提供了一种确定给定数据集的最佳组件数的自然方法。生成模型本质上是数据集的概率分布，因此我们可以简单地评估模型下数据的可能性，使用交叉验证来避免过度拟合。校正过度拟合的另一种方法是使用一些分析标准来调整模型可能性，例如 Akaike information criterion (AIC) 或 Bayesian information criterion (BIC) 。Scikit-Learn的GMM估计器实际上包含计算这两者的内置方法，因此在这种方法上操作非常容易。

让我们看看在moon数据集中，使用AIC和BIC函数确定GMM组件数量：

最佳的聚类数目是使得AIC或BIC最小化的值，具体取决于我们希望使用的近似值。 AIC告诉我们，我们上面选择的16个组件可能太多了：大约8-12个组件可能是更好的选择。与此类问题一样，BIC建议使用更简单的模型。

注意重点：这个组件数量的选择衡量GMM作为密度估算器的效果，而不是它作为聚类算法的效果。我鼓励您将GMM主要视为密度估算器，并且只有在简单数据集中保证时才将其用于聚类。

我们刚刚看到了一个使用GMM作为数据生成模型的简单示例，以便根据输入数据定义的分布创建新样本。在这里，我们将运行这个想法，并从我们以前使用过的标准数字语料库中生成新的手写数字。

首先，让我们使用Scikit-Learn的数据工具加载数字数据：

接下来让我们绘制前100个，以准确回忆我们正在看的内容：

我们有64个维度的近1,800位数字，我们可以在这些位置上构建GMM以产生更多。 GMM可能难以在如此高维空间中收敛，因此我们将从数据上的可逆维数减少算法开始。在这里，我们将使用一个简单的PCA，要求它保留99％的预测数据方差：

结果是41个维度，减少了近1/3，几乎没有信息丢失。根据这些预测数据，让我们使用AIC来计算我们应该使用的GMM组件的数量：

似乎大约110个components最小化了AIC；我们将使用这个模型。我们迅速将其与数据拟合并确保它已收敛合：

现在我们可以使用GMM作为生成模型在这个41维投影空间内绘制100个新点的样本：

最后，我们可以使用PCA对象的逆变换来构造新的数字：

大部分结果看起来像数据集中合理的数字！

考虑一下我们在这里做了什么：给定一个手写数字的样本，我们已经模拟了数据的分布，这样我们就可以从数据中生成全新的数字样本：这些是“手写数字”，不是单独的出现在原始数据集中，而是捕获混合模型建模的输入数据的一般特征。这种数字生成模型可以证明作为贝叶斯生成分类器的一个组成部分非常有用，我们将在下一节中看到。

GMM估计是用于解决内生性问题的一种方法，除此之外还有TSLS两阶段最小二乘回归。

 

如果存在异方差则GMM的效率会优于TSLS，但通常情况下二者结论表现一致，很多时候研究者会认为数据或多或少存在异方差问题，因而可直接使用GMM估计。

 

内生变量是指与误差项相关的解释变量。对应还有一个术语叫‘外生变量’，其指与误差项不相关的解释变量。产生内生性的原因通常在三类，分别说明如下：

 

内生性问题的判断上 ，通常是使用Durbin-Wu-Hausman检验（SPSSAU在两阶段最小二乘回归结果中默认输出），当然很多时候会结合自身理论知识和直观专业性判断是否存在内生性问题。如果假定存在内生性问题时，直接使用两阶段最小二乘回归或者GMM估计即可。一般不建议完全依照检验进行判断是否存在内生性，结合检验和专业理论知识综合判断较为可取。

 

内生性问题的解决上 ，通常使用工具变量法，其基本思想在于选取这样一类变量（工具变量），它们的特征为：工具变量与内生变量有着相关(如果相关性很低则称为弱工具变量)，但是工具变量与被解释变量基本没有相关关系。寻找适合的工具变量是一件困难的事情，解决内生性问题时，大量的工作用于寻找适合的工具变量。

 

关于 引入工具变量的个数上 ，有如下说明：

过度识别和恰好识别是可以接受的，但不可识别这种情况无法进行建模，似想用一个工具变量去标识两个内生变量，这是不可以的。

 

工具变量引入时，有时还需要对工具变量外生性进行检验（过度识别检验），针对工具变量外生性检验上，SPSSAU提供Hansen J检验。特别提示，只有过度识别时才会输出此两个检验指标。

 

GMM估计类型参数说明如下：

 

本案例引入Mincer(1958)关于工资与受教育年限研究的数据。案例数据中包括以下信息，如下表格：

                                             

数据共有12项，其中编号为1，5，7，8，12共五项并不在考虑范畴。

 

本案例研究‘受教育年限’对于‘Ln工资’的影响。明显的，从理论上可能出现‘双向因果关系’即‘受教育年限’很可能是内生变量。那么可考虑使用‘母亲受教育年限’和‘成绩’这两项数据作为工具变量。同时研究时纳入3个外生变量，分别是‘婚姻’，‘是否大城市‘和’当前单位工作年限’。使用两阶段最小二乘TSLS回归进行解决内生性问题。

 

本案例研究时，工具变量为2个，内生变量为1个，因而为过度识别，可以正常进行TSLS回归。

关于内生性的检验Durbin-Wu-Hausman检验(SPSSAU在两阶段最小二乘回归TSLS结果中默认输出)，其用于检验是否真的为内生变量；如果说检验不通过（接受原假设），那么说明没有内生变量存在，可直接使用OLS回归即可。当然即使没有内生性，一般也可以使用TSLS回归或GMM估计，没有内生性问题时，OLS回归和TSLS回归，或GMM估计结论通常一致；

关于过度识别检验上，SPSSAU提供Hansen J检验，原理上此过度识别检验仅在‘过度识别’时才会输出，即工具变量个数>内生变量个数时，才会输出。

本案例分别将被解释变量，内生变量，工具变量和外生变量纳入对应的模型框中，如下：

 

SPSSAU共输出5类表格，分别是研究变量类型表格，GMM估计模型分析结果表格，GMM估计模型分析结果-简化格式表格，模型汇总(中间过程)表格和过度识别检验(overidentifying restrictions)。说明如下：

 

上一表格展示本次研究时涉及的各变量属性，包括被解释变量，内生变量，工具变量和外生变量组成情况。

 

上表格列出GMM估计的最终结果，首先模型通过Wald 卡方检验(Wald χ² =272.418，p=0.000<0.05)，意味着模型有效。同时R方值为0.341，意味着内生和外生变量对于工资的解释力度为34.1%。具体查看内生和外生变量对于被解释变量‘工资’的影响情况来看：

 

受教育年限的回归系数值为0.112(p=0.000<0.01)，意味着受教育年限会对工资产生显著的正向影响关系。

婚姻(已婚为1)的回归系数值为0.168(p=0.000<0.01)，意味着相对未婚群体来讲，已婚群体的工资水平明显会更高。

是否大城市(1为大城市)的回归系数值为0.145(p=0.000<0.01)，意味着相对来讲，大城市样本群体，他们的工资水平明显会更高。

当前单位工作年限的回归系数值为0.036(p=0.000<0.01)，意味着当前单位工作年限会对工资产生显著的正向影响关系。

 

总结分析可知： 受教育年限, 婚姻，是否大城市, 当前单位工作年限全部均会对工资产生显著的正向影响关系。

上表格展示模型的基础指标值，包括模型有效检验wald卡方值（此处提供wald卡方非F检验），R值，Root MSE等指标值。

过度识别检验用于检验工具变量是否为外生变量，本次研究涉及工具变量为2个，分别是‘母亲受教育年限’和‘成绩’。从上表可知，过度识别Hansen J检验显示接受原假设（p=0.868>0.05），说明无法拒绝‘工具变量外生性’这一假定，模型良好。

 

特别提示：

对GMM估计效果的检查至少包括：1，diff-in-hansen检验；2，AR test的残差序列自相关。但是从估计结果中，我没弄清楚应该如何去做这两方面。我的疑问在下文有特殊颜色的字体中。1，diff-in-hansen检验估计结果如下：Difference-in-Hansen tests of exogeneity of instrument subsets: GMM instruments for levelsHansen test excluding group: chi2(24) = 30.44 Prob >chi2 = 0.171Difference (null H = exogenous): chi2(3)= -0.00 Prob >chi2 = 1.000 iv(x3 x4)Hansen test excluding group: chi2(25) = 30.30 Prob >chi2 = 0.213Difference (null H = exogenous): chi2(2)= 0.13 Prob >chi2 = 0.935复制代码结果分2个部分：一是检验gmm()中工具变量，我们需要通过“Hansen test excluding group”的p值或是“Difference (null H = exogenous)”的p值来判断？另外，由于null H是假设工具变量外生，但gmm()中的变量是内生/弱外生变量，是否应该是p值小于显著性水平时候，才说明该部分工具变量通过检验？还是根据其他文献中，这里的p值必须大于显著性水平才通过检验？二是检验iv()中的严格外生的工具变量，按我理解，应该是要求p大于显著性水平，使得不能拒绝变量外生的原假设。

---