8-无监督学习技术

理论部分

无监督学习算法的用途

• 聚类
  • 目标是将相似的实例分组到集群中。聚类是很好的工具，用于数据分析、客户细分、推荐系统、搜索引擎、图像分割、半监督学习、降维等
• 异常检测
  • 目的是学习“正常”数据看起来是什么样的，然后将其用于异常检测情况，例如生产线上的缺陷产品或时间序列中的新趋势
• 密度估算
  • 这是估计生成数据集的随机过程的概率密度函数(PDF)的任务，密度估算通常用于异常检测：位于非常低密度区域的实例很可能是异常，它对于数据分析和可视化也很有用

聚类用途

• 客户细分
  • 你可以根据客户的购买记录和他们在网站上的活动对客户进行聚类，这对于了解你的客户是谁以及他们的需求很有用，因此你可以针对每个细分客户条调整产品和营销活动
• 数据分析
  • 在分析数据集时，运行聚类算法然后分别分析每个集群
• 降维技术
  • 数据集聚类后，通常可以测量每个实例与每个集群的相似度(相似度是衡量一个实例和一个集群的相似程度)。然后可以将每个实例的特征向量x替换为其集群相似度向量。如果有k个集群，则词向量为k维。词向量的维度通常比原始特征向量低得多，但它可以保留足够的信息以进行进一步处理
• 异常检测(也称为离群值检测)
  • 对所有集群具有低相似度的任何实例都可能是异常。异常检测在检测制造生产线中的缺陷或欺诈检测中特别有用
• 半监督学习
  • 如果你只有几个标签，则可以执行聚类并架构标签传播到同一集群中的所有实例。该技术可以大大增加可用于后续有监督学习算法的标签数量，从而提高其性能
• 搜索引擎
  • 一些搜索引擎可让你搜索与参考图像相似的图像。要构建这样的系统，首先要对数据库中的所有图像应用聚类算法，相似的图像最终会出现在同一集群中。然后，当用户提供参考图像时，你需要做的就是使用训练好的聚类模型找到该图像的集群，然后可以简单地从该集群中返回所有的图像
• 分割图像
  • 通过根据像素的颜色对像素进行聚类，然后用其聚类的平均颜色替换每个像素的颜色，可以显著减少图像中不同颜色的数量。图像分割用于许多物体检测和跟踪系统中，因为它可以更轻松地检测每个物体的轮廓

K-Means

• 你必须要指定这个算法必须要找到的集群数k
• 每个实例都会被分配给k个集群之一
• 当集群具有不同的直径时，K-Means算法的性能不是很好，因为将实例分配给某个集群时，它所关心的只是与中心点的距离
• 与其将每个实例分配给一个单独的集群(称为硬聚类)，不如为每个实例赋予每个集群一个分数(称为软聚类)会更有用。分数可以是实例与中心点之间的距离。相反，它可以是相似性分数(或远近程度)，例如高斯径向基函数
• 如果你有一个高维数据集并以这种方式进行转换，那么最终将得到一个k维数据集：这种转换是一种非常有效的非线性降维技术
• 原理
  • 只要从随机放置中心点开始(例如，随机选择k个实例并将其位置用作中心点)。然后标记实例，更新中心点，标记实例，更新中心点，以此类推，直到中心点停止移动。这个算法能保证在有限数量的步骤内收敛(通常很小)，不会永远震荡
  • 尽管算法保证会收敛，但它可能不会收敛到正确解(即可能会收敛到局部最优值)：这个取决于中心点的初始化
• 中心点初始化方法
  • 如果你碰巧知道了中心点应该在哪里，则可以将超参数init设置为包含中心点列表的Numpy数组，并将n_init设置为1
  • 另一种解决方案是使用不同的随机初始化多次运行算法，并保留最优解。如何确切知道哪种解答是最优解呢？它使用性能指标，这个指标称为模型的惯性，即每个实例与其最接近的中心点之间的均方距离，选择惯性最小的解
  • K-Means++
    • 这是一种更智能的初始化步骤，该步骤倾向于选择彼此相距较远的中心点，这一改进使得K-Means算法收敛到次优解的可能性很小

[1、取一个中心点c^{(1)}，从数据集中随机选择一个中心点\ 2、对于每个实例x^{(i)}，计算D(x^{(i)})=min||x^{(i)}-μ||^2_2,i=1,2,...,k_{selected}\ 选择一个新的数据点作为新的聚类中心，选择的原则是D(x)较大的点被选择作为\ 中心点的概率较大\ 3、重复上一步，直到选择了所有k个中心点 ]

• 加速的K-Means和小批量K-Means
  • Elkan通过利用三角不等式(即一条直线是两点之间的最短距离)并通过跟踪实例和中心点之间的上下限距离来实现这一目的
  • 小批量K-Means在每次迭代中使用小批量K-Means稍微移动中心点，而不是在每次迭代中使用完整的数据集。这将算法的速度提高了3到4倍，并且可以对不容纳内存的大数据进行聚类。尽管小批量K-Means算法比常规K-Means算法快得多，但其惯性通常略差一些，尤其是随着集群的增加
• 寻找最佳聚类数
  • 惯性
    • 尝试选择k时，惯性不是一个好的性能指标，因为随着k的增加，惯性会不断降低。实际上，集群越多，每个实例将越接近于其最近的中心点，因此惯性将越低。当绘制惯性作为集群k的函数曲线时，曲线通常包含一个称为“肘”的拐点，如果我们不知道更好的选择，则该拐点是一个不错的选择：较低的值变化比较大，而较大的值没有多大帮助，有时可以把集群再细分成两半
  • 轮廓分数
    • 轮廓分数是所有实例的平均轮廓系数。实例的轮廓系数等于(b-a)/max(a,b)，其中a是与同一集群中其它实例的平均距离(即集群内平均距离)，b是平均最近集群距离(即到下一个最近集群实例的平均距离，定义为使b最小的那个，不包括私立自身的集群)。轮廓系数在-1和+1之间变化。接近+1的系数表示该实例很好地位于其自身的集群中，并且原理其他集群，而接近0的系数表示该实例接近一个集群的边界，接近-1的系数意味着该实例可能已分配给错误的集群
  • 轮廓图
    • 当你绘制每个实例的轮廓系数时，可以获得更加丰富的可视化效果，这些轮廓系数按它们分配的集群和系数值进行排序，这称为轮廓图。每个轮廓图包含一个刀形的聚类。形状的高度表示集群包含的实例数，宽度表示集群中实例的排序轮廓系数(越宽越好)。虚线表示平均轮廓系数。垂直虚线表示每个集群的轮廓分数。当集群中大多数的实例的系数均低于此分数时，则该集群比较糟糕，因为这意味着其实例太接近其他集群了
• K-Means的局限
  • 必须多次运行该算法才能避免次优解
  • 需要指定集群数，这很麻烦
  • 当集群具有不同的大小、不同的密度或非球形时，K-Means的表现也不佳
  • 在运行K-Means之前，请先缩放输入特征，否则集群可能会变形，这样K-Means的性能会很差。缩放特征并不能保证所有的集群都很好，但是通常可以改善结果

DBSCAN

• 该算法说明了一种基于局部密度估计的不同的方法，这种方法允许算法识别任意形状的聚类
• 该算法将集群定义为高密度的连续区域，原理如下：
  • 对于每个实例，该算法都会计算在距它一小段距离ε内有多少个实例，该区域称为实例的ε-邻域
  • 如果一个实例在其ε邻域中至少包含min_samples个实例(包括自身)，则该实例被视为核心实例，换句话说，核心实例是位于密集区域中的实例
  • 核心实例附近的所有实例都属于同一集群。这个邻域可能包括其他核心实例，因此，一长串相邻的黑心实例形成一个集群
  • 任何不是核心实例且邻居中没有实例的实例都被视为异常
• 尽管DBSCAN类具有fit_predict()方法，但是它没有predict()方法。换句话说，它无法预测新实例属于哪个集群。之所以做出这个决定，是因为不同的分类算法可以更好地完成不同的任务，因此作者决定让用选择使用哪种算法
  • 我们可以仅在核心实例不上训练分类器，但是我们也可以选择在所有实例或除异常之外的所有实例上训练分类器，这取决于你的任务。
  • 直接引入最大距离的概念，在这种情况下，距离两个集群均较远的两个实例被归类为异常
• DBSCAN是一种非常简单但功能强大的算法，能够识别任何数量的任何形状的集群。它对异常值具有鲁棒性，并且只有两个超参数(eps和min_samples)。但是，如果密度在整个集群中变化很大，则可能无法正确识别所有的集群

其他聚类方法

• 聚集聚类
  • 在每次迭代中，聚集聚类连接最近的一对集群(从单个实例开始)。如果为每对合并的集群绘制一个带有分支的树，将得到一个二叉树集群，其中叶子是各个实例。这种方法可以很好地扩展到大量实例或集群，它可以识别各种形状的集群，生成灵活且丰富的集群树，而不强迫你选择特定的集群规模，并且可以与任何成对的距离一起使用
• BIRCH
  • BIRCH(使用层次结构的平衡迭代约简和聚类)算法是专门为非常大的数据集设计的，并且只要特征数量不太大(<20)，它可以逼批量处理的K-Means更快，结果相似。在训练期间，它构建了一个树结构，该树结构仅包含足以将每个新实例快速分配给集群的信息，而不必将所欲实例存储在树种，这种方法使它在处理大型数据集时使用有限的内存
• 均值漂移
  • 此算法从开始在每个实例上居中放置一个圆圈。然后对于每个圆，它会计算位于其中的所有实例的均值，然后移动圆，使其以均值为中心。接下来，迭代次均值移动步骤，直到所有的圆都停止移动(直到每个圆都以其包含的实例的均值为中心)为止。均值漂移将圆向较高密度的方向移动，直到每个圆都找到了局部最大密度。
  • 均值漂移与DBSCAN有相同的特性，例如找到具有任何形状的任意数量的集群
  • 不适用于大型数据集
• 相似性传播
  • 此算法使用投票机制，实例对相似的实例进行投票并将其作为代表，一旦该算法收敛，每个代表及其投票者将组成一个集群
  • 不适用于大型数据集
• 谱聚类
  • 该算法采用实例之间的相似度矩阵，并由其创建低维嵌入(即降维)，然后在该低维空间中使用另一中聚类算法。谱聚类可以识别复杂的集群结构，也可以用于分割图(例如，识别社交网络上的朋友集群)。它不能很好地扩展到大量实例，并且当集群的大小不同时，它的效果也不好

高斯混合模型

• 高斯混合模型(GMM)是一种概率模型，它假定实例是由多个未知的高斯分布的混合生成的。从单个高斯分布生成的所有实例都形成一个集群，通常看起来像一个椭圆。每个集群都可以具有不同的椭圆形状、大小、密度和方向
• 假定数据集X是通过以下概率过程生成的

[对于每个实例，从k个集群中随机选择一个集群。选择第j个集群的概率是由集群的权重phi^{(j)}定义。\ 第i个实例选择的集群的索引记为z^{(i)} \ 如果z^{(i)}=j，表示第i个实例已分配给第j个集群，则从高斯分布中随机采样该实例的位置x^{(i)}，\ 其中均值mu^{(j)}和协方差矩阵{sum}^{(j)}。这记为x^{(i)}～N(mu^{(j)},{sum}^{(j)})\ 给定数据集X，首先估计出权重phi及所有分布参数mu^{(1)}至mu^{(k)}和{sum}^{(1)}至{sum}^{(k)} ]

• GaussianMixture类依赖于期望最大化算法(EM)，该算法与K-Means算法有很多相似之处：它随机初始化集群参数，然后重复两个步骤直到收敛，首先将实例分配给集群(这称为期望步骤)，然后更新集群(这称为最大化步骤)。
  • EM使用软集群分配而不是应分配，对于每个实例，在期望步骤中，该算法(基于当前的集群参数)估计它属于每个集群的概率。然后，在最大化步骤中，使用数据集中所有实例来更新每个聚类，并通过每个实例属于该集群的估计概率对其进行加权。这些概率称为实例集群的职责。在最大化步骤中，每个集群的更新主要收到其最负责的实例的影响
  • EM最终可能会收敛到较差的解，因此它需要运行几次，仅仅保留最优解
  • 当存在多个维度、多个集群或很少实例时，EM可能难以收敛到最佳解决方案。需要通过限制算法必须学习的参数数量来降低任务的难度。一种方法是限制集群可能具有的形状和方向的范围，这可以通过对协方差矩阵施加约束来实现
    • spherical
      • 所有集群都必须是球形的，但它们可以具有不同的直径(即不同的方差)
    • diag
      • 集群可以采用任何大小的任意椭圆形，但是椭圆形的轴必须平行于坐标轴(即协方差矩阵必须是对角线)
    • tied
      • 所有集群必须具有相同的椭圆形状、大小和方向(即所有集群共享相同的协方差矩阵)

使用高斯混合进行异常检测

• 异常检测(也称为离群值检测)是检测严重偏离标准的实例的任务。这些实例称为异常或离群值，而正常实例称为内值。使用高斯混合模型进行异常检测非常简单：位于低密度区域的任何实例都可以被视为异常。你必须定义需要使用的密度阈值。
• 高斯混合模型会尝试拟合所有数据，包括异常数据。因此，如果你有太多异常值，这将使模型的“正常性”产生偏差，某些离群值可能被错误地视为正常数据。如果发生这种情况，你可以尝试拟合模型一次，使用它来检测并去除最极端的离群值，然后再次将模型拟合到清理后的数据集上
• 选择聚类数
  • 对于高斯混合模型，无法使用惯性或轮廓分数作为度量标准，因为当集群不是球形或具有不同大小时，它们不可靠。相反，你可以尝试找到最小化理论信息量准则的模型
  • BIC和AIC都对具有更多要学习的参数(例如，更多的集群)的模型进行惩罚，并奖励非常拟合数据的模型。它们常常最终会选择相同的墨西哥，当它们不同时，BIC选择的模型往往比AIC选择的模型更简单(参数更少)，但往往不太拟合数据(对于较大的数据集尤其如此)

[贝叶斯信息准则(BIC)和赤池信息准则(AIC)：\ BIC=log(m)p-2log(hat L)\ AIC=2p-2log(hat L)\ 在这些等式中:\ m是总实例数，p是模型学习的参数数量，hat L是模型的似然函数的最大值 ]

似然函数

• 给定具有一些参数θ的统计模型，用“概率”一词描述未来的结果x的合理性(知道参数θ)，而用“似然”一词表示描述在知道结果x之后，一组特定参数值θ的合理性
• PDF是x的函数(固定了θ)，而似然函数是θ的函数(固定了x)
• 给定数据集X，一项常见的任务是尝试估计模型参数最有可能的值。为此，你必须找到给定X的最大化似然函数的值。如果存在关于θ的先验概率分布g，则可以通过最大化ζ(θ|x)g(θ)来考虑，而不仅仅是最大化ζ(θ|x)。这称为最大后验(MAP)估计。由于MAP会约束参数值，因此你可以将其视为MLE的正则化版本

贝叶斯高斯混合模型

• 该模型可以不用手动搜索集群的最佳数目，它为不必要的集群赋予等于(或接近于)零的权重。将集群数据n_components设置为一个你有充分理由相信的值，该值大于最佳集群的数量，算法会自动消除不必要的集群
• 在此模型中，集群参数(包括权重、均值和协方差矩阵)不再被视为固定的模型参数，而是被视为潜在的随机变量，像集群分配，因此，z现在同时包括集群参数和集群分配
• Beta分布通常用于对值在固定范围内的随机变量进行建模。在这种情况下，范围是0到1
  • 突破性过程(SBP)：假设Φ=[0.3, 0.6, 0.5, ...]，那么分配30%的实例到集群0，剩余实例的60%分配给进群1，然后将剩余实例的50%分配给集群2，以此类推
  • 在这种数据集中，新实例相比小集群更可能加入大集群(例如人们更可能迁移到较大的城市)。如果α集中度高，则Φ值很可能接近0，并且SBP会生成许多集群。相反，如果集中度低，则Φ值很可能接近1，并且集群很少。最后，使用Wishart分布对协方差矩阵进行采样：参数d和V控制集群形状的分布

• 贝叶斯定理告诉我们，在观察到一些数据X之后，如何更新潜在变量的概率分布。它计算后验分布p(z|X)，这是给定X的z的条件概率。

[贝叶斯定理：\ p(z|X) = 后验=frac{概率*先验}{证据}=frac{p(X|z)p(z)}{p(X)}\ 概率分布p(z)中有关潜在变量z的先验知识称为先验\ 不幸的是，在高斯混合模型(以及许多其他问题)中，分母p(X)很难处理，\因为它需要对z的所有可能值进行积分，\ 这需要考虑所有可能的集群参数和参数分配的组合：\ p(X) = ∫p(X|z)p(z)dz\ 解决方法之一是变分推理，它选择具有变分参数λ的分不族q(z;lambda)，\然后优化这些参数以使q(z)成为\ p(z|X)的良好近似值。这是通过找到最小化从q(z)到p(z|X)的KL散度的lambda值来实现的，\记为D_{KL}(q||p)\ 可以将KL散度方程重写为证据的对数(log p(X))减去证据的下届(ELBO)。\由于证据的对数不依赖于q，\ 它是一个常数项，因此，使KL散度最小化只需使ELBO最大化\ 判断p、q两个分布逼近的程度用KL散度表示，\ KL(q(z)||p(z|X))=∫_zq(z)logfrac{q(z)}{p(z|X)}dz,\ KL散度越小，逼近程度越大。但KL散度难以直接求，根据贝叶斯公式：\ frac{p(z,X)}{q(z)}frac{q(z)}{p(z|X)}=p(X)\ 取对数并求关于z的期望：\ ∫_zq(z)logfrac{p(z, X)}{q(z)}dz+∫q(z)logfrac{q(z)}{p(z|X)}dz=∫_zq(z)logp(x)dz\ =ELBO+KL(q(z)||p(z|X))\ 其中，p(x)是客观固定的，∫_zq(z)logp(x)dz=logp(x)也是固定的\ 所以这里KL散度的大小完全由ELBO决定，只需最大化ELBO，\既可达到最小化KL散度的目的 ]

• 一种最大化ELBO的更简单方法称为黑盒随机变分推理(BBSVI)：在每次迭代中，从q中抽取一些样本，然后使用它们来估计ELBO关于变分参数λ的梯度，然后将其用于梯度上升步骤。这种方法使得可以将贝叶斯推理用于任何类型的模型(假设它是可微分的)，甚至是深度神经网络。将贝叶斯推理与深度神经网络结合使用称为贝叶斯深度学习

其他用于异常检测和新颖性检测的算法

• PCA
  • 如果你将正常实例的重建误差与异常的重建误差进行比较，则后者通常会大得多
• Fast-MCD(最小方差协定)
  • 假定正常实例是根据一个单个高斯分布生成的，它还假定数据集被不是由该高斯分布生成的异常值所污染
• 隔离森林
  • 该算法构建一个随机森林，其中每个决策树都是随机生长的：在每个节点上，它随机选择一个特征，然后选择一个随机阈值来架构数据集分成两部分。数据集以这种方式逐渐切成分支，直到所有实例最终都与其他实例隔离开来。异常通擦与其他实例相距甚远，因此平均而言，与正常实例相比，它们倾向于用更少的步骤被隔离
• 局部离群因子(LOF)
  • 它将给定实例周围的实例密度与其相邻实例周围的密度进行比较。异常值通常比k个最近的邻居更孤立
• 单类SVM
  • 由于只有单类实例，因此单类SVM算法会尝试将高维空间中的实例与原点分开。在原始空间中，这将对应于找到一个包含所有实例的小区域，如果新实例不在该区域内，则为异常

代码部分

引入

import sys
assert sys.version_info >= (3, 5)

import sklearn
assert sklearn.__version__ >= '0.20'

import numpy as np
import os

np.random.seed(42)

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)

PROJECT_ROOT_DIR = '.'
CHAPTER_ID = 'unsupervised_learning'
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, 'images', CHAPTER_ID)
os.makedirs(IMAGES_PATH, exist_ok=True)

def save_fig(fig_id, tight_layout=True, fig_extension='png', resolution=300):
    path = os.path.join(IMAGES_PATH, fig_id + '.' + fig_extension)
    print('Saving figure', fig_id)
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format=fig_extension, dpi=resolution)

分类与聚类

from sklearn.datasets import load_iris

data = load_iris()
X = data.data
y = data.target
data.target_names  # array(['setosa', 'versicolor', 'virginica'], dtype='<U10')

plt.figure(figsize=(9, 3.5))

plt.subplot(121)
plt.plot(X[y==0, 2], X[y==0, 3], 'yo', label='Iris setosa')
plt.plot(X[y==1, 2], X[y==1, 3], 'bs', label='Iris versicolor')
plt.plot(X[y==2, 2], X[y==2, 3], 'g^', label='Iris virginica')
plt.xlabel('Petal length', fontsize=14)
plt.ylabel('Petal width', fontsize=14)
plt.legend(fontsize=12)

plt.subplot(122)
plt.scatter(X[:, 2], X[:, 3], c='k', marker='.')
plt.xlabel('Petal length', fontsize=14)
plt.tick_params(labelleft=False)

save_fig('classification_vs_clustering_plot')
plt.show()

from sklearn.mixture import GaussianMixture
from scipy import stats

y_pred = GaussianMixture(n_components=3, random_state=42).fit(X).predict(X)
mapping = {}
for class_id in np.unique(y):
    # 寻找数组或者矩阵每行/每列中最常出现成员以及出现的次数
    mode, _ = stats.mode(y_pred[y==class_id])
    mapping[mode[0]] = class_id
mapping  # {1: 0, 2: 1, 0: 2}

y_pred = np.array([mapping[cluster_id] for cluster_id in y_pred])

plt.plot(X[y_pred==0, 2], X[y_pred==0, 3], 'yo', label='Cluster 1')
plt.plot(X[y_pred==1, 2], X[y_pred==1, 3], 'bs', label='Cluster 2')
plt.plot(X[y_pred==2, 2], X[y_pred==2, 3], 'g^', label='Cluster 3')
plt.xlabel('Petal length', fontsize=14)
plt.ylabel('Petal width', fontsize=14)
plt.legend(loc='upper left', fontsize=12)
plt.show()

np.sum(y_pred == y)  # 145

np.sum(y_pred == y) / len(y_pred)  # 0.9666666666666667

K-Means

from sklearn.datasets import make_blobs

blob_centers = np.array([
    [0.2, 2.3],
    [-1.5, 2.3],
    [-2.8, 1.8],
    [-2.8, 2.8],
    [-2.8, 1.3]
])
blob_std = np.array([0.4, 0.3, 0.1, 0.1, 0.1])

X, y = make_blobs(n_samples=2000, centers=blob_centers, cluster_std=blob_std, random_state=7)

def plot_clusters(X, y=None):
    plt.scatter(X[:, 0], X[:, 1], c=y, s=1)
    plt.xlabel('$x_1$', fontsize=14)
    plt.ylabel('$x_2$', fontsize=14, rotation=0)

plt.figure(figsize=(8, 4))
plot_clusters(X)
save_fig('blobs_plot')
plt.show()

from sklearn.cluster import KMeans

k = 5
kmeans = KMeans(n_clusters=k, random_state=42)
y_pred = kmeans.fit_predict(X)

y_pred  # array([4, 0, 1, ..., 2, 1, 0])

y_pred is kmeans.labels_  # True

kmeans.cluster_centers_
'''
array([[-2.80389616,  1.80117999],
       [ 0.20876306,  2.25551336],
       [-2.79290307,  2.79641063],
       [-1.46679593,  2.28585348],
       [-2.80037642,  1.30082566]])
'''

kmeans.labels_  # array([4, 0, 1, ..., 2, 1, 0])

X_new = np.array([[0, 2], [3, 2], [-3, 3], [-3, 2.5]])
kmeans.predict(X_new)  # array([1, 1, 2, 2])

决策边界

def plot_data(X):
    plt.plot(X[:, 0], X[:, 1], 'k.', markersize=2)
    
def plot_centroids(centroids, weights=None, circle_color='w', cross_color='k'):
    if weights is not None:
        centroids = centroids[weights > weights.max() / 10]
    plt.scatter(centroids[:, 0], centroids[:, 1], marker='o', s=35, linewidths=8, color=circle_color, zorder=10, alpha=0.9)
    plt.scatter(centroids[:, 0], centroids[:, 1], marker='x', s=2, linewidths=12, color=cross_color, zorder=11, alpha=1)

def plot_decision_boundaries(clusterer, X, resolution=1000, show_centroids=True, show_xlabels=True, show_ylabels=True):
    mins = X.min(axis=0) - 0.1
    maxs = X.max(axis=0) + 0.1
    xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution), np.linspace(mins[1], maxs[1], resolution))
    Z = clusterer.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    
    plt.contourf(Z, extent=(mins[0], maxs[0], mins[1], maxs[1]), cmap='Pastel2')
    plt.contour(Z, extent=(mins[0], maxs[0], mins[1], maxs[1]), linewidths=1, color='k')
    plot_data(X)
    if show_centroids:
        plot_centroids(clusterer.cluster_centers_)
    if show_xlabels:
        plt.xlabel('$x_1$', fontsize=14)
    else:
        plt.tick_params(labelbottom=False)
    if show_ylabels:
        plt.ylabel('$x_2$', fontsize=14, rotation=0)
    else:
        plt.tick_params(labelleft=False)

plt.figure(figsize=(8, 4))
plot_decision_boundaries(kmeans, X)
save_fig('voronoi_plot')
plt.show()

硬聚类与软聚类

kmeans.transform(X_new)
'''
array([[2.81093633, 0.32995317, 2.9042344 , 1.49439034, 2.88633901],
       [5.80730058, 2.80290755, 5.84739223, 4.4759332 , 5.84236351],
       [1.21475352, 3.29399768, 0.29040966, 1.69136631, 1.71086031],
       [0.72581411, 3.21806371, 0.36159148, 1.54808703, 1.21567622]])
'''

# np.tile将x沿着给定轴平铺给定次数
np.linalg.norm(np.tile(X_new, (1, k)).reshape(-1, k, 2) - kmeans.cluster_centers_, axis=2)
'''
array([[2.81093633, 0.32995317, 2.9042344 , 1.49439034, 2.88633901],
       [5.80730058, 2.80290755, 5.84739223, 4.4759332 , 5.84236351],
       [1.21475352, 3.29399768, 0.29040966, 1.69136631, 1.71086031],
       [0.72581411, 3.21806371, 0.36159148, 1.54808703, 1.21567622]])
'''

K-Means原理

kmeans_iter1 = KMeans(n_clusters=5, init='random', n_init=1, algorithm='full', max_iter=1, random_state=0)
kmeans_iter2 = KMeans(n_clusters=5, init='random', n_init=1, algorithm='full', max_iter=2, random_state=0)
kmeans_iter3 = KMeans(n_clusters=5, init='random', n_init=1, algorithm='full', max_iter=3, random_state=0)
kmeans_iter1.fit(X)
kmeans_iter2.fit(X)
kmeans_iter3.fit(X)

plt.figure(figsize=(10, 8))

plt.subplot(321)
plot_data(X)
plot_centroids(kmeans_iter1.cluster_centers_, circle_color='r', cross_color='w')
plt.ylabel('$x_2$', fontsize=14, rotation=0)
plt.tick_params(labelbottom=False)
plt.title('Update the centroids (initially randomly)', fontsize=14)

plt.subplot(322)
plot_decision_boundaries(kmeans_iter1, X, show_xlabels=False, show_ylabels=False)
plt.title('Label the instances', fontsize=14)

plt.subplot(323)
plot_decision_boundaries(kmeans_iter1, X, show_centroids=False, show_xlabels=False)
plot_centroids(kmeans_iter2.cluster_centers_)

plt.subplot(324)
plot_decision_boundaries(kmeans_iter2, X, show_xlabels=False, show_ylabels=False)

plt.subplot(325)
plot_decision_boundaries(kmeans_iter2, X, show_centroids=False)
plot_centroids(kmeans_iter3.cluster_centers_)

plt.subplot(326)
plot_decision_boundaries(kmeans_iter3, X, show_ylabels=False)

save_fig('kmeans_algorithm_plot')
plt.show()

K-Means的不稳定性，不同的聚类结果

def plot_clusterer_comparison(clusterer1, clusterer2, X, title1=None, title2=None):
    clusterer1.fit(X)
    clusterer2.fit(X)
    
    plt.figure(figsize=(10, 3.2))
    
    plt.subplot(121)
    plot_decision_boundaries(clusterer1, X)
    if title1:
        plt.title(title1, fontsize=14)
    
    plt.subplot(122)
    plot_decision_boundaries(clusterer2, X, show_ylabels=False)
    if title2:
        plt.title(title2, fontsize=14)

kmeans_rnd_init1 = KMeans(n_clusters=5, init='random', n_init=1, algorithm='full', random_state=2)
kmeans_rnd_init2 = KMeans(n_clusters=5, init='random', n_init=1, algorithm='full', random_state=5)

plot_clusterer_comparison(kmeans_rnd_init1, kmeans_rnd_init2, X, 'Solution 1', 'Solution 2 (with a different random init)')

save_fig('kmeans_variablity_plot')
plt.show()

惯性

kmeans.inertia_  # 211.5985372581683

X_dist = kmeans.transform(X)
np.sum(X_dist[np.arange(len(X_dist)), kmeans.labels_]**2)  # 211.59853725816856

多次执行，保留最好结果

kmeans_rnd_init1.inertia_  # 219.43539442771407

kmeans_rnd_init2.inertia_  # 211.5985372581683

kmeans_rnd_10_inits = KMeans(n_clusters=5, init='random', n_init=10, algorithm='full', random_state=2)
kmeans_rnd_10_inits.fit(X)

plt.figure(figsize=(8, 4))
plot_decision_boundaries(kmeans_rnd_10_inits, X)
plt.show()

K-Means++

KMeans()
good_init = np.array([[-3, 3], [-3, 2], [-3, 1], [-1, 2], [0, 2]])
kmeans = KMeans(n_clusters=5, init=good_init, n_init=1, random_state=42)
kmeans.fit(X)
kmeans.inertia_  # 211.5985372581683

加速KMeans

# 指定循环执行50次
%timeit -n 50 KMeans(algorithm='elkan', random_state=42).fit(X)  # 107 ms ± 15 ms per loop (mean ± std. dev. of 7 runs, 50 loops each)

%timeit -n 50 KMeans(algorithm='full', random_state=42).fit(X)  # 281 ms ± 25.6 ms per loop (mean ± std. dev. of 7 runs, 50 loops each)

小批量KMeans

from sklearn.cluster import MiniBatchKMeans

minibatch_kmeans = MiniBatchKMeans(n_clusters=5, random_state=42)
minibatch_kmeans.fit(X)

minibatch_kmeans.inertia_  # 211.93186531476786

import urllib.request
from sklearn.datasets import fetch_openml

mnist = fetch_openml('mnist_784', version=1, as_frame=False)
mnist.target = mnist.target.astype(np.int64)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(mnist['data'], mnist['target'], random_state=42)

filename = 'my_mnist.data'
X_mm = np.memmap(filename, dtype='float32', mode='write', shape=X_train.shape)
X_mm[:] = X_train

minibatch_kmeans = MiniBatchKMeans(n_clusters=10, batch_size=10, random_state=42)
minibatch_kmeans.fit(X_mm)

def load_next_batch(batch_size):
    # 表示从a中随机选取size个数;
    # replace 代表的意思是抽样之后还放不放回去，如果是False的话，那么每一次挑选出来的数都不一样，
    # 如果是True的话， 有可能会出现重复的，因为前面的抽的放回去了
    return X[np.random.choice(len(X), batch_size, replace=False)]

np.random.seed(42)

k = 5
n_init = 10
n_iterations = 100
batch_size = 100
init_size = 500
evaluate_on_last_n_iters = 10
best_kmeans = None
for init in range(n_init):
    minibatch_kmeans = MiniBatchKMeans(n_clusters=k, init_size=init_size)
    X_init = load_next_batch(init_size)
    minibatch_kmeans.partial_fit(X_init)
    
    minibatch_kmeans.sum_interia_ = 0
    for iteration in range(n_iterations):
        X_batch = load_next_batch(batch_size)
        minibatch_kmeans.partial_fit(X_batch)
        if iteration >= n_iterations - evaluate_on_last_n_iters:
            minibatch_kmeans.sum_interia_ += minibatch_kmeans.inertia_
    if (best_kmeans is None or minibatch_kmeans.sum_interia_ < best_kmeans.sum_interia_):
        best_kmeans = minibatch_kmeans

best_kmeans.score(X)  # -211.63940111556585

%timeit KMeans(n_clusters=5, random_state=42).fit(X)  # 56.7 ms ± 1.69 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit MiniBatchKMeans(n_clusters=5, random_state=42).fit(X)  # 20.9 ms ± 664 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

from timeit import timeit

times = np.empty((100, 2))
inertias = np.empty((100, 2))
for k in range(1, 101):
    kmeans_ = KMeans(n_clusters=k, random_state=42)
    minibatch_kmeans = MiniBatchKMeans(n_clusters=k, random_state=42)
    print('
{}/{}'.format(k, 100), end='')
    times[k-1, 0] = timeit('kmeans_.fit(X)', number=10, globals=globals())
    times[k-1, 1] = timeit('minibatch_kmeans.fit(X)', number=10, globals=globals())
    inertias[k-1, 0] = kmeans_.inertia_
    inertias[k-1, 1] = minibatch_kmeans.inertia_

plt.figure(figsize=(10, 4))

plt.subplot(121)
plt.plot(range(1, 101), inertias[:, 0], 'r--', label='K-Means')
plt.plot(range(1, 101), inertias[:, 1], 'b.-', label='Mini-batch K-Means')
plt.xlabel('$k$', fontsize=16)
plt.title('Inertia', fontsize=14)
plt.legend(fontsize=14)
plt.axis([1, 100, 0, 100])

plt.subplot(122)
plt.plot(range(1, 101), times[:, 0], 'r--', label='K-Means')
plt.plot(range(1, 101), times[:, 1], 'b.-', label='Mini-batch K-Means')
plt.xlabel('$k$', fontsize=16)
plt.title('Training time (seconds)', fontsize=14)
plt.axis([1, 100, 0, 6])

save_fig('minibatch_kmeans_vs_kmeans')
plt.show()

寻找最优类数

kmeans_k3 = KMeans(n_clusters=3, random_state=42)
kmeans_k8 = KMeans(n_clusters=8, random_state=42)
plot_clusterer_comparison(kmeans_k3, kmeans_k8, X, '$k=3$', '$k=8$')
save_fig('bad_n_clusters_plot')
plt.show()

kmeans_k3.inertia_  # 653.2167190021556

kmeans_k8.inertia_  # 119.1198341610289

kmeans_per_k = [KMeans(n_clusters=k, random_state=42).fit(X) for k in range(1, 10)]
inertias = [model.inertia_ for model in kmeans_per_k]

plt.figure(figsize=(8, 3.5))
plt.plot(range(1, 10), inertias, 'bo-')
plt.xlabel('$k$', fontsize=14)
plt.ylabel('Inertia', fontsize=14)
plt.annotate('Elbow', xy=(4, inertias[3]), xytext=(0.55, 0.55), textcoords='figure fraction',
            fontsize=16, arrowprops=dict(facecolor='black', shrink=0.1))
plt.axis([1, 8.5, 0, 1300])
save_fig('inertia_vs_k_plot')
plt.show()

plot_decision_boundaries(kmeans_per_k[4-1], X)
plt.show()

轮廓系数

from sklearn.metrics import silhouette_score

silhouette_score(X, kmeans.labels_)  # 0.655517642572828

silhouette_scores = [silhouette_score(X, model.labels_) for model in kmeans_per_k[1:]]

plt.figure(figsize=(8, 3))
plt.plot(range(2, 10), silhouette_scores, 'bo-')
plt.xlabel('$k$', fontsize=14)
plt.ylabel('Silhouette score', fontsize=14)
plt.axis([1.8, 8.5, 0.55, 0.7])
save_fig('silhouette_score_vs_k_plot')
plt.show()

from sklearn.metrics import silhouette_samples
from matplotlib.ticker import FixedLocator, FixedFormatter

plt.figure(figsize=(11, 9))
for k in (3, 4, 5 ,6):
    plt.subplot(2, 2, k-2)
    y_pred = kmeans_per_k[k-1].labels_
    silhouette_coefficients = silhouette_samples(X, y_pred)
    padding = len(X) // 30
    pos = padding
    ticks = []
    for i in range(k):
        coeffs = silhouette_coefficients[y_pred == i]
        coeffs.sort()
        
        color = mpl.cm.Spectral(i / k)
        plt.fill_betweenx(np.arange(pos, pos + len(coeffs)), 0, coeffs, facecolor=color, edgecolor=color, alpha=0.7)
        ticks.append(pos + len(coeffs) // 2)
        pos += len(coeffs) + padding
        
    plt.gca().yaxis.set_major_locator(FixedLocator(ticks))
    plt.gca().yaxis.set_major_formatter(FixedFormatter(range(k)))
    if k in (3, 5):
        plt.ylabel('Cluster')
    if k in (5, 6):
        plt.gca().set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
        plt.xlabel('Silhouette Coefficient')
    else:
        plt.tick_params(labelbottom=False)
    plt.axvline(x=silhouette_scores[k - 2], color='red', linestyle='--')
    plt.title('$k={}$'.format(k), fontsize=16)

save_fig('silhouette_analysis_plot')
plt.show()

K-Means的局限

X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y = np.r_[y1, y2]

plot_clusters(X, y)

kmeans_good = KMeans(n_clusters=3, init=np.array([[-1.5, 2.5], [0.5, 0], [4, 0]]), n_init=1, random_state=42)
kmeans_bad = KMeans(n_clusters=3, random_state=42)
kmeans_good.fit(X)
kmeans_bad.fit(X)

plt.figure(figsize=(10, 3.2))

plt.subplot(121)
plot_decision_boundaries(kmeans_good, X)
plt.title('Inertia={:.1f}'.format(kmeans_good.inertia_), fontsize=14)

plt.subplot(122)
plot_decision_boundaries(kmeans_bad, X, show_ylabels=False)
plt.title('Inertia={:.1f}'.format(kmeans_bad.inertia_), fontsize=14)

save_fig('bad_kmeans_plot')
plt.show()

利用聚类分割图像

# Download the ladybug image
images_path = os.path.join(PROJECT_ROOT_DIR, "images", "unsupervised_learning")
os.makedirs(images_path, exist_ok=True)
DOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml2/master/"
filename = "ladybug.png"
print("Downloading", filename)
url = DOWNLOAD_ROOT + "images/unsupervised_learning/" + filename
urllib.request.urlretrieve(url, os.path.join(images_path, filename))

from matplotlib.image import imread

image = imread(os.path.join(images_path, filename))
image.shape  # (533, 800, 3)

X = image.reshape(-1, 3)
kmeans = KMeans(n_clusters=8, random_state=42).fit(X)
segmented_img = kmeans.cluster_centers_[kmeans.labels_]
segmented_img = segmented_img.reshape(image.shape)

segmented_imgs = []
n_colors = (10, 8, 6, 4, 2)
for n_clusters in n_colors:
    kmeans = KMeans(n_clusters=n_clusters, random_state=42).fit(X)
    segmented_img = kmeans.cluster_centers_[kmeans.labels_]
    segmented_imgs.append(segmented_img.reshape(image.shape))

plt.figure(figsize=(10, 5))
plt.subplots_adjust(wspace=0.05, hspace=0.1)

plt.subplot(231)
plt.imshow(image)
plt.title('Original image')
plt.axis('off')

for idx, n_clusters in enumerate(n_colors):
    plt.subplot(232 + idx)
    plt.imshow(segmented_imgs[idx])
    plt.title('{} colors'.format(n_clusters))
    plt.axis('off')

save_fig('image_segmentation_diagram', tight_layout=False)
plt.show()

利用聚类进行预处理

from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

X_digits, y_digits = load_digits(return_X_y = True)
X_train, X_test, y_train, y_test = train_test_split(X_digits, y_digits, random_state=42)
log_reg = LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=5000, random_state=42)
log_reg.fit(X_train, y_train)
log_reg_score = log_reg.score(X_test, y_test)
log_reg_score  # 0.9688888888888889

from sklearn.pipeline import Pipeline

pipeline = Pipeline([
    ('kmeans', KMeans(n_clusters=50, random_state=42)),
    ('log_reg', LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=5000, random_state=42))
])
pipeline.fit(X_train, y_train)

pipeline_score = pipeline.score(X_test, y_test)
pipeline_score  # 0.9777777777777777

1 - (1 - pipeline_score) / (1 - log_reg_score)  # 0.28571428571428414

from sklearn.model_selection import GridSearchCV

param_grid = dict(kmeans__n_clusters=range(2, 100))
grid_clf = GridSearchCV(pipeline, param_grid, cv=3, verbose=2)
grid_clf.fit(X_train, y_train)

grid_clf.best_params_  # {'kmeans__n_clusters': 88}

grid_clf.score(X_test, y_test)  # 0.9822222222222222

聚类用于半监督学习

n_labeled = 50
log_reg = LogisticRegression(multi_class='ovr', solver='lbfgs', random_state=42)
log_reg.fit(X_train[:n_labeled], y_train[:n_labeled])
log_reg.score(X_test, y_test)  # 0.8333333333333334

k = 50
kmeans = KMeans(n_clusters=k, random_state=42)
X_digits_dist = kmeans.fit_transform(X_train)
representative_digit_idx = np.argmin(X_digits_dist, axis=0)
X_representativ_digits = X_train[representative_digit_idx]

plt.figure(figsize=(8, 2))
for index, X_representative_digit in enumerate(X_representativ_digits):
    plt.subplot(k // 10, 10, index + 1)
    plt.imshow(X_representative_digit.reshape(8, 8), cmap='binary', interpolation='bilinear')
    plt.axis('off')

save_fig('representative_iamges_diagram', tight_layout=False)
plt.show()

y_train[representative_digit_idx]
'''
array([4, 8, 0, 6, 8, 3, 7, 7, 9, 2, 5, 5, 8, 5, 2, 1, 2, 9, 6, 1, 1, 6,
       9, 0, 8, 3, 0, 7, 4, 1, 6, 5, 2, 4, 1, 8, 6, 3, 9, 2, 4, 2, 9, 4,
       7, 6, 2, 3, 1, 1])
'''

y_representative_digits = np.array([
    4, 8, 0, 6, 8, 3, 7, 7, 9, 2,
    5, 5, 8, 5, 2, 1, 2, 9, 6, 1,
    1, 6, 9, 0, 8, 3, 0, 7, 4, 1, 
    6, 5, 2, 4, 1, 8, 6, 3, 9, 2,
    4, 2, 9, 4, 7, 6, 2, 3, 1, 1
    ])

log_reg = LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=5000, random_state=42)
log_reg.fit(X_representativ_digits, y_representative_digits)
log_reg.score(X_test, y_test)  # 0.9222222222222223

y_train_propagated = np.empty(len(X_train), dtype=np.int32)
for i in range(k):
    y_train_propagated[kmeans.labels_==i] = y_representative_digits[i]

log_reg = LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=5000, random_state=42)
log_reg.fit(X_train, y_train_propagated)

log_reg.score(X_test, y_test)  # 0.9333333333333333

percentile_closest = 75

X_cluster_dist = X_digits_dist[np.arange(len(X_train)), kmeans.labels_]
for i in range(k):
    in_cluster = (kmeans.labels_ == i)
    cluster_dist = X_cluster_dist[in_cluster]
    cutoff_distance = np.percentile(cluster_dist, percentile_closest)
    above_cutoff = (X_cluster_dist > cutoff_distance)
    X_cluster_dist[in_cluster & above_cutoff] == -1

partially_propagated = (X_cluster_dist != -1)
X_train_partially_propagated = X_train[partially_propagated]
y_train_partially_propagated = y_train_propagated[partially_propagated]

log_reg = LogisticRegression(multi_class='ovr', solver='lbfgs', max_iter=5000, random_state=42)
log_reg.fit(X_train_partially_propagated, y_train_partially_propagated)

log_reg.score(X_test, y_test)  # 0.9333333333333333

np.mean(y_train_partially_propagated == y_train[partially_propagated])  # 0.9391239792130661

DBSCAN

from sklearn.datasets import make_moons

X, y = make_moons(n_samples=1000, noise=0.05, random_state=42)

from sklearn.cluster import DBSCAN

dbscan = DBSCAN(eps=0.05, min_samples=5)
dbscan.fit(X)
dbscan.labels_[:10]  # array([ 0,  2, -1, -1,  1,  0,  0,  0,  2,  5], dtype=int64)

len(dbscan.core_sample_indices_)  # 808

dbscan.core_sample_indices_[:10]  # array([ 0,  4,  5,  6,  7,  8, 10, 11, 12, 13], dtype=int64)

dbscan.components_[:3]
'''
array([[-0.02137124,  0.40618608],
       [-0.84192557,  0.53058695],
       [ 0.58930337, -0.32137599]])
'''

np.unique(dbscan.labels_)  # array([-1,  0,  1,  2,  3,  4,  5,  6], dtype=int64)

dbscan2 = DBSCAN(eps=0.2)
dbscan2.fit(X)

def plot_dbscan(dbscan, X, size, show_xlabels=True, show_ylabels=True):
    # 实现构造一个矩阵x_new，其维度与矩阵x一致，并为其初始化为全0
    core_mask = np.zeros_like(dbscan.labels_, dtype=bool)
    core_mask[dbscan.core_sample_indices_] = True
    anomalies_mask = dbscan.labels_ == -1
    non_core_mask = ~(core_mask | anomalies_mask)
    
    cores = dbscan.components_
    anomalies = X[anomalies_mask]
    non_cores = X[non_core_mask]
    
    plt.scatter(cores[:, 0], cores[:, 1], c=dbscan.labels_[core_mask], marker='o', s=size, cmap='Paired')
    plt.scatter(cores[:, 0], cores[:, 1], marker='*', s=20, c=dbscan.labels_[core_mask])
    plt.scatter(anomalies[:, 0], anomalies[:, 1], c='r', marker='x', s=100)
    plt.scatter(non_cores[:, 0], non_cores[:, 1], c=dbscan.labels_[non_core_mask], marker='.')
    if show_xlabels:
        plt.xlabel('$x_1$', fontsize=14)
    else:
        plt.tick_params(labelbotto=False)
    if show_ylabels:
        plt.ylabel('$x_2$', fontsize=14, rotation=0)
    else:
        plt.tick_params(labelleft=False)
    plt.title('eps={:.2f}, min_samples={}'.format(dbscan.eps, dbscan.min_samples), fontsize=14)

plt.figure(figsize=(9, 3.2))
plt.subplot(121)
plot_dbscan(dbscan, X, size=100)
plt.subplot(122)
plot_dbscan(dbscan2, X, size=600, show_ylabels=False)

save_fig('dbscan_plot')
plt.show()

from sklearn.neighbors import KNeighborsClassifier

dbscan = dbscan2
knn = KNeighborsClassifier(n_neighbors=50)
knn.fit(dbscan.components_, dbscan.labels_[dbscan.core_sample_indices_])

X_new = np.array([[-0.5, 0], [0, 0.5], [1, -0.1], [2, 1]])
knn.predict(X_new)  # array([1, 0, 1, 0], dtype=int64)

knn.predict_proba(X_new)
'''
array([[0.18, 0.82],
       [1.  , 0.  ],
       [0.12, 0.88],
       [1.  , 0.  ]])
'''

plt.figure(figsize=(6, 3))
plot_decision_boundaries(knn, X, show_centroids=False)
plt.scatter(X_new[:, 0], X_new[:, 1], c='b', marker='+', s=200, zorder=10)
save_fig('cluster_classification_plot')
plt.show()

y_dist, y_pred_idx = knn.kneighbors(X_new, n_neighbors=1)
y_pred = dbscan.labels_[dbscan.core_sample_indices_][y_pred_idx]
y_pred[y_dist > 0.2] = -1
y_pred.ravel()  # array([-1,  0,  1, -1], dtype=int64)

其他聚类算法

# Spectral Clustering
from sklearn.cluster import SpectralClustering

sc1 = SpectralClustering(n_clusters=2, gamma=100, random_state=42)
sc1.fit(X)
sc2 = SpectralClustering(n_clusters=2, gamma=1, random_state=42)
sc2.fit(X)
np.percentile(sc1.affinity_matrix_, 95)  # 0.04251990648936265

def plot_spectral_clustering(sc, X, size, alpha, show_xlables=True, show_ylabels=True):
    plt.scatter(X[:, 0], X[:, 1], marker='o', s=size, c='gray', cmap='Paired', alpha=alpha)
    plt.scatter(X[:, 0], X[:, 1], marker='o', s=30, c='w')
    plt.scatter(X[:, 0], X[:, 1], marker='.', s=10, c=sc.labels_, cmap='Paired')
    if show_xlables:
        plt.xlabel('$x_1$', fontsize=14)
    else:
        plt.tick_params(labelbottom=False)
    if show_ylabels:
        plt.ylabel('$x_2$', fontsize=14, rotation=0)
    else:
        plt.tick_params(labelleft=False)
    plt.title('RBF gamma={}'.format(sc.gamma), fontsize=14)

plt.figure(figsize=(9, 3.2))
plt.subplot(121)
plot_spectral_clustering(sc1, X, size=500, alpha=0.1)

plt.subplot(122)
plot_spectral_clustering(sc2, X, size=4000, alpha=0.01, show_ylabels=False)

plt.show()

# Agglomerative Clustering
from sklearn.cluster import AgglomerativeClustering

X = np.array([0, 2, 5, 8.5]).reshape(-1, 1)
agg = AgglomerativeClustering(linkage='complete').fit(X)

def learned_parameters(estimator):
    return [attrib for attrib in dir(estimator) if attrib.endswith('_') and not attrib.startswith('_')]

learned_parameters(agg)
'''
['children_',
 'labels_',
 'n_clusters_',
 'n_connected_components_',
 'n_features_in_',
 'n_leaves_']
'''

agg.children_
'''
array([[0, 1],
       [2, 3],
       [4, 5]])
'''

高斯混合模型

X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y= np.r_[y1, y2]

from sklearn.mixture import GaussianMixture

gm = GaussianMixture(n_components=3, n_init=10, random_state=42)
gm.fit(X)

gm.weights_  # array([0.39025715, 0.40007391, 0.20966893])

gm.means_
'''
array([[ 0.05131611,  0.07521837],
       [-1.40763156,  1.42708225],
       [ 3.39893794,  1.05928897]])
'''

gm.covariances_
'''
array([[[ 0.68799922,  0.79606357],
        [ 0.79606357,  1.21236106]],

       [[ 0.63479409,  0.72970799],
        [ 0.72970799,  1.1610351 ]],

       [[ 1.14833585, -0.03256179],
        [-0.03256179,  0.95490931]]])
'''

gm.converged_  # True

gm.n_iter_  $ 

gm.predict(X)  # array([0, 0, 1, ..., 2, 2, 2], dtype=int64)

gm.predict_proba(X)
'''
array([[9.76741808e-01, 6.78581203e-07, 2.32575136e-02],
       [9.82832955e-01, 6.76173663e-04, 1.64908714e-02],
       [7.46494398e-05, 9.99923327e-01, 2.02398402e-06],
       ...,
       [4.26050456e-07, 2.15512941e-26, 9.99999574e-01],
       [5.04987704e-16, 1.48083217e-41, 1.00000000e+00],
       [2.24602826e-15, 8.11457779e-41, 1.00000000e+00]])
'''

X_new, y_new = gm.sample(6)
X_new
'''
array([[-0.86944074, -0.32767626],
       [ 0.29836051,  0.28297011],
       [-2.8014927 , -0.09047309],
       [ 3.98203732,  1.49951491],
       [ 3.81677148,  0.53095244],
       [ 2.84104923, -0.73858639]])
'''

y_new  # array([0, 0, 1, 2, 2, 2])

gm.score_samples(X)
'''
array([-2.60768954, -3.57110232, -3.32987086, ..., -3.51347241,
       -4.39798588, -3.80746532])
'''

resolution = 100
grid = np.arange(-10, 10, 1 / resolution)
xx, yy = np.meshgrid(grid, grid)
X_full = np.vstack([xx.ravel(), yy.ravel()]).T
pdf = np.exp(gm.score_samples(X_full))
pdf_probas = pdf * (1 / resolution) ** 2
pdf_probas.sum()  # 0.9999999999215021

from matplotlib.colors import LogNorm

def plot_gaussian_mixture(clusterer, X, resolution=1000, show_ylabels=True):
    mins = X.min(axis=0) - 0.1
    maxs = X.max(axis=0) + 0.1
    xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution), np.linspace(mins[1], maxs[1], resolution))
    Z = -clusterer.score_samples(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    plt.contourf(xx, yy, Z, norm=LogNorm(vmin=1.0, vmax=30.0), levels=np.logspace(0, 2, 12))
    plt.contour(xx, yy, Z, norm=LogNorm(vmin=1.0, vmax=30.0), levels=np.logspace(0, 2, 12), linewidths=1, colors='k')
    Z = clusterer.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    plt.contour(xx, yy, Z, linewidths=2, colors='r', linestyle='dashed')
    plt.plot(X[:, 0], X[:, 1], 'k.', markersize=2)
    plot_centroids(clusterer.means_, clusterer.weights_)
    plt.xlabel('$x_1$', fontsize=14)
    if show_ylabels:
        plt.ylabel('$x_2$', fontsize=14, rotation=0)
    else:
        plt.tick_params(labelleft=False)

plt.figure(figsize=(8, 4))
plot_gaussian_mixture(gm, X)
save_fig('gaussian_mixatures_plot')
plt.show()

gm_full = GaussianMixture(n_components=3, n_init=10, covariance_type='full', random_state=42)
gm_tied = GaussianMixture(n_components=3, n_init=10, covariance_type='tied', random_state=42)
gm_spherical = GaussianMixture(n_components=3, n_init=10, covariance_type='spherical', random_state=42)
gm_diag = GaussianMixture(n_components=3, n_init=10, covariance_type='diag', random_state=42)
gm_full.fit(X)
gm_tied.fit(X)
gm_spherical.fit(X)
gm_diag.fit(X)

def comare_gaussian_mixtures(gm1, gm2, X):
    plt.figure(figsize=(9, 4))
    plt.subplot(121)
    plot_gaussian_mixture(gm1, X)
    plt.title('covariance_type="{}"'.format(gm1.covariance_type), fontsize=14)
    
    plt.subplot(122)
    plot_gaussian_mixture(gm2, X, show_ylabels=False)
    plt.title('covariance_type="{}"'.format(gm2.covariance_type), fontsize=14)

comare_gaussian_mixtures(gm_tied, gm_spherical, X)
save_fig('covariance_type_plot')
plt.show()

comare_gaussian_mixtures(gm_full, gm_diag, X)
plt.tight_layout()
plt.show()

利用高斯混合模型进行异常检测

densities = gm.score_samples(X)
density_threshold = np.percentile(densities, 4)
anomalies = X[densities < density_threshold]

plt.figure(figsize=(8, 4))
plot_gaussian_mixture(gm, X)
plt.scatter(anomalies[:, 0], anomalies[:, 1], color='r', marker='*')
plt.ylim(top=5.1)
save_fig('mixture_anomaly_detection_plot')
plt.show()

模型选择

gm.bic(X)  # 8189.747000497186

gm.aic(X)  # 8102.521720382148

n_clusters = 3
n_dims = 2
n_params_for_weights = n_clusters - 1
n_params_for_means = n_clusters * n_dims
n_params_for_covariance = n_clusters * n_dims * (n_dims + 1) // 2
n_params = n_params_for_weights + n_params_for_means + n_params_for_covariance
max_log_likelihood = gm.score(X) * len(X)
bic = np.log(len(X)) * n_params - 2 * max_log_likelihood
aic = 2 * n_params - 2 * max_log_likelihood

bic, aic  # (8189.747000497186, 8102.521720382148)

n_params  # 17

gms_per_k = [GaussianMixture(n_components=k, n_init=10, random_state=42).fit(X) for k in range(1, 11)]

bics = [model.bic(X) for model in gms_per_k]
aics = [model.aic(X) for model in gms_per_k]

plt.figure(figsize=(8, 3))
plt.plot(range(1, 11), bics, 'bo-', label='BIC')
plt.plot(range(1, 11), aics, 'gp--', label='AIC')
plt.xlabel('$k$', fontsize=14)
plt.ylabel('Information Criterion', fontsize=14)
plt.axis([1, 9.5, np.min(aics) - 50, np.max(aics) + 50])
plt.annotate('Minimum', xy=(3, bics[2]), xytext=(0.35, 0.6), textcoords='figure fraction', fontsize=14, arrowprops=dict(facecolor='black', shrink=0.1))
plt.legend()
save_fig('aic_bic_vs_k_plot')
plt.show()

min_bic = np.infty
for k in range(1, 11):
    for covariance_type in ('full', 'tied', 'spherical', 'diag'):
        bic = GaussianMixture(n_components=k, n_init=10, covariance_type=covariance_type, random_state=42).fit(X).bic(X)
        if bic < min_bic:
            min_bic = bic
            best_k = k
            best_covariance_type = covariance_type

best_k  # 3

best_covariance_type  # 'full'

变分贝叶斯高斯混合模型

from sklearn.mixture import BayesianGaussianMixture

bgm = BayesianGaussianMixture(n_components=10, n_init=10, random_state=42)
bgm.fit(X)

np.round(bgm.weights_, 2)  # array([0.4 , 0.21, 0.4 , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ])

plt.figure(figsize=(8, 5))
plot_gaussian_mixture(bgm, X)
plt.show()

bgm_low = BayesianGaussianMixture(n_components=10, max_iter=1000, n_init=1, weight_concentration_prior=0.01, random_state=42)
bgm_high = BayesianGaussianMixture(n_components=10, max_iter=1000, n_init=1, weight_concentration_prior=10000, random_state=42)
nn = 73
bgm_low.fit(X[:nn])
bgm_high.fit(X[:nn])

np.round(bgm_low.weights_, 2)  # array([0.52, 0.48, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ])

np.round(bgm_high.weights_, 2)  # array([0.01, 0.18, 0.27, 0.11, 0.01, 0.01, 0.01, 0.01, 0.37, 0.01])

plt.figure(figsize=(9, 4))
plt.subplot(121)
plot_gaussian_mixture(bgm_low, X[:nn])
plt.title('weight_concentration_prior=0.01', fontsize=14)
plt.subplot(122)
plot_gaussian_mixture(bgm_high, X[:nn], show_ylabels=False)
plt.title('weight_concentration_prior=10000', fontsize=14)
save_fig('mixture_concentration_prior_plot')
plt.show()

X_moons, y_moons = make_moons(n_samples=1000, noise=0.05, random_state=42)
bgm = BayesianGaussianMixture(n_components=10, n_init=10, random_state=42)
bgm.fit(X_moons)

plt.figure(figsize=(9, 3))
plt.subplot(121)
plot_data(X_moons)
plt.xlabel('$x_1$', fontsize=14)
plt.ylabel('$x_2$', fontsize=14, rotation=0)
plt.subplot(122)
plot_gaussian_mixture(bgm, X_moons, show_ylabels=False)
save_fig('moons_vs_bgm_plot')
plt.show()

似然函数

from scipy.stats import norm

xx = np.linspace(-6, 4, 101)
ss = np.linspace(1, 2, 101)
XX, SS = np.meshgrid(xx, ss)
ZZ =  2 * norm.pdf(XX - 1.0, 0, SS) + norm.pdf(XX + 4.0, 0, SS)
ZZ = ZZ / ZZ.sum(axis=1)[:, np.newaxis] / (xx[1] - xx[0])

from matplotlib.patches import Polygon

plt.figure(figsize=(8, 4.5))
x_idx = 85
s_idx = 30

plt.subplot(221)
plt.contourf(XX, SS, ZZ, cmap='GnBu')
plt.plot([-6, 4], [ss[s_idx], ss[s_idx]], 'k-', linewidth=2)
plt.plot([xx[x_idx], xx[x_idx]], [1, 2], 'b-', linewidth=2)
plt.xlabel(r'$x$')
plt.ylabel(r'$	heta$', fontsize=14, rotation=0)
plt.title(r'Model $f(x; 	heta)$', fontsize=14)

plt.subplot(222)
plt.plot(ss, ZZ[:, x_idx], 'b-')
max_idx = np.argmax(ZZ[:, x_idx])
max_val = np.max(ZZ[:, x_idx])
plt.plot([ss[max_idx], ss[max_idx]], [0, max_val], 'r:')
plt.plot([0, ss[max_idx]], [max_val, max_val], 'r:')
plt.text(1.01, max_val + 0.005, r'$hat {L}$', fontsize=14)
plt.text(ss[max_idx] + 0.01, 0.055, r'$hat {	heta}$', fontsize=14)
plt.text(ss[max_idx] + 0.01, max_val - 0.02, r'$Max$', fontsize=12)
plt.axis([1, 2, 0.05, 0.15])
plt.xlabel(r'$	heta$', fontsize=14)
plt.grid(True)
plt.text(1.99, 0.135, r'$=f(x=2.5;	heta)$', fontsize=14, ha='right')
plt.title(r'Likelihood function $mathcal{L}(	heta|x=2.5)$', fontsize=14)

plt.subplot(223)
plt.plot(xx, ZZ[s_idx], 'k-')
plt.axis([-6, 4, 0, 0.025])
plt.xlabel(r'$x$', fontsize=14)
plt.grid(True)
plt.title(r'PDF $f(x;	heta=1.3)$', fontsize=14)
verts = [(xx[41], 0)] + list(zip(xx[41:81], ZZ[s_idx, 41:81])) + [(xx[80], 0)]
poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')
plt.gca().add_patch(poly)

plt.subplot(224)
plt.plot(ss, np.log(ZZ[:, x_idx]), 'b-')
max_idx = np.argmax(np.log(ZZ[:, x_idx]))
max_val = np.max(np.log(ZZ[:, x_idx]))
plt.plot(ss[max_idx], max_val, 'r.')
plt.plot([ss[max_idx], ss[max_idx]], [-5, max_val], 'r.')
plt.plot([0, ss[max_idx]], [max_val, max_val], 'r:')
plt.axis([1, 2, -2.4, -2])
plt.xlabel(r'$	heta$', fontsize=14)
plt.text(ss[max_idx] + 0.01, max_val - 0.05, r'$Max$', fontsize=12)
plt.text(ss[max_idx] + 0.01, -2.39, r'$hat {	heta}$', fontsize=14)
plt.text(1.01, max_val + 0.02, r'$log \, hat{L}$', fontsize=14)
plt.grid(True)
plt.title(r'$log \, mathcal{L}(	heta|x=2.5)$', fontsize=14)

save_fig('likelihood_function_plot')
plt.show()