gmm函数python python gmssl

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

[译] 高斯混合模型 --- python教程

本文翻译自

上一节中探讨的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对象的逆变换来构造新的数字：

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

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

相关问答

Q1: 高斯混合模型（GMM）

最近在实际工作中用到了高斯混合模型（Gaussian Mixture Model），遂写一篇笔记来整理记录相关知识点，以便复查巩固。

简单回顾一下本科概率论讲过的高斯模型。

高斯模型是一种常用的变量分布模型，又称正态分布，在数理统计领域有着广泛的应用。

当样本数据 X 是一维数据（Univariate）时，高斯分布遵从下方概率密度函数（Probability Density Function）(下文简称pdf)如下：其中为数据均值（期望），为数据标准差（Standard deviation）。

当样本数据 X 是多维数据（Multivariate）时，高斯分布pdf为：

其中，为数据均值（期望），为协方差（Covariance），描述各维变量之间的相关度，D 为数据维度。

高斯混合模型可以看作是由 K 个单高斯模型组合而成的模型，这 K 个子模型是混合模型的隐变量（Hidden variable）。一般来说，一个混合模型可以使用任何概率分布，这里使用高斯混合模型是因为高斯分布具备很好的数学性质以及良好的计算性能。

先来看一组数据。

所以，混合高斯模型并不是什么新奇的东西，它的本质就是融合几个单高斯模型，来使得模型更加复杂，从而产生更复杂的样本。理论上，如果某个混合高斯模型融合的高斯模型个数足够多，它们之间的权重设定得足够合理，这个混合模型可以拟合任意分布的样本。

对于单高斯模型，我们可以用最大似然法（Maximum likelihood）估算参数的值这里我们假设了每个数据点都是独立的（Independent），似然函数由概率密度函数（PDF）给出。

由于每个点发生的概率都很小，乘积会变得极其小，不利于计算和观察，因此通常我们用 Maximum Log-Likelihood 来计算（因为 Log 函数具备单调性，不会改变极值的位置，同时在 0-1 之间输入值很小的变化可以引起输出值相对较大的变动）：对其进行求导并令导数为0，所求出的参数就是最佳的高斯分布对应的参数。

所以最大化似然函数的意义就是：通过使得样本集的联合概率最大来对参数进行估计，从而选择最佳的分布模型。

对于高斯混合模型，Log-Likelihood 函数是：如何计算高斯混合模型的参数呢？这里我们无法像单高斯模型那样使用最大似然法来求导求得使 likelihood 最大的参数，因为对于每个观测数据点来说，事先并不知道它是属于哪个子分布的（hidden variable），因此 log 里面还有求和，对于每个子模型都有未知的，直接求导无法计算。需要通过迭代的方法求解。

EM 算法是一种迭代算法，1977 年由 Dempster 等人总结提出，用于含有隐变量（Hidden variable）的概率模型参数的最大似然估计。

每次迭代包含两个步骤：

这里不具体介绍一般性的 EM 算法，（通过 Jensen 不等式得出似然函数的下界 Lower bound，通过极大化下界做到极大化似然函数，有log(E(x))=E(log(x))），只介绍怎么在高斯混合模型里应用从来推算出模型参数。

通过 EM 迭代更新高斯混合模型参数的方法（我们有样本数据和一个有个子模型的高斯混合模型，想要推算出这个高斯混合模型的最佳参数）：

至此，我们就找到了高斯混合模型的参数。需要注意的是，EM 算法具备收敛性，但并不保证找到全局最大值，有可能找到局部最大值。解决方法是初始化几次不同的参数进行迭代，取结果最好的那次。

Q2: 单高斯模型SGM & 高斯混合模型GMM

在gmm函数python了解高斯混合模型之前gmm函数python，我们先来看看什么是高斯分布gmm函数python，高斯分布大家应该都比较熟悉了，就是我们平时所说的正态分布，也叫高斯分布。正态分布是一个在数学、物理及工程等领域都非常重要的概率分布，在统计学的许多方面有着重大的影响力。

正态分布的特点

集中性：正态曲线的高峰位于正中央，即均数所在的位置。

对称性：正态曲线以均数为中心，左右对称，曲线两端永远不与横轴相交。

均匀变动性：正态曲线由均数所在处开始，分别向左右两侧逐渐均匀下降。

若随机变量服从一个数学期望为、方差为的正态分布，记为。其中期望值决定了其位置，标准差决定了分布的幅度。当 = 0， = 1时，正态分布是标准正态分布。

正态分布有极其广泛的实际背景，生产与科学实验中很多随机变量的概率分布都可以近似地用正态分布来描述。例如，在生产条件不变的情况下，产品的强力、抗压强度、口径、长度等指标；同一种生物体的身长、体重等指标；同一种种子的重量；测量同一物体的误差；弹着点沿某一方向的偏差；某个地区的年降水量；以及理想气体分子的速度分量，等等。一般来说，如果一个量是由许多微小的独立随机因素影响的结果，那么就可以认为这个量具有正态分布（见中心极限定理）。从理论上看，正态分布具有很多良好的性质，许多概率分布可以用它来近似；还有一些常用的概率分布是由它直接导出的，例如对数正态分布、t分布、F分布等。

高斯模型有单高斯模型（SGM）和混合高斯模型（GMM）两种。

概率密度函数服从上面的正态分布的模型叫做单高斯模型，具体形式如下：

当样本数据是一维数据（Univariate）时，高斯模型的概率密度函数为：

其中：为数据的均值，为数据的标准差。

当样本数据是多维数据（Univariate）时，高斯模型的概率密度函数为：

其中：为数据的均值，为协方差，d为数据维度。

高斯混合模型（GMM）是单高斯概率密度函数的延伸，就是用多个高斯概率密度函数（正态分布曲线）精确地量化变量分布，是将变量分布分解为若干基于高斯概率密度函数（正态分布曲线）分布的统计模型。

用通俗一点的语言解释就是，个单高斯模型混合在一起，生成的模型，就是高斯混合模型。这个子模型是混合模型的隐变量（Hidden variable）。一般来说，一个混合模型可以使用任何概率分布，这里使用高斯混合模型是因为高斯分布具备很好的数学性质以及良好的计算性能。

GMM是工业界使用最多的一种聚类算法。它本身是一种概率式的聚类方法，假定所有的样本数据X由K个混合多元高斯分布组合成的混合分布生成。

高斯混合模型的概率密度函数可以表示为：

其中：

是观察数据属于第个子模型的概率，；

是第个的单高斯子模型的概率密度函数，或

，具体函数见上方单高斯模型的概率密度函数。

参数估计有多种方法，有矩估计、极大似然法、一致最小方差无偏估计、最小风险估计、同变估计、最小二乘法、贝叶斯估计、极大验后法、最小风险法和极小化极大熵法等。最基本的方法是最小二乘法和极大似然法。

极大似然估计的思想是：随机试验有多个可能的结果，但在一次试验中，有且只有一个结果会出现，如果在某次试验中，结果w出现了，则认为该结果发生的概率最大。

1）写出似然函数：

假设单个样本的概率函数为 ,对每个样本的概率函数连乘，就可以得到样本的似然函数

2）对似然函数取对数：

目的是为了让乘积变成加法，方便后续运算

3）求导数，令导数为0，得到似然方程：

和在同一点取到最大值，所以可以通过对求导，令导数为零，实现同个目的

4）解似然方程，得到的参数即为所求

对于单高斯模型，可以使用极大似然估计（MLE）来求解出参数的值。

单高斯模型的对数似然函数为：

上式分别对和求偏导数，然后令其等于0，可以得到对应的参数估计值：

如果依然按照上面的极大似然估计方法求参数

GMM的对数似然函数为：

对上式求各个参数的偏导数，然后令其等于0，并且还需要附件一个条件：。

我们会发现，直接求导无法计算出参数。所以我们需要用其它方式去解决参数估计问题，一般情况下我们使用的是迭代的方法，用期望最大算法（Expectation Maximization，EM）进行估计。

EM算法的具体原理以及示例见我的另外一篇文章。

Q3: GMM或2SLS操作问题求助

比如，在微观层面，如果面板的观测值是时序相关的，用GMM估计的动态面板就是一种最自然的解决办法；在宏观研究中，我们经常将理论模型推衍出的一阶条件作为GMM估计的矩条件（moment conditions），理论因而能够得到数据的检验。不过，GMM估计涉及到的矩条件和工具变量的选择，经常让人头疼得要命。这篇短文就是讨论GMM估计中矩条件选择的问题。我不是研究计量经济学的，很多最基本的东西都不懂，下面这些观点大都来自Victor Chernozhukov和Whitney Newey两位老师，引述的不对的地方，请大伙儿指出来。所谓矩条件，就是一个同时含有随机变量和待估计参数的式子，经济理论告诉我们，它的期望等于0。矩条件最常见的形式是：E{工具变量*残差}=0。GMM估计就是在一个限定的范围内寻找参数，使这个我们在理论上认为正确的等式填入数据后尽可能接近于0。按照我的理解，GMM不仅是一种估计方法，还是一个计量经济学有经典框架，我们能想到的大多数经典估计方法，OLS、GLS、2SLS、MD、QR、MLE、QMLE等等，都可以写成GMM的形式。另一个与之匹敌的经典框架是极值估计（extreme estimation）。粗略地说，两者的差别在于：前者是寻找参数，使矩条件尽可能被满足；后者是寻找参数，最大化或最小化一个目标函数（求极值）。简而易见的是，两种方法在算术上基本是等价的，因为任何一个极值函数的一阶条件都是矩条件，而GMM中的目标函数——矩条件经验期望的二次型——本身又是一个极值函数。但是，两者在算法上并不等价。老朱在评论林文夫（Fumio Hayashi）教授那本有名的教科书时说，“GMM的概念很优美，也可以应用到很多问题上。一般化的概念虽然适用性广，还是有代价的。”这里的代价，我猜测，就是指GMM经常算不出来——由于矩条件本身的特点，GMM的目标函数经常是接近锯齿状的（piecewise constant），在这种情况下，GMM会陷在局部最优里，达不到全局最优。分位数回归（quantile regression）就是一个这样的例子。不过，GMM在很多时候还是有用的，而且算起来特别快，所以还是有讨论的必要。首先我们要问，在GMM估计中，矩条件是不是最多越好？在大样本下，基本上是这样——矩条件越多，GMM估计的渐近效率就越高。说“基本上”是因为：第一，这些矩条件必须都是成立的；第二，矩条件的数目相对于样本数要趋向于0。如果矩条件数与样本数是等阶的，会造成“过度拟合”的问题。形象地说，我们本来要用工具变量来应付内生性问题，但是工具变量太多了，以至于几乎把内生变量完全拟合了出来，那么即使工具变量是外生的，也会导致估计量不一致。而在小样本下，过多的矩条件会造成可怕的高阶偏误，并且矩条件非线性的程度越高，偏误就越大。要注意的是，这个问题与“弱工具变量问题”并不等价。即使这些工具变量整体上不弱，甚至每个都不弱，过多的矩条件还是会造成严重的小样本偏误。这里其实涉及到GMM估计的高阶偏误问题，其核心是由GMM目标函数的“非线性”特征造成的。当矩条件的数目很多，矩条件本身又是非线性的时候，这个问题就愈加严重；但即使没有“矩条件过多”的问题，GMM仍然存在不可忽视的小样本偏误。在实践中，我们可以用几种方法来减轻这一问题。第一，检验矩条件（或工具变量）是否成立。 “过度识别检定”（overidentification test，OIT）可以被用检验某一组矩条件是否成立，前提是去除待检验的矩条件后，剩余的矩条件数目仍大于等于待估计参数的维度。第二，选择“最有效率”的工具变量。给定一个理论上有效的工具变量Zi，我们可以通过简便的方法找出Zi的某种最优的函数形式f*(Zi)，把f*(Zi)放入矩条件会使得估计量的渐近方差比放入其他f(Zi)要小。这么做可以尽可能地利用Zi中的信息，而不必将不同函数形式的Zi写成并列的几个矩条件。第三，在所有成立的矩条件中选择一组最优的矩条件。用任意组合的矩条件进行估计，看其中哪一组矩条件得到的估计量的“经验均方误”最小。最后，我们还可以用Fuller、HFUL等K-class估计量或LIML、CUE等经验似然估计量来进行估计，然后用Bekker估计量校正其标准误。这些估计量可以在很大程度上减小偏误，即使无法完全消除它。

Q4: Python hmmlearn中的混淆矩阵是怎么表示的

hmmlearn这个库有三种模型，分别是Gaussian，Multinomial和GMMHMM。这三种模型对应的就是三种emission matrix（即混淆矩阵，也就是隐状态到观察态的概率）。Gaussian就是说混淆矩阵是一个高斯分布，即观察态是连续的。Multinomiual就是说混淆矩阵事一个Multibimiual distribution，即观察态势离散的。GMMHMM则是说混淆矩阵是遵循gaussinan mixture 分布，也是连续的。

题主问如何把混淆矩阵输入到模型里面。首先你要确定你的混淆矩阵的类型。对于Gaussian类型，就是把你希望的 mean和variance值放到模型里面。我就直接把文档里面的例子搬过来，例子里是建立了一个高斯分布的隐马尔科夫模型。

import numpy as np

from hmmlearn import hmm

#一个隐马尔科夫模型由（p向量，状态转移矩阵，混淆矩阵）来定义。

startprob = np.array([0.6, 0.3, 0.1])

# 定义初始状态的概率

transmat = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2], [0.3, 0.3, 0.4]])#定义转移矩阵的概率

means = np.array([[0.0, 0.0], [3.0, -3.0], [5.0, 10.0]])

#定义混淆矩阵的均值

covars = np.tile(np.identity(2), (3, 1, 1))# 定义混淆矩阵的方差

model = hmm.GaussianHMM(3, "full", startprob, transmat)# 定义一个混淆矩阵为高斯分布的隐马尔科夫模型。这里‘full’的意思就是说你输入的方差矩阵每个元素都给出了，不是一个只是对角线上的元素为0的矩阵

model.means_ = means