深度学习方法拟合高斯核模型 高斯拟合优点

一、k-means聚类的缺点2维k-means模型的本质是，它以每个簇的中心为圆心，簇中点到簇中心点的欧氏距离最大值为半径画一个圆。这个圆硬性的将训练集进行截断。而且，k-means要求这些簇的形状必须是圆形的。因此，k-means模型拟合出来的簇(圆形)与实际数据分布(可能是椭圆)差别很大，经常出现多个圆形的簇混在一起，相互重叠。总的来说，k-means存在两个缺点，使得它对许多数据集(特别是低维数据集)的拟合效果不尽如人意：

1. 类的形状不够灵活，拟合结果与实际相差较大，精度有限。
2. 样本属于每一个簇的概率是定性的，只有是与否，不能输出概率值。应用中缺少鲁棒性。

二、高斯混合模型高斯混合模型(GMM)可以看做是k-means模型的一个优化。它既是一种工业界常用的技术手段，也是一种生成式模型。高斯混合模型试图找到多维高斯模型概率分布的混合表示，从而拟合出任意形状的数据分布。在最简单的场景中，GMM可以用与k-means相同的方式进行聚类。

import matplotlib.pyplot as pltimport seaborn as sns; sns.set()import numpy as np#产生实验数据from sklearn.datasets.samples_generator import make_blobsX, y_true = make_blobs(n_samples=700, centers=4,                       cluster_std=0.5, random_state=2019)X = X[:, ::-1] #方便画图from sklearn.mixture import GaussianMixture as GMMgmm = GMM(n_components=4).fit(X) #指定聚类中心个数为4labels = gmm.predict(X)plt.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis')
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np

#产生实验数据
from sklearn.datasets.samples_generator import make_blobs
X, y_true = make_blobs(n_samples=700, centers=4,
                       cluster_std=0.5, random_state=2019)
X = X[:, ::-1] #方便画图

from sklearn.mixture import GaussianMixture as GMM
gmm = GMM(n_components=4).fit(X) #指定聚类中心个数为4
labels = gmm.predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis')

它使用EM算法进行迭代：1.选择位置和初始形状2.循环直至收敛：    E步骤：对于每个点，为每个点分别计算由该混合模型内的每个分量生成的概率。    M步骤：调整模型参数以最大化模型生成这些参数的可能性。该算法保证该过程内的参数总会收敛到一个局部最优解。

三、GMM中的概率模型

事实上在GMM算法中，有一个隐含的概率模型。可以通过其得到簇分配结果的概率。打印前十个点分别属于四个类的概率。

probs = gmm.predict_proba(X)print(probs[:10].round(2))
probs = gmm.predict_proba(X)
print(probs[:10].round(2))

因为GMM模型并不是通过硬截断进行分割类别，而是通过高斯平滑模型进行估计的。所以将每个点的概率进行可视化时，散点图并不是严格成椭圆形状的。

size = probs.max(1)plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size)
size = probs.max(1)
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size)

如果允许使用全部的协方差类型，则可以拟合任意形状的分布。为了更好的展示GMM模型的拟合结果，首先需要构造一个画椭圆的函数。在网上找到的代码因为一些API有改动，重新更新了一版。

from matplotlib.patches import Ellipse#给定的位置和协方差画一个椭圆def draw_ellipse(position, covariance, ax=None, **kwargs):    ax = ax or plt.gca()    #将协方差转换为主轴    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)    #画出椭圆    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=4, cmap='viridis', zorder=2)    else:        ax.scatter(X[:, 0], X[:, 1], s=4, 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)        
from matplotlib.patches import Ellipse

#给定的位置和协方差画一个椭圆
def draw_ellipse(position, covariance, ax=None, **kwargs):
    ax = ax or plt.gca()
    #将协方差转换为主轴
    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)

    #画出椭圆
    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=4, cmap='viridis', zorder=2)
    else:
        ax.scatter(X[:, 0], X[:, 1], s=4, 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)

下面使用椭圆形来拟合数据。

rng = np.random.RandomState(13)X_stretched = np.dot(X, rng.randn(2, 2))gmm = GMM(n_components=4, covariance_type='full', random_state=42)plot_gmm(gmm, X_stretched)
rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))
gmm = GMM(n_components=4, covariance_type='full', random_state=42)
plot_gmm(gmm, X_stretched)

四、GMM模型的组件下面考虑一个特殊的分布形式。如下图所示

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

如果使用两个高斯分布进行拟合，则得到的结果如下。

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

即使是椭圆形状，仍有一部分点被错误的归类为另一个分布。这时，如果使用更多的高斯分布进行归纳，则可以得到更好的效果。

gmm10 = GMM(n_components=10, covariance_type='full', random_state=0)plot_gmm(gmm10, Xmoon, label=False

这里使用了10个高斯分布。 但是并不是为了得到10个聚类结果。 而是通过10个分布进行集成得到最终的归纳结果。 也就是说，GMM模型的本质并不是聚类，而是得到一个，能够生成当前样本形式的分布。因此可以使用前面10个高斯分布集成的生成模型，来生成服从当前分布形式的200个新的样本。

Xnew = gmm10.sample(200)[0]plt.scatter(Xnew[:, 0], Xnew[:, 1])

五、最优组件个数的确定为了确定最优组件的个数，需要

使用一些分析标准来调整模型可能性。模型中封装了Akaike information criterion (AIC)  或 Bayesian information criterion (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')
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')

最佳的聚类数目是使得AIC或BIC最小化的值，具体取决于我们希望使用的近似值。AIC告诉我们，我们上面选择的10个组件差不多就是最优解了。而BIC则倾向于使用更简单的模型。 六、总结

GMM模型因其优秀的聚类表现，以及可以生产样本的强大功能，在风控领域的应用非常广泛。如对反欺诈中的欺诈样本抓取与生成、模型迭代中的幸存者偏差等问题都有一定的作用。 比如说在反欺诈检测模型中，可以先通过GMM模型对欺诈样本进行聚类，后将聚类后得到的簇作为同一类欺诈手段，后续只针对不同的簇进行建模，在实践中对欺诈客户的召回有很好的效果。