高斯分布（二元高斯分布）原创2021-05-31 14:47·DataStrategy

The clustering model most closely related to statistics is based on distribution models. Clusters can then easily be defined as objects belonging most likely to the same distribution. A convenient property of this approach is that this closely resembles the way artificial data sets are generated: by sampling random objects from a distribution. --wiki

在聚类模型中与统计学最相关的是基于分布的模型。然后可以将集群定义为最有可能属于同一分布的对象。这样做有一个好处就是，这与人工数据集的生成方式非常地相似，通过从分布中采样随机对象来进行。

上一期谈到的k-means聚类模型相对来说是比较简单和容易理解的。也正是因为这种简单，它在应用方面会面临一些挑战。尤其是，k均值的非概率性质及其使用简单的到集群中心距离来分配该集群成员的方法导致了k-means在现实应用场景中表现得很差。

例如，我们有一些简单的数据点，k-means算法可以快速地给那些聚类以某种方式分类，并且分类结果和我们肉眼判别的结果很接近：

%matplotlib inlineimport matplotlib.pyplot as pltimport seaborn as sns; sns.set()import numpy as np

#generate the same datafrom sklearn.datasets.samples_generator import make_blobsX, y_true = make_blobs(n_samples=400, centers=4,cluster_std=0.60, random_state=0)X = X[:, ::-1] # flip axes for better plotting

# Plot the data with k-means labelsfrom sklearn.cluster import KMeanskmeans = KMeans(4, random_state=0)labels = kmeans.fit(X).predict(X)plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')

从直观上来看，我们可能期望有些点在分配给某一个聚类的时候，比其他的点要更确定一些，例如，在两个聚类中间的那些点在分类上会显得有些重叠。因此我们可能不能完全自信地确定分类那些点到附近的哪个类别中。而且，k-means模型没有内在的概率度量，或者，集群分配的不确定性。为此，我们必须考虑对模型进行泛化。

对于k-means模型的另外一种看待方式是，它放置一个圆（或者更高维度，一个超球体）在每个聚类的中心，其半径由集群中最远的点到这个聚类的中心。

在训练集中，该半径充当训练集中的聚类分配的硬截止（可以联想一下支持向量机SVM中的硬边界和软边界）：该圆以外的任何点均不被视为集群的成员。我们可以使用以下函数可视化此集群模型。

from sklearn.cluster import KMeansfrom scipy.spatial.distance import cdist

def plot_kmeans(kmeans, X, n_clusters=4, rseed=0, ax=None):    labels = kmeans.fit_predict(X)    # plot the input data    ax = ax or plt.gca()    ax.axis('equal')    ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)    # plot the representation of the k-means model    centers = kmeans.cluster_centers_    radii = [cdist(X[labels == i], [center]).max()for i, center in enumerate(centers)]    for c, r in zip(centers, radii):        ax.add_patch(plt.Circle(c, r, fc='#CCCCCC', lw=3, alpha=0.5, zorder=1))

kmeans = KMeans(n_clusters=4, random_state=0)plot_kmeans(kmeans, X)

对上面这部分代码简单解析一下：

scipy.spatial.distance.cdist (XA, XB, metric = 'euclidean', *args, **kwargs)

计算每两个输入集之间的距离。

参数：

XA：需要计算距离的矩阵之一，XA矩阵。

XB: 需要计算距离的矩阵之一，XB矩阵。

Metric: 使用的距离测量方式。

这个测量方式function可以是 ‘braycurtis’, ‘canberara’, ‘euclidean’等很多很多。

返回:

返回的是y 矩阵，为 XA和XB 矩阵的距离矩阵, metric dis(u=XA[i], v=XB[j]) 被计算存储在距离矩阵中的ij

Radii= [cdist(X[labels == i], [center]).max() for i, center in enumerate(centers)]

计算每个聚类中，所有的点到聚类中心的距离的最大值。

Matplotlib.patches.Circle(xy, radius=5, **kwargs)

用来画圆， xy，radius分别为圆心和半径两个参数。

对于k-means一个重要的特点是，那些聚类模型是圆形的，k-means没有内置的方法来解析长方形或椭圆形聚类。因此，例如，如果我们使用同样的数据并且转换它，集群分配变得混乱。

rng = np.random.RandomState(13)X_stretched = np.dot(X, rng.randn(2, 2))kmeans = KMeans(n_clusters=4, random_state=0)plot_kmeans(kmeans, X_stretched)

通过肉眼观察，我们发现那些变形的聚类并不是圆形，因此圆形的cluster拟合度比较差。尽管如此，k-means处理这种情况不够灵活并且试图强制将数据拟合成四个圆形聚类。

这会导致被分配的集群混合再一起，其中生成的圆圈重叠：在上图的右下角尤其明显。

K-means 存在集群形状缺乏灵活性，以及缺乏概率集群分配概念 这两个缺点，意味着对于许多数据集（尤其是低维度数据集）其性能可能低于期望。

你可能会想通过generalizing （泛化）k-means模型来解决那些缺点：例如，通过比较每个点到所有聚类中心的距离来测量集群分配中的不确定性，而不是只聚焦在最近的聚类中心上。

你也可以想象使用椭圆形的聚类边界而不是圆形，用以解决非圆形的边界。事实证明这是高斯混合模型两个基本要素。

高斯混合模型试图找到一个最佳模拟任何输入数据集的多维高斯概率分布的混合。简单地举个例，GMM可以使用和k-means相同的方式来查找聚类。

from sklearn.mixture import GaussianMixture as GMMgmm = GMM(n_components=4).fit(X)labels = gmm.predict(X)plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')

GMM在内部原理上是包含一个概率模型，它也可能找到概率聚类分类--在scikit-learn中，我们使用predict_proba 方法。它会返回一个尺寸为[n_samples,n_clusters]的矩阵，用来测试所给聚类中的任何点的可能性。

我们可以可视化这种不确定性，比如，使每个点的尺寸与预测的确定性成正比。如下图，可以精确地看到在聚类的边界的地方的点的大小反应了这种属于哪个类别的不确定性。

size = 50 * probs.max(1) ** 2 # square emphasizes differencesplt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size)

高斯混合模型和k-means非常地像，它使用了一个expectation-maximization（期望最大化）方法，定性地执行了下面的操作：

1.选择开始猜测的位置和形状。

2.重复直到收敛：

A. E-step: 对于每个点，找到权重，并且编码每个成员（sample）在每个聚类中的可能性。

B. M-step: 对于每个聚类，更新它的位置，归一化和基于所有数据点上的形状，使用所有的权重。

其结果是，每个聚类都不和硬边界相关，而是一个平滑的高斯模型。 就像k-means期望最大化方法一样。这个算法有时候也会错过全局最优解，因此，实际上，使用了多个随机初始化。

下面这个function通过基于GMM输出画椭圆来帮助我们可视化GMM聚类的位置和形状。

from matplotlib.patches import Ellipsedef draw_ellipse(position, covariance, ax=None, **kwargs):# """Draw an ellipse with a given position and covariance"""    ax = ax or plt.gca()    # Convert covariance to principal axes    if covariance.shape == (2, 2):        U, s, Vt = np.linalg.svd(covariance)        angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))        width, height = 2 * np.sqrt(s)    else:        angle = 0        width, height = 2 * np.sqrt(covariance)        # Draw the ellipse    for nsig in range(1, 4):        ax.add_patch(Ellipse(position, nsig * width, nsig * height,angle, **kwargs))

def plot_gmm(gmm, X, label=True, ax=None):    ax = ax or plt.gca()    labels = gmm.fit(X).predict(X)    if label:        ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)    else:        ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)    ax.axis('equal')    w_factor = 0.2 / gmm.weights_.max()    for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):        draw_ellipse(pos, covar, alpha=w * w_factor)

gmm = GMM(n_components=4, random_state=42)plot_gmm(gmm, X)

同样地，我们可以使用GMM方法来拟合我们变形的数据集；在允许计算协方差的情况下，该模型甚至可以拟合非常长的，延申的聚类。

Np.linalg.svd(covariance)

Singular value decomposition 奇异值分解。

具体的原理和详细例子可以访问下面这篇：

PCA(Principal Components Analysis)主成分分析

经过这个奇异值分解，返回三个矩阵：unitary array(s); vector(s); unitary array(s) 用来画椭圆。

然后使用

matplotlib.patches.Ellipse(xy,width, height,angle=0,**kwargs)

画椭圆：

Xy: 椭圆中心的坐标。

Width: 水平轴的总长度。

Height: 垂直轴的总长度。

Angle: float, default: 0 ; 逆时针方向旋转的角度。

Sklearn.mixture.GMM(n_components=1, covariance_type=’diag’, random_state=None, thresh=None, tol=0.001,...)

高斯混合模型。

展示高斯混合模型概率分布，这可以轻松评估GMM分布的参数，从中采样，及对GMM分布的参数进行最大似然估计。

初始化参数，以使每个混合成分具有零均值和相同性协方差。

N_components: 混合成分数，自定义为1。

Covariance_type: 描述要使用的协方差参数类型。

必须是 ‘spherical’, ‘tied’,’diag’,’full’中的一个，默认为 ‘diag’.

个人理解这里是对协方差转换的不同方式。用于高斯混合模型中画出那个椭圆形区域，可以参见下面这个小例子，这个小例子中，无论使用的是哪种方法，里面的点的分类预测的聚类类别结果是一样的，只是在画图的时候，画出了不同形状的椭圆，比如参数选为spherical时，它画出来的区域为圆形，它们的聚类的中心点，无论选择的是哪个参数，是不会变化的，只是聚类的边界形状的一个定义。也可以由下面的连接所画的图形看到：

https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html

下面是比较了不同的’spherical’, ‘diagonal’, ‘full’ 和’tied’ 协方差矩阵，使用的是iris数据库。

Min_covar: 协方差对角矩阵的下限，以防止过拟合。默认值为：0.001

Tol: 收敛阈值，当对数似然概率的平均增益低于此阈值时，EM迭代将停止。该值默认值为1e-3.

N_iter: EM执行迭代的次数。

N_int: 执行的初始化的次数，保留最佳结果。

Params: 用来控制哪些参数在训练的过程中需要更新；这里可以是wmc的任何组合；w: weights 权重；m：means 均值；c: covars 协方差；默认值为’wmc’。

Init_params: 控制哪些参数在初始化的过程中，需要要更新；这里可以是wmc的任何组合；w: weights 权重；m：means 均值；c: covars 协方差；默认值为’wmc’。

选择协方差类型：Covariance_type这个超参数控制每个集群形状的自由度；对于任何给定的问题都需要小心谨慎地设置这个参数，covariance_type的默认值为’diag’，这意味着可以独立设置沿每个维度的簇的大小，将生成的椭圆约束为与轴对齐。

一个稍微简单和快速的模型是covariance_type=’spherical’, 它限制了簇的形状，使得所有尺寸相等（可以理解为椭圆的长轴和短轴相等）。产生的聚类结果将与k-means产生的聚类结果有相似的特征的，尽管它们不完全等效。当使用covariance_type=’full’时，模型更为复杂且计算量更大，尤其是随着维数的增加，设置’full’参数时，它允许每个聚类建模为具有任意方向的椭圆。

对于单个聚类，这三种选择的视觉上呈现如下图。

尽管GMM通常被归类为聚类算法，但从根本上讲它是用于密度估计的算法。也就是说，GMM适合用来拟合某些技术上来说不是聚类模型的数据，它通常是描述数据分布的生成概率模型。

下面使用scikit-leanr的make_moons 功能来举个例子：

from sklearn.datasets import make_moonsXmoon, ymoon = make_moons(200, noise=.05, random_state=0)plt.scatter(Xmoon[:, 0], Xmoon[:, 1])

如果我们试图使用2个component GMM 来fit这些数据，结果并不是特别有用：

gmm2 = GMM(n_components=2, covariance_type='full', random_state=0)plot_gmm(gmm2, Xmoon)

但是如果我们使用很多个components并且忽略这个聚类的label,我们找到一个fit更接近于输入的数据：

gmm2 = GMM(n_components=16, covariance_type='full', random_state=0)plot_gmm(gmm2, Xmoon)

gmm16 = GMM(n_components=16, covariance_type='full', random_state=0)plot_gmm(gmm16, Xmoon, label=False)

由上面两个图可以看出，不忽略label的情况下，就可以看出这16个cluster，也就是这16个不同的椭圆大致的分布区域。忽略label的话，可以看到图上面的小椭圆分布大致有两个连续的区域。

这里由16个高斯混合而成，但是并没有找到数据分开的聚类，而是对于所有输入的数据的整体分布进行建模。这是分布的生成模型，意味着GMM给我们提供了产生 生成新的随机数据与输入数据分布相似 的方法。比如，使用这个16-component GMM 产生400个新点，如下图。

Xnew = gmm16.sample(400)plt.scatter(Xnew[0][:,0], Xnew[0][:,1])

GMM().sample([n_samples, random_state]):从模型中随机产生样本。

GMM 是一种便利的对于任意多维数据分布的灵活建模方法。

那么对于随意输入的数据使用多少个components 才合适呢？

GMM 是一个生成模型，生成模型本质上是数据集的概率分布，因此我们可以使用交叉验证来避免过度拟合；从而简单地评估模型下数据的可能性。矫正过度拟合的另外一个方法是使用一些分析标准（例如赤池信息准则（AIC），或贝叶斯信息准则（BIC））来调整模型（如似然度）。

这两个参数在之前的文章也有提到过，详细解析可以关注gzh找到下面这篇文章查看：

时间序列分析预测未来Ⅱ SARIMA 实操手册（滚动条滑到页面中部的位置查看）

Scikit-lean的GMM算子实际上包括AIC和BIC这两种内置方法。让我们来看一下不同components下面的AIC和BIC值。

n_components = np.arange(1, 21)models = [GMM(n, covariance_type='full', random_state=0).fit(Xmoon)for n in n_components]plt.plot(n_components, [m.bic(Xmoon) for m in models], label='BIC')plt.plot(n_components, [m.aic(Xmoon) for m in models], label='AIC')plt.legend(loc='best')plt.xlabel('n_components')

最佳集群数(component的数量）是使AIC或BIC最小化的值，从AIC的曲线可以看出，选择16个components还是有点多的，大概8-12个components是更好的选择。BIC建议使用更简单的模型，即components的数量更少。

需要注意的是：components数量的选择可以测量GMM作为一个密度估算工作性能有多好，而不是作为一个聚类算法。所以先将GMM视为一个density estimator, 并仅在有必要的简单数据集时才将它用于聚类。

At last,分享一个来自University of California 的video：