Python-无监督学习实用手册-二-Python 无监督学习实用手册（二）五、软聚类和高斯混合模型在本章中，我们将

Python 无监督学习实用手册（二）

原文：Hands-on unsupervised learning with Python

协议：CC BY-NC-SA 4.0

五、软聚类和高斯混合模型

在本章中，我们将讨论软聚类的概念，它允许我们获得数据集的每个样本相对于定义的聚类配置的隶属度。也就是说，考虑从 0%到 100%的范围，我们想知道x_IT3 在多大程度上属于一个集群。极值为 0，表示 x _i 完全在集群的域外，为 1 (100%)，表示 x _i 完全分配给单个集群。所有中间值都意味着两个或多个不同簇的部分区域。因此，与硬聚类相反，在这里，我们感兴趣的不是确定一个固定的赋值，而是一个具有概率分布(或概率本身)相同属性的向量。这种方法允许对边界样本进行更好的控制，并帮助我们找到生成过程的合适近似，数据集应该从该生成过程中提取。

特别是，我们将讨论以下主题:

模糊 C 均值
高斯混合
AIC 和 BIC 作为绩效指标
贝叶斯高斯混合(简介)
生成(半监督)高斯混合

技术要求

本章将介绍的代码需要以下内容:

Python 3.5+ ( Anaconda 发行版强烈推荐)
以下库:
- SciPy 0.19+
- NumPy 1.10+
- 学习 0.20+
- Scikit-fuzzy 0.2
- 熊猫 0.22+
- Matplotlib 2.0+
- seaborn 0.9+

这些例子可以在 GitHub 资源库中找到，网址为 https://GitHub . com/packktpublishing/HandsOn-Unsupervised-Learning-with-Python/tree/master/chapter 05。

软聚类

在第 4 章、行动中的层次聚类中分析的所有算法都属于硬聚类方法家族。这意味着一个给定的样本总是被分配给一个单独的聚类。另一方面，软聚类的目的是将每个样本， x _i 与一个向量相关联，一般表示 x _i 属于每个聚类的概率:

或者，输出可以解释为隶属向量:

在形式上，这两个版本没有区别，但通常情况下，当算法不是明确基于概率分布时，采用后者。然而，出于我们的目的，我们总是将 c(x _i ) 与概率联系起来。这样，读者就会被激励去思考用来获取数据集的数据生成过程。一个明显的例子是将这些向量解释为与构成数据生成过程近似值的特定贡献相关的概率， p _data 。例如，采用概率混合，我们可以决定近似 p _数据 如下:

因此，该过程被分解成(独立的)分量的加权和，输出是每一个分量的 x _i 的概率。当然，我们通常期望每个样本都有一个主导成分，但是通过这种方法，我们对所有边界点都有很大的了解，这些边界点会受到小的扰动，可以分配给不同的聚类。由于这个原因，当输出可以被馈送到另一个可以利用整个概率向量的模型(例如，神经网络)中时，软聚类非常有帮助。例如，推荐者可以首先使用软聚类算法来分割用户，然后处理向量，以便基于显式反馈找到更复杂的关系。一个常见的场景是通过回答诸如“这个结果与你相关吗？”或者“你想看到更多像这样的结果吗？”因为答案是由用户直接提供的，所以它们可以被用在监督或强化学习模型中，该模型的输入基于软自动分割(例如，基于购买历史或详细的页面视图)。通过这种方式，可以通过改变原始分配的效果来轻松管理边缘用户(由于不同集群提供的大量贡献，这可能完全不重要)，同时可以稍微修改对其他具有强成员资格(例如，概率接近 1)的用户的推荐，以提高他们的回报。

我们现在可以开始讨论 Fuzzy c-means，这是一种非常灵活的算法，它将针对 k-means 讨论的概念扩展到软聚类场景。

模糊 c 均值

我们将提出的第一个算法是基于软分配的 k-means 的变体。 Fuzzy c-means 这个名字来源于一个模糊集的概念，它是经典二进制集(即在这种情况下，一个样本可以属于一个单一的聚类)基于代表整个集合不同区域的不同子集的叠加而扩展到集合。例如，基于一些用户年龄的集合可以具有与三个不同(并且部分重叠)年龄范围相关联的程度young、adult和senior:18-35 岁、28-60 岁和> 50 岁。因此，例如，一个 30 岁的用户在不同程度上既是young又是adult(考虑到界限，确实是一个边缘用户)。关于这类集合和所有相关运算的更多细节，我建议阅读《概念和模糊逻辑》一书，麻省理工学院出版社，2011 年。出于我们的目的，我们可以将包含 m 个样本的数据集 X 划分为 k 个重叠的聚类，这样每个样本总是与每个聚类相关联，根据隶属度，w_ijT18】(该值绑定在 0 和 1 之间)。如果 w _ij = 0 ，则表示 x _i 完全在集群C_jT30】之外，反之，则w_ijT35】= 1表示对集群C_jT40】的硬分配。所有中间值代表部分成员资格。当然，由于显而易见的原因，样本的所有隶属度之和必须归一化为 1(类似于概率分布)。这样，一个样本总是属于所有聚类的并集，而将一个聚类分成两个或更多的子聚类总是会产生一个连贯的结果，就成员而言。**

该算法基于广义惯性的优化，S_fT3:

在之前的公式中，μ_jT3】是星团C_jT7】的质心，而 m (m > 1) 是重加权指数系数。当 m ≈ 1 时，重量不受影响。数值越大，如w_ij∑(0，1) ，其重要性成比例降低。可以选择这样的系数来比较不同值的结果和期望的模糊程度。事实上，在每次迭代之后(完全等同于 k 均值)，使用以下公式更新权重:**

如果 x _i 接近质心， μ _j ，则和变得接近 0，权重递增(当然，为了避免数值不稳定，在分母上加了一个小常数，所以永远不能等于 0)。当 m > > 1 时，指数变得接近 0，总和中的所有项都趋于 1。这意味着对特定集群的偏好减弱，并且 w _ij ≈ 1/k 对应于均匀分布。因此，更大的 m 意味着更平坦的划分，在不同的分配之间没有明显的差异(除非样本非常接近质心)，而当 m ≈ 1 时，单个主导权重将几乎等于 1，而其他权重将接近 0(即分配是困难的)。

质心以类似于 k 均值的方式更新(换句话说，目标是最大化分离和内部凝聚力):

重复该过程，直到质心和权重变得稳定。收敛后，可以用一个特定的度量来评估结果，称为归一化的邓恩分割系数，定义如下:

这样的系数被限制在 0 和 1 之间。当 P _C ≈ 0 时，表示 w _C ≈ 1/k ，暗示分布平坦，模糊程度高。另一方面，当 P _C ≈ 1 时，则 w _C ≈ 1 ，表示几乎是硬性作业。所有其他值都与模糊程度成正比。因此，给定一项任务，数据科学家可以根据期望的结果立即评估算法的执行情况。在某些情况下，硬分配更可取，因此， P _C 可以被认为是在切换到例如标准 k 均值之前要执行的检查。其实当 P _C ≈ 1 (而且这样的结果是预期的结果)的时候，再用 Fuzzy c-means 就没有意义了。相反，小于 1 的值(例如， P _C = 0.5 )告知我们，由于存在许多边界样本，可能会有非常不稳定的硬性分配。

现在，让我们将模糊 c 均值算法应用于 scikit-learn 提供的简化 MNIST 数据集。该算法由 Scikit-Fuzzy 库提供，它实现了所有最重要的模糊逻辑模型。第一步是加载和标准化样本，如下所示:

from sklearn.datasets import load_digits

digits = load_digits()
X = digits['data'] / 255.0
Y = digits['target']

X数组包含 1，797 个展平样本， x ∈ ℜ ⁶⁴ ，对应灰度 8 × 8 图像(其值在 0 和 1 之间归一化)。我们要分析不同 m 系数(1.05 到 1.5 之间的 5 个统一值)的行为，并检查样本的权重(在我们的例子中，我们将使用X[0])。因此，我们调用 Scikit-Fuzzy cmeans函数，设置c=10(簇的数量)和两个收敛参数，error=1e-6和maxiter=20000。而且，出于重现性的原因，我们也会随机设置一个标准seed=1000。输入数组应包含列形式的样本；因此，我们需要转置它，如下所示:

from skfuzzy.cluster import cmeans

Ws = []
pcs = []

for m in np.linspace(1.05, 1.5, 5):
    fc, W, _, _, _, _, pc = cmeans(X.T, c=10, m=m, error=1e-6, maxiter=20000, seed=1000)
    Ws.append(W)
    pcs.append(pc)

前面的片段执行不同类型的聚类，并将相应的权重矩阵W和划分系数pc附加到两个列表中。在分析特定配置之前，显示测试样本(代表数字 0)的最终重量(对应于每个数字)会很有帮助:

Weights (in an inverse logarithmic scale) for the sample X[0], corresponding to different m values

由于极值往往差异很大，我们选择了使用对数反比标度(即 -log(w _0j ) 而不是 w ₀ _j )。当 m = 1.05 、 P _C 约为 0.96，所有重量(除了与 C ₂ 对应的重量)都很小时(记住如果 -log(w) = 30 ，那么 w = e ^-30 )。这样的配置清楚地显示了具有主导成分的非常硬的聚类( C ₂ )。上图后续三个地块继续呈现优势，但是，当 m 增加(而P_CT35】减少)时，优势成分和次优势成分的差异越来越小。这种效果证实了模糊度的增加，对于 m > 1.38 达到最大值。事实上，当 m = 1.5 时，即使 P _C ≈ 0.1 ，所有的权重几乎是相同的，测试样本也不容易被分配到优势聚类。正如我们之前所讨论的，我们现在知道像 k-means 这样的算法可以很容易地找到硬划分，因为平均来说，对应于不同数字的样本彼此之间有很大的不同，并且欧几里德距离足以将它们分配到正确的质心。在这个例子中，我们希望保持适度的模糊性；因此，我们选择了 m = 1.2 (对应P**_C≈0.73):

fc, W, _, _, _, _, pc = cmeans(X.T, c=10, m=1.2, error=1e-6, maxiter=20000, seed=1000)
Mu = fc.reshape((10, 8, 8))

Mu数组包含质心，如下图所示:

Centroids corresponding to m = 1.2 and P_C ≈ 0.73

可以看到，所有不同的数字都被选中，不出所料，第三组(由C₂T5】表示)对应数字 0。现在我们来检查一下X[0]对应的权重(同样， W 是换位的，所以存储在W[:, 0]中):

print(W[:, 0])

输出如下:

[2.68474857e-05 9.14566391e-06 9.99579876e-01 7.56684450e-06
 1.52365944e-05 7.26653414e-06 3.66562441e-05 2.09198951e-05
 2.52320741e-04 4.41638611e-05]

即使作业不是特别辛苦，集群C_T5【2】T6的统治力还是显而易见的。第二个潜在赋值是 C ₈ ，对应数字 9(比例约为 4000)。这样的结果与数字的形状绝对一致，而且考虑到最大重量和第二个重量之间的差异，很明显大多数样本几乎不会被赋值(也就是说，就像在 k 均值中一样)，即使有 P _C ≈ 0.75 。为了检查硬分配(使用权重矩阵上的argmax函数获得)的性能，并考虑到我们知道基本事实，我们可以使用adjusted_rand_score，如下所示:

from sklearn.metrics import adjusted_rand_score

Y_pred = np.argmax(W.T, axis=1)

print(adjusted_rand_score(Y, Y_pred))

上一个片段的输出如下:

0.6574291419247339

这样的值证实了大多数样本已经被成功硬分配。作为补充练习，让我们找出权重标准差最小的样本:

im = np.argmin(np.std(W.T, axis=1))

print(im)
print(Y[im])
print(W[:, im])

输出如下:

414
8
[0.09956437 0.05777962 0.19350572 0.01874303 0.15952518 0.04650815
 0.05909216 0.12910096 0.17526108 0.06091973]

样本*X【414】*代表一个数字(8)，如下图截图所示:

Plot of the sample, X[414], corresponding to the weight vector with the smallest standard deviation

在这种情况下，有三个优势集群: C ₈ 、 C ₄ 、 C ₇ (按降序排列)。不幸的是，它们都不对应数字 8，数字 8 与C₅T15 相关联。不难理解，这样的错误主要是由于数字下半部分的畸形，这产生了更类似于 9 的结果(这样的错误分类也可能发生在人类身上)。然而，低标准偏差和缺乏明显的主导成分应该告诉我们，决策不容易做出，样本具有属于三个主要类别的特征。一个更复杂的监督模型可以很容易地避免这个错误，但是考虑到我们正在执行一个无监督的分析，并且我们仅仅使用基础事实来进行评估，结果并不是那么负面。我建议你用其他 m 值来检验结果，并尝试找出一些潜在的组成规则(即 8 个数字的大部分被柔和地赋给 C _i 和 C _j ，这样我们就可以假设对应的质心编码了部分共有的特征，例如被所有 8 个和 9 个数字共享)。

我们现在可以讨论高斯混合的概念，这是一种非常广泛使用的建模数据集分布的方法，其特征是由低密度区域包围的密集斑点。

高斯混合

高斯混合是最著名的软聚类方法之一，有几十种具体的应用。它可以被认为是 k-means 之父，因为它的工作方式非常相似；但是，与该算法相反，给定样本 x _i ∈ X 和 k 聚类(其被表示为高斯分布)，它提供概率向量*【p(X_I**∈C₁)，...，p(x_I∈C_k)】*。

更一般地说，如果数据集 X 是从数据生成过程 p _数据T5】中采样的，则高斯混合模型基于以下假设:

换句话说，数据生成过程由多元高斯分布的加权和来近似。这种分布的概率密度函数如下:

每个多元高斯的每个分量的影响取决于协方差矩阵的结构。下图显示了二元高斯分布的三种主要可能性(结果可以很容易地扩展到 n 维空间):

Full covariance matrix (left); diagonal covariance (center); circular/spherical covariance (right)

从现在开始，我们将一直考虑全协方差矩阵的情况，这允许实现最大的表达能力。很容易理解，当这样的分布是完全对称的(即协方差矩阵是圆形/球形)时，伪聚类的形状与 k 均值相同(当然，在高斯混合中，聚类没有边界，但在固定数量的标准差后，总是有可能剪切高斯)。相反，当协方差矩阵不是对角的或具有不同的方差时，影响不再对称(例如，在二元的情况下，一个分量可能比另一个分量显示更大的方差)。在这两种情况下，高斯混合允许我们计算实际概率，而不是测量样本之间的距离， x _i ，和均值向量， μ _j (如 k 均值)。下图显示了一个单变量混合的示例:

Example of a univariate Gaussian mixture

在这种情况下，每个样本在每个高斯下总是有一个非空概率，其影响由其均值和协方差矩阵决定。例如，对应于 x 位置 2.5 的点既可以属于中心高斯，也可以属于右侧高斯(而左侧高斯的影响最小)。正如本章开头所解释的，任何软聚类算法都可以通过选择影响最大的组件(argmax)而转变为硬聚类算法。

您将立即理解，在这种特定情况下，通过对角协方差矩阵，argmax提供了一条附加信息(被 k 均值完全丢弃)，可以在进一步的处理步骤中使用(也就是说，推荐器应用程序可以提取所有聚类的主要特征，并根据相对概率重新加权它们)。

高斯混合模型的电磁算法

完整的算法(在掌握机器学习算法中有充分描述，作者:Bonaccorso G .，Packt Publishing，2018)比 k-means 稍微复杂一点，需要更深层次的数学知识。由于这本书的范围更实用，我们只讨论主要步骤，没有提供正式证据。

让我们首先考虑一个数据集， X ，包含 n 个样本:

给定 k 分布，我们需要找到权重，*w_jT5】，以及每个高斯 *(μ _j ，σ_j)的参数，条件如下:

最后一个条件对于保持与概率定律的一致性是必要的。如果我们将所有参数分组为一个集合， θ _j = (w _j 、μ _j 、σ_j)，我们可以定义样本x_IT13】在jT16】th 高斯下的概率，如下:

类似的，我们可以引入一个伯努利分布，z_IT3】j= p(j | x_I，θ_j)∞B(p)，这是jT12】th 高斯生成样本x_IT17】的概率。换句话说，给定一个样本， x _i ，z _ij 将等于 1，概率 p(j|x _i ，θ _j ) ，否则为 0。

此时，我们可以计算整个数据集的联合对数似然，如下所示:

在前面的公式中，我们利用了指数指标表示法，它依赖于这样一个事实，即z_ijT3 只能是 0 或 1。因此，当 z _ij = 0 时，意味着样品 x _i 还没有由 j ^th 高斯生成，产品中对应的项变为 1(即 x ⁰ = 1 )。反之，当 z _ij = 1 时，该项等于x_IT27】和jT30】th 高斯的联合概率。因此，联合对数似然是模型生成整个数据集的联合概率，假设每个 x _i ∈ X 是独立同分布 ( IID )。要解决的问题是最大似然估计 ( MLE )，或者换句话说，寻找最大化L(θ；x，Z) 。但是变量 z _ij 并没有被观察到(或者潜伏)，所以不可能直接最大化可能性，因为我们不知道它们的值。**

解决这一问题的最有效方法是采用 em 算法(由 Dempster A. P，Laird N. M .和 Rubin D. B .提出，通过 EM 算法从不完整数据中获得最大似然，皇家统计学会杂志，系列 B. 39 (1)，1977)。完整的解释超出了本书的范围，但我们想提供主要步骤。首先要做的是使用概率链规则，以便将前面的表达式转换为条件概率的总和(这可以很容易地管理):

这两种可能性现在很简单。术语 p(x _i |j，θ _j ) 是jT12】th 高斯下的x_IT9】的概率，而 p(j|θ _j ) 只是jT20】th 高斯下的概率，相当于权重，为了消除潜在的变量，EM 算法以迭代的方式进行，由两个步骤组成。第一个(称为期望步骤，或 E 步骤)是计算没有潜在变量的可能性的代理。如果我们将整个参数集表示为 θ ，并且在迭代 t 时计算的同一个参数集表示为θ_tT37】，我们可以计算出以下函数:

Q( θ|θ _t ) 是关于变量 z _ij 的联合对数似然的期望值，并被调节到数据集 X 和迭代时的参数集 t 。这种操作的效果是消除潜在的变量(这些变量被求和或积分出来)，并产生实际对数似然的近似值。不难想象，第二步(称为最大化-步，或M-步)的目标是最大化 Q(θ|θ _t ) ，生成一个新的参数集， θ _t+1 。重复该过程，直到参数变得稳定，并且有可能证明最终参数集对应于最大似然估计。跳过所有中间步骤，假设最佳参数设置为 θ _f ，最终结果如下:

为了清楚起见，概率 p(j|x _i ，θ _f ) 可以通过使用贝叶斯定理来计算:

比例性可以通过归一化所有项来消除，使得它们的和等于 1(满足概率分布的要求)。

现在，让我们考虑一个使用 scikit-learn 的实际例子。由于目标纯粹是说教，我们使用了一个易于可视化的二维数据集:

from sklearn.datasets import make_blobs

nb_samples = 300
nb_centers = 2

X, Y = make_blobs(n_samples=nb_samples, n_features=2, center_box=[-1, 1], centers=nb_centers, cluster_std=[1.0, 0.6], random_state=1000)

数据集是从两个具有不同标准差(1.0 和 0.6)的高斯分布中采样生成的，如下图所示:

Dataset for the Gaussian mixture example

我们的目标是使用高斯混合模型和 k 均值，并比较最终结果。正如我们预期的两个组成部分，数据生成过程的近似如下:

我们现在可以用n_components=2训练一个GaussianMixture实例。默认的协方差类型是完整的，但是可以通过设置covariance_type参数来更改该选项。允许的值是full、diag、spherical和tied(这迫使算法对所有高斯函数使用共享的单个协方差矩阵):

from sklearn.mixture import GaussianMixture

gm = GaussianMixture(n_components=2, random_state=1000)
gm.fit(X)
Y_pred = gm.fit_predict(X)

print('Means: \n{}'.format(gm.means_))
print('Covariance matrices: \n{}'.format(gm.covariances_))
print('Weights: \n{}'.format(gm.weights_))

上一个片段的输出如下:

Means: 
[[-0.02171304 -1.03295837]
 [ 0.97121896 -0.01679101]]

Covariance matrices: 
[[[ 0.86794212 -0.18290731]
  [-0.18290731  1.06858097]]

 [[ 0.44075382  0.02378036]
  [ 0.02378036  0.37802115]]]

Weights: 
[0.39683899 0.60316101]

因此，最大似然估计产生了两个分量，其中稍占优势的一个分量(即 w ₂ = 0.6 )。为了知道高斯轴的方向，我们需要计算协方差矩阵的归一化特征向量(这个概念将在第 7 章、降维和分量分析中充分解释):

import numpy as np

c1 = gm.covariances_[0]
c2 = gm.covariances_[1]

w1, v1 = np.linalg.eigh(c1)
w2, v2 = np.linalg.eigh(c2)

nv1 = v1 / np.linalg.norm(v1)
nv2 = v2 / np.linalg.norm(v2)

print('Eigenvalues 1: \n{}'.format(w1))
print('Eigenvectors 1: \n{}'.format(nv1))

print('Eigenvalues 2: \n{}'.format(w2))
print('Eigenvectors 2: \n{}'.format(nv2))

输出如下:

Eigenvalues 1: 
[0.75964929 1.17687379]
Eigenvectors 1: 
[[-0.608459   -0.36024664]
 [-0.36024664  0.608459  ]]

Eigenvalues 2: 
[0.37002567 0.4487493 ]
Eigenvectors 2: 
[[ 0.22534853 -0.6702373 ]
 [-0.6702373  -0.22534853]]

在两个二元高斯(一旦被截断并从顶部观察，可以想象成椭圆)中，主要成分是第二个(即第二列，对应于最大的特征值)。椭圆的偏心率由特征值之间的比值决定。如果这样的比例等于 1，形状就是圆，高斯是完全对称的；否则，它们会沿着一个轴拉伸。主要部件与 x 轴之间的角度(度)如下:

import numpy as np

a1 = np.arccos(np.dot(nv1[:, 1], [1.0, 0.0]) / np.linalg.norm(nv1[:, 1])) * 180.0 / np.pi
a2 = np.arccos(np.dot(nv2[:, 1], [1.0, 0.0]) / np.linalg.norm(nv2[:, 1])) * 180.0 / np.pi

之前的公式是基于主成分 v ₁ 与 x -versor、e₀T9】(即【1，0】)之间的点积:

在显示最终结果之前，使用 k 均值对数据集进行聚类将会很有帮助:

from sklearn.cluster import KMeans

km = KMeans(n_clusters=2, random_state=1000)
km.fit(X)
Y_pred_km = km.predict(X)

聚类结果如下图所示:

Gaussian mixture result (left) with the shapes of three horizontal sections; k-means result (right)

正如预期的那样，两种算法产生非常相似的结果，主要的差异是由于高斯分布的非对称性。特别地，对应于数据集左下部分的伪聚类在两个方向上具有较大的方差，并且对应的高斯是主导的。为了检查混合物的行为，让我们计算三个样本点的概率( (0，-2)； (1，-1)—一个临界样本；和 (1，0) ，使用predict_proba()方法:

print('P([0, -2]=G1) = {:.3f} and P([0, -2]=G2) = {:.3f}'.format(*list(gm.predict_proba([[0.0, -2.0]]).squeeze())))
print('P([1, -1]=G1) = {:.3f} and P([1, -1]=G2) = {:.3f}'.format(*list(gm.predict_proba([[1.0, -1.0]]).squeeze())))
print('P([1, 0]=G1) = {:.3f} and P([1, 0]=G2) = {:.3f}'.format(*list(gm.predict_proba([[1.0, 0.0]]).squeeze())))

上一个块的输出如下:

P([0, -2]=G1) = 0.987 and P([0, -2]=G2) = 0.013
P([1, -1]=G1) = 0.354 and P([1, -1]=G2) = 0.646
P([1, 0]=G1) = 0.068 and P([1, 0]=G2) = 0.932

我邀请读者通过使用其他协方差类型来重复这个例子，然后用 k 均值来比较所有的硬赋值。

用 AIC 和 BIC 评估高斯混合模型的性能

由于高斯混合模型是一个概率模型，找到最佳的组件数量需要一种不同于前面几章分析的方法。最广泛使用的技术之一是阿卡克信息标准 ( AIC )，其基于信息论(首次在阿卡克首次提出，统计模型识别的新观点，IEEE 自动控制交易， 19 (6))。如果概率模型具有n_pT9】参数(即必须学习的单个值)并达到最大负对数似然， L _opt ，则 AIC 定义如下:

这种方法有两个重要的含义。第一个是关于价值本身；AIC 越小，分数越高。事实上，考虑到奥卡姆剃刀原理，模型的目标是以最小的参数数量实现最优似然。第二个含义与信息论严格相关(我们不讨论细节，这些细节在数学上很重要)，特别是数据生成过程和一般概率模型之间的信息损失。有可能证明 AIC 的渐近最小化(即当样本数趋于无穷大时)等价于信息损失的最小化。考虑基于不同数量的分量的几个高斯混合( n _p 是所有权重、均值和协方差参数的总和)，具有最小 AIC 的配置对应于以最高精度再现数据生成过程的模型。AIC 的主要局限在于小数据集。在这种情况下，AIC 倾向于在大量参数下达到最小值，这与奥卡姆剃刀原理形成对比。然而，在大多数现实生活中，AIC 提供了一个有用的相对度量，可以帮助数据科学家排除许多配置，只分析最有希望的配置。

当需要强制参数数量保持较低时，可以采用贝叶斯信息准则 ( BIC ，定义如下:

上式中， n 为样本数(例如 n = 1000 且采用自然对数，罚值约为 6.9)；因此，BIC 几乎等同于 AIC，对参数的数量有更强的惩罚。然而，即使 BIC 倾向于选择较小的模型，结果通常也不如 AIC 可靠。BIC 的主要优点是当 n → ∞ 时，数据生成过程 p _数据T9】与模型p_mT13】之间的 kulback-Leibler 散度(BIC 最小)趋于 0:**

由于当两个分布相同时，库尔巴克-莱布勒散度为零，所以前一个条件意味着 BIC 趋向于渐近地选择精确再现数据生成过程的模型。

现在，让我们考虑前面的例子，检查 AIC 和 BIC 的不同数量的组件。Scikit-learn 将这些措施作为GaussianMixture课程的方法(aic()和bic())纳入其中。此外，我们还想计算每个模型所达到的最终对数似然。这可以通过将通过score()方法获得的值(每样本对数似然平均值乘以样本数)相乘来实现，如下所示:

from sklearn.mixture import GaussianMixture

n_max_components = 20

aics = []
bics = []
log_likelihoods = []

for n in range(1, n_max_components + 1):
 gm = GaussianMixture(n_components=n, random_state=1000)
 gm.fit(X)
 aics.append(gm.aic(X))
 bics.append(gm.bic(X))
 log_likelihoods.append(gm.score(X) * nb_samples)

结果图如下图所示:

AICs, BICs, and log-likelihoods for Gaussian mixtures with the number of components in the range (1, 20)

在这种情况下，我们知道数据集是由两个高斯分布生成的，但是让我们假设我们没有这条信息。AIC 和 BIC 的(当地)最低气温都是 n _c = 2 。然而，当 BIC 不断变得越来越大时，AIC 的伪全球最小值为 n _c = 18 。因此，如果我们信任 AIC，我们应该选择 18 个分量，这相当于用许多高斯子超分割数据集，具有小的方差。另一方面， n _c = 2 与 n _c = 18 相比其他数值的差异并不是很大，所以我们也可以更倾向于前者的配置，考虑到要简单得多。BIC 证实了这一选择。事实上，即使也有对应于 n _c = 18 的局部最小值，其值也比 n _c = 2 的 BIC 值大得多。正如我们之前解释的，这种行为是由于 BIC 强加的样本量的额外惩罚。由于每个二元高斯需要一个权重变量，两个均值变量，四个协方差矩阵变量，对于 n _c = 2 ，我们得到n_p= 2(1+2+4)= 14，对于 n _c = 18 ，我们得到 n _p = 18(1 由于有 300 个样本，BIC 被罚 log(300) ≈ 5.7 ，导致 BIC 增加约 350 个。当 n _c 变大时，对数似然性增加(因为在极端情况下，每个点可以被认为是由具有零方差的单个高斯生成的，相当于狄拉克δ)，参数的数量在模型选择过程中起着主要作用。

没有任何额外的惩罚，更大的模型很可能被选为最佳选择，但是在聚类过程中，我们也需要执行最大分离原则。这种情况部分与较少的组件有关，因此，BIC 应成为最佳方法。总的来说，我建议比较这两个标准，试图找到对应于 AIC 和 BIC 之间最大一致的n_cT3。此外，还应该考虑基本的背景知识，因为许多数据生成过程都有明确定义的行为，并且可以通过排除所有不现实的值来限制潜在组件的范围。我请读者用 n _c = 18 重复前面的例子，画出所有高斯人，比较某些特定点的概率。

基于贝叶斯高斯混合的组件选择

贝叶斯高斯混合模型是基于变分框架的标准高斯混合模型的扩展。这个题目相当高深，需要透彻的数学描述，超出了本书的范围(你可以在 Nasios N .和 Bors A. G .，高斯混合模型的变分学习，IEEE Transactions On Systems，Man，and 控制论， 36/ 4，08/2006 中找到)。然而，在我们讨论主要属性之前，理解主要概念和差异将是有帮助的。假设我们有一个数据集， X ，以及一个用向量 θ 参数化的概率模型。前面几节你看到了概率p(X |**【θ】)就是可能性 L(θ|X) ，它的最大化导致了一个生成概率最大的 X 的模型。但是，我们不对参数施加任何约束，它们的最终值完全取决于 X 。如果我们引入贝叶斯定理，我们会得到以下结果:

左边是给定数据集的参数后验概率，我们知道它与参数的似然性乘以先验概率成正比。在一个标准的最大似然估计中，我们只使用 p(X|θ) ，但是，当然，我们也可以包含一段关于 θ (根据概率分布)的先验知识，并最大化 p(θ|X) 或比例代理函数。然而，一般来说， p(θ|X) 是难以处理的，之前的 p(θ) 往往很难定义，因为没有足够的关于高度可能区域的知识。为此，最好将参数建模为用 η (所有特定参数的集合，如平均值、系数等)参数化的概率分布，并引入变分后验、q(θ| X；η )近似真实分布。

这样的工具是一种称为变分贝叶斯推理的技术的关键要素(你可以在前面提到的论文中找到进一步的细节)，它允许我们轻松找到最佳参数，而不需要与实际的 p(θ|X) 一起工作。特别是，在高斯混合模型中，有三组不同的参数，每组参数都用适当的分布建模。在这种情况下，我们倾向于不讨论这些选择的细节，但是理解其基本原理是有用的。

在贝叶斯框架中，给定一个似然性， p(X|θ) ，一个概率密度函数， p(θ) ，属于同一个后验族， p(θ|X) ，被称为共轭先验。在这种情况下，程序显然被简化了，因为可能性的影响仅限于修改先前的参数。为此，由于似然性是正态的，为了对均值建模，我们可以采用正态分布(它是相对于均值的共轭先验)，对于协方差矩阵，我们可以采用威沙特分布(它是相对于协方差矩阵的逆的共轭先验)。在本讨论中，没有必要熟悉所有这些分布(除了正态分布)，但记住它们是共轭先验是有帮助的，因此，给定参数的初始猜测，可能性的作用是调整它们，以便在给定数据集的情况下最大化它们的联合概率。

由于混合的权重被归一化，因此它们的和必须总是等于 1，并且我们希望只自动选择大量分量的子集，因此我们可以使用狄利克雷分布，它具有稀疏的有用特性。换句话说，给定一组权重 w ₁ ，w₂T5，...，以及 w _n ，狄利克雷分布倾向于保持大部分权重的概率相当低，而非空权重的较小子组决定了主要贡献。狄利克雷过程提供了另一种选择，狄利克雷过程是一种产生概率分布的特殊随机过程。在这两种情况下，目标都是调整单个参数(称为重量浓度参数，该参数增加或减少具有稀疏分布(或简单地说，狄利克雷分布的稀疏性)的概率。

Scikit-learn 实现了贝叶斯高斯混合(通过BayesianGaussianMixture类)，它可以基于狄利克雷过程和分布。在本例中，我们将保留默认值(process)并检查不同浓度值的行为(weight_concentration_prior参数)。也可以调整高斯的平均值作为均值，调整威沙特的自由度作为逆协方差。然而，在没有任何特定的先验知识的情况下，很难设置这些值(我们假设我们不知道均值可能位于何处或协方差矩阵的结构)，因此，最好保留从问题结构中导出的值。因此，平均值(高斯)将等于 X 的平均值(可以用mean_precision_prior参数控制位移；值 < 1.0 倾向于将单个均值向 X 的均值移动，而较大的值会增加位移)，自由度(Wishart)的数量被设置为等于特征的数量(维度为 X )。在许多情况下，这些参数是由学习过程自动调整的，没有必要改变它们的初始值。相反，weight_concentration_prior可以调整，以增加或减少活性成分的数量(即，其重量不接近零或比其他成分低得多)。

在这个例子中，我们将生成 500 个二维样本，使用 5 个部分重叠的高斯分布(特别是，其中 3 个共享一个非常大的重叠区域):

from sklearn.datasets import make_blobs

nb_samples = 500
nb_centers = 5

X, Y = make_blobs(n_samples=nb_samples, n_features=2, center_box=[-5, 5], 
                  centers=nb_centers, random_state=1000)

先从大重量浓度参数(1000)开始，最大成分数等于5。在这种情况下，我们期望找到大量(可能是5)的活性成分，因为狄利克雷过程不能实现高度稀疏:

from sklearn.mixture import BayesianGaussianMixture

gm = BayesianGaussianMixture(n_components=5, weight_concentration_prior=1000, 
                             max_iter=10000, random_state=1000)
gm.fit(X)

print('Weights: {}'.format(gm.weights_))

上一个片段的输出如下:

Weights: [0.19483693 0.20173229 0.19828598 0.19711226 0.20803253]

正如所料，所有组件的重量大致相同。为了得到进一步的确认，我们可以检查(通过argmax功能)有多少样本几乎没有分配给每个样本，如下所示:

Y_pred = gm.fit_predict(X)

print((Y_pred == 0).sum())
print((Y_pred == 1).sum())
print((Y_pred == 2).sum())
print((Y_pred == 3).sum())
print((Y_pred == 4).sum())

输出如下:

因此，平均来说，所有高斯人产生的点数是相同的。最终配置如下图所示:

Final configuration, with five active components

这种模式一般可以接受；然而，让我们假设我们知道潜在原因(即生成高斯分布)的数量可能是 4，而不是 5。我们可以尝试的第一件事是保持原有的最大组分数，降低重量浓度参数(即 0.1)。如果近似能够使用较少的高斯分布成功生成 X ，我们应该找到一个零权重:

gm = BayesianGaussianMixture(n_components=5, weight_concentration_prior=0.1, 
                             max_iter=10000, random_state=1000)

gm.fit(X)

print('Weights: {}'.format(gm.weights_))

现在输出如下:

Weights: [3.07496936e-01 2.02264778e-01 2.94642240e-01 1.95417680e-01 1.78366038e-04]

可以看到，第五个高斯的权重比其他高斯小得多，可以完全丢弃(我邀请你检查是否有几个样本几乎没有分配给它)。新配置包含四个活动组件，如下图所示:

Final configuration, with four active components

正如可以看到的，模型已经执行了组件数量的自动选择，并且它已经将较大的右斑点分成两部分，这两部分几乎是正交的。即使使用大量初始组件(例如，10 个)训练模型，该结果也保持不变。作为练习，我建议用其他值重复该示例，检查权重之间的差异。贝叶斯高斯混合模型非常强大，因为它们能够避免过度拟合。事实上，虽然标准高斯混合模型将通过减少它们的协方差来使用所有的分量，如果必要的话(以便覆盖密集的区域)，但是这些模型利用了狄利克雷过程/分布的属性，以避免激活太多的分量。例如，通过检查模型可实现的最小组件数量，可以很好地洞察潜在的数据生成过程。在没有任何其他先验知识的情况下，这样的值是最终配置的一个很好的候选值，因为较少数量的组件也会产生较低的最终可能性。当然，也可以将 AIC/BIC 与这种方法结合使用，以获得另一种形式的确认。然而，与标准高斯混合模型的主要区别在于，它可能包含来自专家的先验信息(例如，均值和协方差方面的原因结构)。为此，我邀请您通过更改mean_precision_prior的值来重复该示例。例如，可以将mean_prior参数设置为不同于 X 平均值的值，并调整mean_precision_prior，从而迫使模型基于一些先验知识实现不同的分割(即一个区域中的所有样本都应该由特定的分量生成)。

生成高斯混合

高斯混合模型主要是生成模型。这意味着训练过程的目标是优化参数，以便最大化模型生成数据集的可能性。如果假设是正确的，并且 X 已经从特定的数据生成过程中采样，则最终的近似值必须能够生成所有其他潜在的样本。换句话说，我们假设 x _i ∈ X 为 IDD，X_I∨p_数据；因此，当已经找到最优近似 p ≈ p _数据 时，所有在 p 下概率高的样本x_jT19】也很可能是由 p _数据 生成的。

在这个例子中，我们希望在半监督场景中使用高斯混合模型。这意味着我们有一个包含标记和未标记样本的数据集，我们希望利用标记样本作为基础事实，并找到可以生成整个数据集的最佳混合。当标注非常大的数据集非常困难和昂贵时，这种情况非常常见。为了克服这个问题，可以标记一个均匀采样的子集，并训练一个生成模型，该模型能够生成具有最大可能可能性的剩余样本。

我们将使用更新后的权重、均值和协方差矩阵公式，这些公式在主要部分通过一个简单的过程进行了讨论，如下所示:

所有被标记的样本都被认为是基本事实；所以，如果有 k 个类，我们还需要定义 k 个组件，并将每个类分配给其中一个。因此，如果x_IT5】是标有 y _i = {1，2，...，k} ，对应的概率向量将为 p(x _i ) = (0，0，..., 1, 0, ...，0) ，其中 1 对应于与 y _i 类相关联的高斯。换句话说，我们信任已标记的样本，并强制单个高斯生成具有相同标签的子集。
所有未标记的样本都以标准方式处理，概率向量是通过将权重乘以每个高斯下的概率来确定的。

让我们首先生成一个包含 500 个二维样本的数据集(100已标记，其余的未标记)，具有真实的标记，0和1，以及等于-1的未标记标记:

from sklearn.datasets import make_blobs

nb_samples = 500
nb_unlabeled = 400

X, Y = make_blobs(n_samples=nb_samples, n_features=2, centers=2, cluster_std=1.5, random_state=100)

unlabeled_idx = np.random.choice(np.arange(0, nb_samples, 1), replace=False, size=nb_unlabeled)
Y[unlabeled_idx] = -1

此时，我们可以初始化高斯参数(权重被选择为相等，并且协方差矩阵必须是半正定的。如果读者不熟悉这个概念，我们可以说对称方阵a∈ℜ^n×nT3】是正半定的，如果:

而且所有的特征值都是非负的，特征向量生成一个正交基(这个概念在第 7 章、降维和成分分析讲 PCA 的时候会很有帮助)。

如果协方差矩阵是随机选择的，为了成为半正定的，需要将它们中的每一个乘以它的转置):

import numpy as np

m1 = np.array([-2.0, -2.5])
c1 = np.array([[1.0, 1.0],
               [1.0, 2.0]])
q1 = 0.5

m2 = np.array([1.0, 3.0])
c2 = np.array([[2.0, -1.0],
               [-1.0, 3.5]])
q2 = 0.5

数据集和初始高斯分布如下图所示:

Dataset (the unlabeled samples are marked with an x) and initial configuration

现在，我们可以按照之前定义的规则执行几次迭代(在我们的例子中是 10 次)(当然，也可以检查参数的稳定性，以便停止迭代)。使用 SciPy multivariate_normal类计算每个高斯分布下的概率:

from scipy.stats import multivariate_normal

nb_iterations = 10

for i in range(nb_iterations):
    Pij = np.zeros((nb_samples, 2))

    for i in range(nb_samples):

        if Y[i] == -1:
            p1 = multivariate_normal.pdf(X[i], m1, c1, allow_singular=True) * q1
            p2 = multivariate_normal.pdf(X[i], m2, c2, allow_singular=True) * q2
            Pij[i] = [p1, p2] / (p1 + p2)
        else:
            Pij[i, :] = [1.0, 0.0] if Y[i] == 0 else [0.0, 1.0]

    n = np.sum(Pij, axis=0)
    m = np.sum(np.dot(Pij.T, X), axis=0)

    m1 = np.dot(Pij[:, 0], X) / n[0]
    m2 = np.dot(Pij[:, 1], X) / n[1]

    q1 = n[0] / float(nb_samples)
    q2 = n[1] / float(nb_samples)

    c1 = np.zeros((2, 2))
    c2 = np.zeros((2, 2))

    for t in range(nb_samples):
        c1 += Pij[t, 0] * np.outer(X[t] - m1, X[t] - m1)
        c2 += Pij[t, 1] * np.outer(X[t] - m2, X[t] - m2)

    c1 /= n[0]
    c2 /= n[1]

过程结束时的高斯混合参数如下:

print('Gaussian 1:')
print(q1)
print(m1)
print(c1)

print('\nGaussian 2:')
print(q2)
print(m2)
print(c2)

上一个片段的输出如下:

Gaussian 1:
0.4995415573662937
[ 0.93814626 -4.4946583 ]
[[ 2.53042319 -0.10952365]
 [-0.10952365  2.26275963]]

Gaussian 2:
0.5004584426337063
[-1.52501526  6.7917029 ]
[[ 2.46061144 -0.08267972]
 [-0.08267972  2.54805208]]

不出所料，由于数据集的对称性，权重几乎保持不变，同时更新了均值和协方差矩阵，以最大化可能性。最后的图显示在下面的截图中:

Final configuration, after 10 iterations

可以看到，这两个 Gaussians 都已经成功优化，能够从几个扮演可信向导角色的标注样本开始生成整个数据集。这样的方法非常强大，因为它允许我们在模型中包含一些先验知识，而无需任何修改。然而，由于标记样本具有等于 1 的固定概率，该方法在异常值方面不是非常稳健。如果数据生成过程没有生成样本，或者样本受到噪声的影响，模型可能会导致高斯错误。然而，这种情况通常应该被忽略，因为任何先验知识，当包括在估计中时，总是必须被预先评估，以便检查它是否可靠。这样的步骤是必要的，以避免迫使模型只学习原始数据生成过程的一部分的风险。相反，当标记样本真正代表底层过程时，它们的包含减少了误差并加速了收敛。我邀请读者在介绍一些有噪声的点(例如，(-20，-10))后重复这个例子，并比较几个未标记测试样本的概率。

摘要

在本章中，我们介绍了一些最常见的软聚类方法，重点介绍了它们的属性和特点。模糊 c 均值是经典 k 均值算法的扩展，基于模糊集的概念。集群不被认为是互斥的分区，而是可以与其他一些集群重叠的灵活集合。所有的样本总是被分配给所有的聚类，但是权重向量决定了每个聚类的隶属度。连续的簇可以定义部分重叠的属性；因此，对于两个或多个聚类，给定样本可以具有非空权重。大小决定了它属于每个片段的多少。

高斯混合是一个生成过程，它基于一个假设，即可以用高斯分布的加权和来近似真实的数据生成过程。给定预定义数量的组件，为了最大化可能性，对模型进行训练。我们讨论了如何使用 AIC 和 BIC 作为性能度量，以便找到高斯分布的最佳数量。我们还简要介绍了贝叶斯高斯混合的概念，并研究了包含先验知识如何有助于自动选择一小部分活动组件。在最后一部分，我们讨论了半监督高斯混合的概念，展示了如何使用一些标记样本作为指导，以优化具有大量未标记点的训练过程。

在下一章中，我们将讨论核密度估计的概念及其在异常检测领域的应用。

问题

软聚类和硬聚类的主要区别是什么？
模糊 c 均值可以很容易地处理非凸聚类。这个说法对吗？
高斯混合模型的主要假设是什么？
假设两个模型达到相同的最优对数似然；然而，第一个的 AIC 是第二个的两倍。这是什么意思？
考虑到前面的问题，我们更喜欢哪种模式？
为什么我们要使用狄利克雷分布作为贝叶斯高斯混合模型权重的先验？
假设我们有一个包含 1000 个标记样本的数据集，这些样本的值已经过专家认证。我们从相同的来源收集了 5000 个样本，但我们不想为额外的标签付费。我们可以做些什么来将它们纳入我们的模型？

进一步阅读

理论神经科学，达扬·p .，麻省理工学院出版社，2005
通过 EM 算法不完全数据的最大似然，皇家统计学会杂志，登普斯特 A. P .，莱尔德 N. M .，鲁宾 D. B. ，系列 B. 39 (1)，1977
统计模型识别的新面貌、阿凯克 H. 、 IEEE 自动控制交易，19 (6)
高斯混合模型的变分学习，纳西欧和博思公司， IEEE 系统、人和控制论学报，36/ 4，08/2006
Belohlavek R .，Klir G. J. (编著)，概念与模糊逻辑，麻省理工学院出版社，2011
查佩尔·奥、斯科尔科夫·b、齐恩·a .(编著)半监督学习**麻省理工学院出版社，2010
掌握机器学习算法，博纳科索格，帕克特出版，2018
机器学习算法，第二版，博纳科索格，帕克特出版，2018

六、异常检测

在本章中，我们将讨论无监督学习的一个实际应用。我们的目标是训练模型，这些模型要么能够重现特定数据生成过程的概率密度函数，要么能够识别给定的新样本是内联样本还是外联样本。一般来说，我们可以说我们要追求的具体目标是发现异常，这些异常往往是模型下可能性极小的样本(即给定概率分布 p(x) < < λ 其中 λ 是预定义的阈值)或者离主分布质心相当远的样本。

特别是，本章将包括以下主题:

概率密度函数及其基本性质简介
直方图及其局限性
籽粒密度估算 ( KDE )
带宽选择标准
异常检测的单变量示例
使用 KDD 杯 99 数据集的 HTTP 攻击异常检测示例
一类支持向量机
利用隔离森林进行异常检测

技术要求

本章中的代码要求:

Python 3.5+(强烈推荐 Anaconda 发行版(www.anaconda.com/distributio…)
库:
- SciPy 0.19+
- NumPy 1.10+
- 学习 0.20+
- 熊猫 0.22+
- Matplotlib 2.0+
- seaborn 0.9+

这些例子可以在 GitHub 资源库上获得，网址为 https://GitHub . com/packktpublishing/HandsOn-Unsupervised-Learning-with-Python/tree/master/chapter 06。

概率密度函数

在前面的所有章节中，我们一直假设我们的数据集是从隐式数据生成过程 p _data 中提取的，并且所有的算法都假设 x _i ∈ X 为独立同分布 ( IID )并且均匀采样。我们假设 X 以足够的精度表示 p _数据 ，这样算法就可以在有限的初始知识下学会泛化。相反，在本章中，我们感兴趣的是在没有任何特定限制的情况下直接建模 p _数据 (例如，高斯混合模型通过对分布的结构施加约束来实现这一目标)。在讨论一些非常强大的方法之前，简要回顾一下在可测子集 X ℜn 上定义的一般连续概率密度函数 p(x) 的性质是有帮助的(为了避免混淆，我们将使用 p(x) 表示密度函数，使用 P(x) 表示实际概率):

例如，单变量高斯分布完全由均值 μ 和方差σ²T5 表征:

因此，x∑(a，b) 的概率如下:

即使一个事件在连续空间(例如，高斯)中的绝对概率为零(因为积分具有相同的极值)，概率密度函数也提供了一个非常有用的度量，可以了解一个样本比另一个样本更有可能发生多少。例如:考虑高斯分布 N(0，1) ，密度 p(1) = 0.4 ，而对于 x = 2 ，密度降低到约 0.05 。意思是 1 比 2 大 0.4 / 0.05 = 8 倍。同样的，我们可以设置一个接受阈值 α ，将p(x_I)<**α的所有样本x_IT19】定义为异常(例如，在我们的例子中， α = 0.01 )。这个选择是异常检测过程中至关重要的一步，正如我们将要讨论的，它还必须包括潜在的异常值，然而，这些异常值仍然是常规样本。

在许多情况下，特征向量是用多维随机变量建模的。例如:数据集 X ℜ ³ 可以用联合概率密度函数 p(x，y，z) 来表示。一般情况下，实际概率需要三重积分:

很容易理解，任何使用这种联合概率的算法都会受到复杂性的负面影响。通过假设单个分量的统计独立性，可以获得很强的简化:

不熟悉这个概念的读者可以在考试前想象一群学生。用随机变量建模的特征是学习小时数( x )和已完成课程数( y )，我们想找出给定这些因素的成功概率 p(成功|x，y) (这样的例子是基于条件概率的，但主要概念总是一样的)。我们可以假设，一个学完所有课程的学生，在家需要少学习；然而，这样的选择意味着这两个因素之间的依赖性(和相关性)，这不能再单独评估了。相反，我们可以通过假设没有任何相关性来简化程序，并根据已完成课程的数量和作业的小时数来计算边际成功概率。重要的是要记住，特征之间的独立性不同于随后从分布中提取的样本的独立性。当我们说数据集由 IID 样本组成时，我们指的是概率p(x_I| x_I_-1，x _i _-2 ，...，p ₁ ) =每个样品的 p(x _i ) 。换句话说，我们假设样本之间没有相关性。这样的条件更容易实现，因为它通常足以混洗数据集以消除任何残留的相关性。相反，特征之间的相关性是数据生成过程的一个特殊属性，无法消除。因此，在某些情况下，我们假设独立性，因为我们知道它的影响可以忽略不计，最终结果不会受到严重影响，而在其他情况下，我们将基于整个多维特征向量来训练模型。我们现在可以定义异常的概念，它将在剩下的部分中使用。

异常作为异常值或新奇值

本章的主题是在没有任何监督的情况下自动检测异常。由于模型不是基于标记样本提供的反馈，我们只能依靠整个数据集的属性来找出相似之处并突出不同之处。特别是，我们从一个非常简单但有效的假设出发:常见事件为正常，而不太可能发生的事件一般被视为异常。当然，这个定义意味着我们正在监控的过程正在正常运行，大多数结果被认为是有效的。例如:一个硅加工厂必须把一个晶片切成相等的块。我们知道它们每一个都是 0.2 × 0.2 英寸(约 0.5 × 0.5 厘米)，每一侧的标准偏差为 0.001 英寸。经过 1，000，000 个处理步骤后，已经确定了该度量。我们是否有权将 0.25 × 0.25 英寸芯片视为异常？当然，我们是。实际上，我们假设每条边的长度都建模为高斯分布(非常合理的选择)*μ= 0.2**σ= 0.001；*三个标准差后，概率下降到几乎为零。因此，例如: P(边> 0.23) ≈ 0 这样的尺寸的芯片必须明确视为异常。

显然，这是一个极其简单的例子，不需要任何模型。然而，在现实生活中，密度的结构可能非常复杂，几个高概率区域被低概率区域包围。这就是为什么必须采用更通用的方法来建模整个样本空间。

当然，异常的语义不能标准化，它总是取决于正在分析的具体问题。因此，定义异常概念的一个常见方法是区分异常值和新奇值。前者是包含在数据集中的样本，即使它们与其他样本之间的距离大于平均值。因此，异常值检测过程旨在找出这样的奇怪的样本(例如:考虑到前面的例子，0.25 × 0.25 英寸的芯片如果包含在数据集中显然是异常值)。相反，新颖性检测的目标略有不同，因为在这种情况下，我们假设使用只包含正常样本的数据集；因此，给定一个新的，我们有兴趣理解我们是否可以认为它是从原始数据生成过程中提取的或者是一个异常值(例如:一个新手技术人员问我们这个问题:0.25 × 0.25 英寸的芯片是一个异常值吗？如果我们已经收集了正常芯片的数据集，我们可以使用我们的模型来回答这个问题。

描述这种情况的另一种方式是将样本视为一系列可能受可变噪声影响的值: y(t) = x(t) + n(t) 。当*| | n(t)| |<|<| | x(t)|时，样品可归类为洁净* : y(t) ≈ x(t) 。反之，当*| | n(t)|∑| | x(t)| |*(甚至更大)时，则为异常值，不能代表真实的底层过程 p _数据 。由于噪声的平均幅度通常比信号小得多，*P(| | n(t)|∞| | x(t)| |)*的概率接近于零。因此，我们可以将异常想象成受异常外部噪声影响的正常样本。管理异常和噪声样本之间的真正主要区别通常在于检测真正异常和相应标记样本的能力。事实上，虽然有噪声的信号肯定是被破坏的，并且目标是最小化噪声的影响，但是异常可以被人类非常频繁地识别并正确标记。然而，正如已经讨论过的，在本章中，我们有兴趣找出不依赖于现有标签的发现方法。此外，为了避免混淆，我们总是引用异常，每次都定义数据集的内容(只有内联或内联和外联)和我们分析的目标。在下一节中，我们将简要讨论数据集的预期结构。

数据集的结构

在标准的有监督(通常也是无监督)任务中，数据集应该是平衡的。换句话说，属于每个类的样本数量应该几乎相同。在本章我们将要讨论的任务中，我们假设有非常不平衡的数据集 X (包含 N 个样本):

N _异常值 < < N ，如果存在异常值检测(即数据集部分为污垢；(因此，我们需要找出一种方法来过滤掉所有异常值)
N _异常值 = 0 (或者更现实一点， P(N _异常值 > 0) → 0) ，如果有新颖性检测(也就是我们一般可以信任已有样本，把注意力集中在新样本上)

这些标准的原因很明显:让我们考虑前面讨论的例子。如果在 1000000 个加工步骤后观察到的异常率等于 0.2%，则有 2000 个异常，这对于一个工作过程来说可能是一个合理的值。如果这样的数字大得多，就意味着系统中应该有更严重的问题，这超出了数据科学家的角色。因此，在这种情况下，我们期望数据集包含大量正确的样本和非常少的异常(甚至为零)。在许多情况下，经验法则是反映潜在的数据生成过程，因此，如果专家可以确认，例如，有 0.2%的异常，比率应该是 1000÷2 ，以找出一个现实的概率密度函数。在这种情况下，事实上，找出决定异常值可分辨性的因素更为重要。另一方面，如果我们被要求仅执行新颖性检测(例如:区分有效和恶意网络请求)，数据集必须经过验证，以便不包含异常，但同时反映负责所有可能的有效样本的真实数据生成过程。

事实上，如果正确样本的总体是详尽的，任何与高概率区域的大偏差都足以触发警报。相反，真实数据生成过程的有限区域可能导致假阳性结果(也就是说，没有包含在训练集中的有效样本被错误地识别为异常值)。在最坏的情况下，如果特征被改变(也就是说，异常值被错误地识别为有效样本)，噪声很大的子集也可能确定假阴性。然而，在大多数现实生活中，最重要的因素是样本的数量和收集样本的背景。不言而喻，任何模型都必须用将要测试的相同类型的元素来训练。例如:如果化工厂内的测量是使用低精度仪器进行的，则用高精度仪器收集的测试可能不能代表人口(当然，它们比数据集可靠得多)。因此，在进行分析之前，我强烈建议仔细检查数据的性质，并询问是否所有的测试样本都来自同一个数据生成过程。

我们现在可以引入直方图的概念，这是估计包含观测值的数据集分布最简单的方法。

直方图

找出概率密度函数近似值的最简单方法是基于频率计数。如果我们有一个数据集 X 包含 m 个样本 x _i ∈ ℜ (为简单起见，我们只考虑单变量分布，但该过程对于多维样本完全等效)，我们可以定义 m 和 M 如下:

间隔 (m，M) 可以分成固定数量的 b 箱(其可以具有相同或不同的宽度，表示为 w(b _j ) ，使得n_p(b_j)对应于包含在箱中的样品数量b_jT17】。此时，给定一个测试样本 x _t ，很容易理解，通过检测包含 x _t 的仓，并使用以下公式，可以很容易地获得概率的近似值:

在分析这种方法的利弊之前，让我们考虑一个简单的例子，它基于细分为 10 个不同阶层的人的年龄分布:

import numpy as np

nb_samples = [1000, 800, 500, 380, 280, 150, 120, 100, 50, 30]

ages = []

for n in nb_samples:
    i = np.random.uniform(10, 80, size=2)
    a = np.random.uniform(i[0], i[1], size=n).astype(np.int32)
    ages.append(a)

ages = np.concatenate(ages)

The dataset can be reproduced only using the random seed 1000 (that is, setting np.random.seed(1000)).

ages数组包含了所有的样本，我们想创建一个直方图来初步了解分布。我们将使用 NumPy np.histrogram()功能，它提供了所有必要的工具。首先要解决的问题是找出最佳的箱数。对于标准分布来说，这可能很容易，但是当没有关于概率密度的先验知识时，这可能变得极其困难。原因很简单:因为我们需要用一个逐步函数来逼近一个连续函数，所以仓的宽度决定了最终的精度。例如:如果密度是平坦的(例如:均匀分布)，几个面元就足以达到良好的效果。相反，当有峰值时，在函数的一阶导数较大的区域放置更多(更短)的面元，当导数接近零时放置更小的面元(表示平坦区域)会有所帮助。正如我们将要讨论的，使用更复杂的技术，这个过程变得更容易，而直方图通常基于对最佳箱数的更粗略的计算。具体来说，NumPy 允许设置bins='auto'参数，该参数强制算法根据定义明确的统计方法(基于弗里德曼双精度估计器和斯特奇公式)自动选择数字:

在上式中，四分位数区间 ( IQR )对应于第 75 ^个和第 25 ^个个百分点之间的差值。由于我们对分布没有一个清晰的概念，我们更喜欢依靠自动选择，如下面的代码片段所示:

import numpy as np

h, e = np.histogram(ages, bins='auto')

print('Histograms counts: {}'.format(h))
print('Bin edges: {}'.format(e))

上一个片段的输出如下:

Histograms counts: [177  86 122 165 236 266 262 173 269 258 241 116 458 257 311   1   1   5 6]
Bin edges: [16\.         18.73684211 21.47368421 24.21052632 26.94736842 29.68421053
 32.42105263 35.15789474 37.89473684 40.63157895 43.36842105 46.10526316
 48.84210526 51.57894737 54.31578947 57.05263158 59.78947368 62.52631579
 65.26315789 68\.        ]

因此，该算法定义了 19 个仓，并输出了频率计数和边缘(即最小值为16，最大值为68)。我们现在可以显示直方图:

Histogram of the test distribution

该图证实分布相当不规则，一些区域有被平坦区域包围的峰。如前所述，当查询基于样本属于特定箱的概率时，直方图是有帮助的。例如，在这种情况下，我们可能有兴趣确定一个人的年龄在 48.84 到 51.58 之间的概率(这对应于从 0 开始的第 12 个^{箱)。由于所有的箱柜宽度相同，我们可以简单地用n_p(b₁₂**)(h[12])和 m ( ages.shape[0])之间的比值来近似这样一个值:}

d = e[1] - e[0]
p50 = float(h[12]) / float(ages.shape[0])

print('P(48.84 < x < 51.58) = {:.2f} ({:.2f}%)'.format(p50, p50 * 100.0))

输出如下:

P(48.84 < x < 51.58) = 0.13 (13.43%)

因此，概率的近似值约为 13.5%，这也由直方图的结构所证实。然而，读者应该已经清楚地理解，这样的方法有明显的局限性。第一个也是最明显的是箱子的数量和宽度。事实上，一个小数字产生的粗略结果不能考虑快速振荡。另一方面，大量的数据会被驱动到一个带孔的直方图，因为大部分的箱子都没有样本。因此，考虑到现实生活中可能遇到的所有可能的动态，需要一种更可靠的方法。这是我们将在下一节讨论的内容。

核密度估计(KDE)

直方图不连续问题的解决可以用简单的方法有效地解决。给定一个样本 x _i ∈ X ，可以考虑一个超体积(通常是超立方体或超球)，假设我们正在处理多元分布，其中心是 x _i 。这样一个区域的范围是通过一个称为带宽的常数 h 来定义的(选择这个名称是为了支持值为正的有限区域的含义)。然而，我们现在不是简单地计算属于超体积的样本数量，而是使用平滑核函数 K(x _i 来近似这个值；h) 具有一些重要特征:

此外，出于统计和实际的原因，还需要强制执行以下积分约束(为简单起见，它们仅在单变量情况下显示，但扩展很简单):

在讨论名为核密度估计 ( KDE )的技术之前，展示一下*K()*的一些常见选择是有帮助的。

高斯核

这是最常用的内核之一，其结构如下:

下图显示了图形表示:

Gaussian kernel

鉴于其规律性，高斯核是许多密度估计任务的常见选择。但是，由于该方法不允许混合不同的内核，因此选择时必须考虑所有属性。从统计数据中，我们知道高斯分布可以被认为是峰度的平均参考(峰度与峰值和尾部的重量成正比)。为了最大化内核的选择性，我们需要减少带宽。这意味着即使是最小的振荡也会改变密度，结果是非常不规则的估计。另一方面，当 h 较大(即高斯的方差)时，近似变得非常平滑，并且可能失去捕获所有峰值的能力。因此，在选择最合适的带宽的同时，考虑其他可以自然简化过程的内核也是有帮助的。

epanechnikov 核

这个核被提出来最小化均方误差，并且它还具有非常规则的性质(实际上，它可以被想象成倒抛物线)。公式如下:

引入常数 ε 对内核进行归一化，满足所有要求(类似的方式，可以在范围上扩展内核( -h 、 h )以便与其他函数更加一致)。下图显示了图形表示:

Epanechnikov kernel

h → 0 时，籽粒会变得非常尖。然而，鉴于其数学结构，它将始终保持非常规则；因此，在大多数情况下，没有必要将其用作高斯核的替代(即使后者具有稍大的均方误差)。此外，由于函数在x = h(K(x；h) = 0 对于 |x| > h ，它会导致密度估计快速下降，特别是在边界处，例如高斯函数非常平缓地下降。

指数核

一个非常峰值的内核是指数内核，其一般表达式如下:

与高斯核相反，这个核有非常重的尾部和一个尖锐的峰值。下图显示了一个图表:

Exponential kernel

可以看到，这样的函数适用于模拟密度高度集中在某些特定点周围的非常不规则的分布。另一方面，当数据生成过程非常规则且表面光滑时，误差会变得非常高。可以用来评估内核(和带宽)性能的一个很好的理论度量是平均积分平方误差 ( MISE )，定义如下:

在之前的公式中， p _K (x) 为估算密度， p(x) 为实际密度。遗憾的是， p(x) 未知(否则，我们不需要任何估计)；因此，这样的方法只能用于理论评估(例如:Epanechnikov 核的最优性)。然而，很容易理解，每当内核不能靠近实际表面时，MISE 就会变大。由于指数型的跳跃非常突然，它只能在特定的情况下适用。在所有其他情况下，它的行为导致更大的 MISEs，因此其他内核是更好的。

均匀核

这是最简单且不太平滑的内核函数，其用法类似于构建直方图的标准过程。它等于以下内容:

显然，在带宽限定的范围内，这是一个恒定的步长，只有在估计不需要平滑时才有帮助。

估计密度

一旦选择了一个核函数，就有可能使用 k 近邻方法来构建概率密度函数的完全近似。事实上，给定一个数据集 X (为简单起见， X ∈ ℜ ^m ，因此这些值是实数)，很容易创建例如一个球树(如第 2 章、聚类基础中所讨论的)来以有效的方式划分数据。当数据结构准备好时，可以获得带宽定义的半径内查询点 x _j 的所有邻居。假设这样一套是 X _j = {x ₁ ，...，x _t } 且点数为 N _j 。概率密度的估计如下:

不难证明，如果带宽选择得当(作为邻域包含的样本数的函数)， p _K 在概率上收敛于实际的 p(x) 。换句话说，如果粒度足够大，近似和真实密度之间的绝对误差收敛到零。 p _K (x _j ) 的构建过程如下图所示:

Density estimation of x_j. The Kernel functions are evaluated in each point belonging to the neighborhood of x_j

在这一点上，很自然地会问为什么不对每个查询使用整个数据集，而不是 k-NN 方法？答案很简单，它基于这样的假设:在 x _j 计算的密度函数值可以很容易地使用局部行为进行插值(也就是说，对于多元分布，以 x _j 为中心的球)和远点对估计没有影响。因此，我们可以将计算限制在 X 的较小子集内，避免包含接近于零的贡献。

在讨论如何确定最佳带宽之前，让我们展示一下之前定义的数据集的密度估计(使用 scikit-learn)。由于我们没有任何特定的先验知识，我们将采用不同带宽(0.1、0.5 和 1.5)的高斯核。所有其他参数保持默认值；但是KernelDensity类允许设置度量(默认为metric='euclidean')、数据结构(默认为algorithm='auto'，根据维度在球树和 kd 树之间执行自动选择)以及绝对和相对容差(分别为 0 和 10 ^-8 )。在许多情况下，没有必要更改默认值；但是，对于具有特定特征的非常大的数据集，例如，更改leaf_size参数以提高性能可能会有所帮助(如第 2 章、聚类基础中所述)。此外，默认度量不能适用于所有任务(例如:标准文档显示了一个基于哈弗斯距离的示例，在处理纬度和经度时可以使用该示例)。其他情况下，最好使用超立方体，而不是球(曼哈顿距离就是这种情况)。

让我们从实例化类和拟合模型开始:

from sklearn.neighbors import KernelDensity

kd_01 = KernelDensity(kernel='gaussian', bandwidth=0.1)
kd_05 = KernelDensity(kernel='gaussian', bandwidth=0.5)
kd_15 = KernelDensity(kernel='gaussian', bandwidth=1.5)

kd_01.fit(ages.reshape(-1, 1))
kd_05.fit(ages.reshape(-1, 1))
kd_15.fit(ages.reshape(-1, 1))

此时，可以调用score_samples()方法来获得一组数据点的对数密度估计值(在我们的例子中，我们考虑的是 0.05 增量的范围(10，70))。由于数值为 log(p) ，需要计算 e ^log(p) 才能得到实际概率。

结果图如下图所示:

Gaussian density estimations with bandwidths: 0.1 (top), 0.5 (middle), and 1.5 (bottom)

可以注意到，当带宽很小时(0.1)，由于缺少特定子范围的样本，密度会有强烈的振荡。当 h = 0.5 时，轮廓(由于数据集是单变量的)变得更加稳定，但仍有一些由邻域内部方差引起的残留快速变化。当 h 变大(在我们的例子中为 1.5)时，这种行为几乎完全消除。一个显而易见的问题是:如何确定最合适的带宽？当然，最自然的选择是最小化 MISE 的 h 的值，但是，正如所讨论的，这种方法只能在真实概率密度已知的情况下使用。然而，有几个经验标准已经被证实是非常可靠的。给定完整的数据集 X ∈ ℜ ^m ，第一个数据集基于以下公式:

在我们的案例中，我们获得了以下信息:

import numpy as np

N = float(ages.shape[0])
h = 1.06 * np.std(ages) * np.power(N, -0.2)

print('h = {:.3f}'.format(h))

输出如下:

h = 2.415

因此，建议增加的带宽甚至比我们上次实验中的还要多。因此，第二种方法基于四分位数范围( IQR = Q3 - Q1 或相当于第 75 ^个百分位-第 25 ^个百分位)，它对非常强的内部变化更加稳健:

计算如下:

import numpy as np

IQR = np.percentile(ages, 75) - np.percentile(ages, 25)
h = 0.9 * np.min([np.std(ages), IQR / 1.34]) * np.power(N, -0.2)

print('h = {:.3f}'.format(h))

现在的输出是:

h = 2.051

这个值比前一个值要小得多，说明 p _K (x) 可以用更小的超体积更精确。根据经验，我建议选择带宽最小的方法，即使第二种方法通常在不同的上下文中提供最佳结果。现在让我们使用 h = 2.0 和高斯核、Epanechnikov 核和指数核(我们排除了均匀核，因为最终结果相当于直方图)重新执行估算:

from sklearn.neighbors import KernelDensity

kd_gaussian = KernelDensity(kernel='gaussian', bandwidth=2.0)
kd_epanechnikov = KernelDensity(kernel='epanechnikov', bandwidth=2.0)
kd_exponential = KernelDensity(kernel='exponential', bandwidth=2.0)

kd_gaussian.fit(ages.reshape(-1, 1))
kd_epanechnikov.fit(ages.reshape(-1, 1))
kd_exponential.fit(ages.reshape(-1, 1))

图形输出如下图所示:

Density estimations with bandwidths equal to 2.0 and Gaussian kernel (top), Epanechnikov kernel (middle), and Exponential kernel (bottom)

不出所料，Epanechnikov 核和指数核都比高斯核更具振荡性(因为当 h 小时，它们往往更具峰值)；然而，很明显，中心情节肯定是最准确的(就 MISE 而言)。类似的结果以前已经用高斯核和 h = 0.5 得到，但是，在那种情况下，振荡是极不规则的。如上所述，当值达到带宽边界时，Epanechnikov 内核具有非常强的不连续趋势。这种现象可以通过观察估计值的极值立即理解，估计值几乎垂直下降到零。相反， h = 2 的高斯估计似乎很平滑，没有捕捉到 50 年到 60 年之间的变化。同样的情况也发生在指数内核上，这也显示了它的特殊行为:非常尖锐的极端。在下面的例子中，我们将使用 Epanechnikov 内核；然而，我邀请读者也检查不同带宽的高斯结果。这种选择有一个精确的理由(没有充分的理由不能放弃):我们认为数据集是详尽的，我们想惩罚所有克服自然极端的样本。在所有其他情况下，非常小的剩余概率是优选的；然而，这种选择必须考虑到每一个具体目标。

异常检测

现在让我们应用 Epanechnikov 密度估计来执行异常检测的示例。根据概率密度的结构，我们决定在 p(x) < 0.005 处实施截止。这种情况显示在下面的屏幕截图中:

Epanechnikov density estimation with anomaly cut-off

红点表示样本被归类为异常的年龄限制。让我们计算一些测试点的概率密度:

import numpy as np

test_data = np.array([12, 15, 18, 20, 25, 30, 40, 50, 55, 60, 65, 70, 75, 80, 85, 90]).reshape(-1, 1)

test_densities_epanechnikov = np.exp(kd_epanechnikov.score_samples(test_data))
test_densities_gaussian = np.exp(kd_gaussian.score_samples(test_data))

for age, density in zip(np.squeeze(test_data), test_densities_epanechnikov):
    print('p(Age = {:d}) = {:.7f} ({})'.format(age, density, 'Anomaly' if density < 0.005 else 'Normal'))

上一个片段的输出如下:

p(Age = 12) = 0.0000000 (Anomaly)
p(Age = 15) = 0.0049487 (Anomaly)
p(Age = 18) = 0.0131965 (Normal)
p(Age = 20) = 0.0078079 (Normal)
p(Age = 25) = 0.0202346 (Normal)
p(Age = 30) = 0.0238636 (Normal)
p(Age = 40) = 0.0262830 (Normal)
p(Age = 50) = 0.0396169 (Normal)
p(Age = 55) = 0.0249084 (Normal)
p(Age = 60) = 0.0000825 (Anomaly)
p(Age = 65) = 0.0006598 (Anomaly)
p(Age = 70) = 0.0000000 (Anomaly)
p(Age = 75) = 0.0000000 (Anomaly)
p(Age = 80) = 0.0000000 (Anomaly)
p(Age = 85) = 0.0000000 (Anomaly)
p(Age = 90) = 0.0000000 (Anomaly)

可以看到，函数的突然下降产生了一种垂直分离。一个年龄15的人，几乎到了边界( p(15) ≈ 0.0049 ，而对于上限来说，行为更为激烈。截止年龄约为 58 岁，但年龄为60的样本比年龄为 57 岁的样本可能性小 10 倍左右(这也由初始直方图证实)。由于这只是一个说教的例子，很容易发现异常；然而，如果没有标准化的算法，即使是稍微复杂一点的分布也会产生一些问题。特别是在这种特殊情况下，这是一个简单的单变量分布，异常通常位于尾部。

因此，我们假设给定总体密度估计 p _K (x) :

当考虑包含所有样本(正常样本和异常样本)的数据集时，这种行为通常是不正确的，数据科学家在决定阈值时必须小心。即使可以很明显，通过去除数据集中的所有异常来学习正态分布也是一个好主意，这样可以将异常所在的区域变平( p _K (x) → 0 )。这样，前面的标准仍然有效，并且很容易比较不同的密度来进行区分。

在继续下一个例子之前，我建议通过创建人工洞和为检测设置不同的阈值来修改初始分布。此外，我邀请读者基于例如年龄和身高生成二元分布(例如:基于一些高斯分布的总和)，并创建一个简单的模型，该模型能够检测其参数非常不可能的所有人。

利用 KDD 杯 99 数据集进行异常检测

此示例基于 KDD 杯 99 数据集，该数据集收集了一长串正常和恶意的互联网活动。特别是，我们将重点关注 HTTP 请求的子集，它有四个属性:持续时间、源字节、目标字节和行为(这更像是一个分类元素，但它有助于我们立即访问一些特定的攻击)。由于原始值是零附近的非常小的数字，所有版本(包括 scikit-learn one)都使用公式 log(x + 0.1) 对变量进行了重整(因此，在用新样本模拟异常检测时必须应用它)。当然，逆变换如下:

让我们从使用 scikit-learn 内置函数fetch_kddcup99()加载和准备数据集开始，选择percent10=True将数据限制在原始集合的 10%(非常大)。当然，我邀请读者也用整个数据集和完整的参数列表(包含 34 个数值)进行测试。

在这种情况下，我们还选择了subset='http'，它已经准备好包含非常多的正常连接和一些特定的攻击(如在标准的定期日志中):

from sklearn.datasets import fetch_kddcup99

kddcup99 = fetch_kddcup99(subset='http', percent10=True, random_state=1000)

X = kddcup99['data'].astype(np.float64)
Y = kddcup99['target']

print('Statuses: {}'.format(np.unique(Y)))
print('Normal samples: {}'.format(X[Y == b'normal.'].shape[0]))
print('Anomalies: {}'.format(X[Y != b'normal.'].shape[0]))

输出如下:

Statuses: [b'back.' b'ipsweep.' b'normal.' b'phf.' b'satan.'] Normal samples: 56516 Anomalies: 2209

因此，有四种类型的攻击(其细节在本文中并不重要)具有2209恶意样本和56516正常连接。为了进行密度估计，我们将把这三个分量作为独立的随机变量(这并不完全正确，但它可以是一个合理的起点)进行一些初步考虑，但最终估计是基于完全联合分布的。当我们想要确定最佳带宽时，让我们执行一个基本的统计分析:

import numpy as np

means = np.mean(X, axis=0)
stds = np.std(X, axis=0)
IQRs = np.percentile(X, 75, axis=0) - np.percentile(X, 25, axis=0)

上一个片段的输出如下:

Means: [-2.26381954  5.73573107  7.53879208]
Standard devations: [0.49261436 1.06024947 1.32979463]
IQRs: [0\.         0.34871118 1.99673381]

持续时间(第一个分量)的 IQR 为空；因此，大多数价值观是平等的。让我们绘制一个直方图来证实这一点:

Histogram for the first component (duration)

不出所料，这样的成分并不十分显著，因为只有一小部分样本具有不同的值。因此，在这个例子中，我们将跳过它，只处理源和目标字节。现在让我们按照前面的解释计算带宽:

import numpy as np

N = float(X.shape[0])

h0 = 0.9 * np.min([stds[0], IQRs[0] / 1.34]) * np.power(N, -0.2)
h1 = 0.9 * np.min([stds[1], IQRs[1] / 1.34]) * np.power(N, -0.2)
h2 = 0.9 * np.min([stds[2], IQRs[2] / 1.34]) * np.power(N, -0.2)

print('h0 = {:.3f}, h1 = {:.3f}, h2 = {:.3f}'.format(h0, h1, h2))

输出如下:

h0 = 0.000, h1 = 0.026, h2 = 0.133

排除第一个值，我们需要在h1和h2之间进行选择。由于这些值的幅度不大，我们希望有很大的选择性，我们将设置 h = 0.025 并采用高斯核，这提供了良好的平滑度。包含第一个组件的分割输出(使用 seaborn 可视化库获得，该库包括一个内部 KDE 模块)如下图所示:

Density estimations for normal connections (upper line) and malicious attacks (lower line)

第一行显示正常连接的密度，而较低的是恶意攻击。正如预期的那样，第一个组件(持续时间)在两种情况下几乎相同，可以丢弃。相反，源字节和目标字节表现出非常不同的行为。不考虑对数变换，正常连接平均发送 5 个字节，方差非常低，将潜在范围扩展到区间(4，6)。响应的方差较大，值在 4 到 10 之间，密度从 10 开始非常低。相反，恶意攻击的源字节和目标字节都有两个峰值:较短的字节对应于-2，较高的字节分别对应于约 11 和 9(与正常区域的重叠最小)。即使不考虑完整的联合概率密度，也不难理解大多数攻击发送更多的输入数据，接收更长的响应(而连接持续时间不受强烈影响)。

我们现在可以通过仅选择正常样本(即对应于Y == b'normal.')来训练估计器:

from sklearn.neighbors import KernelDensity

X = X[:, 1:]

kd = KernelDensity(kernel='gaussian', bandwidth=0.025)
kd.fit(X[Y == b'normal.'])

让我们计算正常和异常样本的密度:

Yn = np.exp(kd.score_samples(X[Y == b'normal.']))
Ya = np.exp(kd.score_samples(X[Y != b'normal.']))

print('Mean normal: {:.5f} - Std: {:.5f}'.format(np.mean(Yn), np.std(Yn)))
print('Mean anomalies: {:.5f} - Std: {:.5f}'.format(np.mean(Ya), np.std(Ya)))

输出如下:

Mean normal: 0.39588 - Std: 0.25755
Mean anomalies: 0.00008 - Std: 0.00374

很明显，当例如 p _K (x) < 0.05 (考虑三个标准差)对于一个异常，我们得到p_K(x)∞(0，0.01)) ，而Yn的中位数约为 0.35。这意味着至少一半的样品具有p_K**(x)>0.35。然而，经过简单的计数检查后，我们得到了以下结果:

print(np.sum(Yn < 0.05))
print(np.sum(Yn < 0.03))
print(np.sum(Yn < 0.02))
print(np.sum(Yn < 0.015))

输出如下:

由于有 56，516 个正常样本，我们可以决定选择两个阈值(同时考虑异常异常值):

正常连接:p_K(x)>0.03
中度预警 : 0.03(涉及 3.1%的正常样本可识别为假阳性)
高警戒 : 0.015(这种情况下，只有 1.2%的正常样本可以触发报警)

此外，在第二个警报中，我们发现:

print(np.sum(Ya < 0.015))

输出如下:

因此，只有一个异常样本具有 p _K (x) > 0.015 (有 2，209 个向量)，这证实了这样的选择是合理的。密度直方图也证实了之前的结果:

Histogram of anomaly (left) and normal (right) densities

正态分布的右尾并不惊人，因为异常高度集中在左侧。在这个领域，也有大部分的异常，因此也是最关键的。原因与特定的域(对于不同类型的请求，输入和输出字节可能非常相似)密切相关，在更稳定的解决方案中，有必要考虑进一步的参数(例如:完整的 KDD 杯 99 数据集)。但是，出于教学目的，我们可以定义一个简单的函数(基于之前定义的阈值)来根据源和目标字节量(非对数)检查连接状态:

import numpy as np

def is_anomaly(kd, source, destination, medium_thr=0.03, high_thr=0.015):
    xs = np.log(source + 0.1)
    xd = np.log(destination + 0.1)
    data = np.array([[xs, xd]])

    density = np.exp(kd.score_samples(data))[0]

    if density >= medium_thr:
        return density, 'Normal connection'
    elif density >= high_thr:
        return density, 'Medium risk'
    else:
        return density, 'High risk'

我们现在可以用三个不同的例子来测试这个函数:

print('p = {:.2f} - {}'.format(*is_anomaly(kd, 200, 1100)))
print('p = {:.2f} - {}'.format(*is_anomaly(kd, 360, 200)))
print('p = {:.2f} - {}'.format(*is_anomaly(kd, 800, 1800)))

输出如下:

p = 0.30 - Normal connection
p = 0.02 - Medium risk
p = 0.00000 - High risk

一般来说，还可以考虑源和目标字节密度的二元图:

Bivariate plot of the source and destination bytes densities

前面的截图证实，虽然攻击通常涉及大量的输入字节，但响应与正常的非常相似，即使它们占据了该区域的极端部分。作为练习，我邀请读者用整个 KDD 杯 99 数据集训练一个模型，并找出检测非常危险和中等风险攻击的最佳阈值。

一类支持向量机

schlkopf B，Platt J C，Shawe-Taylor J C，Smola A J 和 Williamson R C 在文章估计高维分布的支持度，神经计算，2001 年 13 月 7 日中提出了单类支持向量机的概念，作为将新奇事物分类为从真实数据生成过程中抽取的样本或离群值的方法。让我们从我们想要实现的目标开始:找到一个无监督的模型，在给定样本 x _i 的情况下，该模型可以产生二进制输出 y _i (传统上，SVMs 结果是双极性的:-1 和+1)，因此，如果 x _i 是内联的 y _i = +1 ，反之， y _i = -1 如果 x _i 是一个离群值(更正确的说法是，在前面提到的论文中，作者假设对于构成训练集的大多数内联者来说，结果是 1 )。乍一看，这似乎是一个经典的监督问题；然而，这并不是因为它不需要有标签的数据集。事实上，给定包含 m 样本x_I∈ℜⁿt41】的数据集 X ，将使用单个固定类训练模型，目标是找到一个分离超平面，使 X 与原点之间的距离最大化。首先，让我们考虑一个简单的线性情况，如下图所示:

Linear one-class SVM scenario: the training set is separated from the origin with the largest margin

该模型被训练以找出超平面的参数，这些参数使到原点的距离最大化。超平面一侧的所有样本预计都是内联的，输出标签为 +1 ，其余所有样本都被认为是外联的，输出标签为 -1 。这个标准看起来很有效，但它只适用于线性可分离的数据集。标准支持向量机通过将数据集(通过函数φ())投影到特征空间 D 来解决这个问题，在特征空间 D 中它获得这样的属性:

特别是，考虑到问题的数学本质，如果选择一个内核，投影在计算上就变得轻量级了。换句话说，我们希望使用具有以下属性的函数:

投影函数φ()的存在性保证在一个非常容易得到的条件(称为 Mercer 条件)下存在(即在一个实子空间中，核必须是正半定的)。这样选择的原因与问题解决所涉及的过程有严格的关系(更详细的解释可以在机器学习算法第二版、博纳科索格、 Packt Publishing ，2018 中找到)。但是，不熟悉 SVMs 的读者不应该担心，因为我们不打算讨论太多的数学细节。要记住的最重要的事情是，不支持任何内核的通用投影会导致计算复杂性的急剧增加(特别是对于大型数据集)。

*K(，)*最常见的选择之一是径向基函数(已经在第三章、高级聚类中进行了分析):

另一个有用的内核是多项式内核:

在这种情况下，指数 c 定义多项式函数的次数，该次数与特征空间的维数成比例。然而，内核及其超参数的选择都是依赖于上下文的，并且没有总是有效的通用规则。因此，对于每一个问题，都需要进行初步分析，通常还需要进行网格搜索，以做出最合适的选择。一旦选择了内核，问题可以用下面的方式表示:

没有完整的讨论(这超出了本书的范围)，我们可以将注意力集中在几个重要的元素上。首先，决策功能如下:

解决方案中涉及的数学过程允许我们简化以下表达式，但是，出于我们的目的，最好保留原始表达式。如果读者具有监督学习的基本知识，他们可以容易地理解权重向量和样本的投影之间的点积 x _i 允许确定 x _i 相对于分离超平面的位置。事实上，如果两个向量之间的角度小于 90°(π/2)，点积是非负的。当角度正好为 90°时，它等于零(也就是说，向量是正交的)，当角度在 90°和 180°之间时，它是负的。该过程如下图所示:

Decision process in SVM

权重向量与分离超平面正交。样本 x _i 被识别为内联，因为点积为正且大于阈值 ρ 。相反， x _j 被标记为异常值，因为决策函数的符号为负。术语*ξ_I(ξ_I≥0)*被称为松弛变量，并被引入以允许在离群值和内联值之间有更灵活的边界。事实上，如果这些变量都等于零(为了简单起见，还有 ρ=1 )，那么施加在优化问题上的条件就变成了:

这意味着所有训练样本必须被视为内联样本，因此，必须选择分离超平面，使得所有 x _i 都在同一侧。然而，通过定义软边界，松弛变量的使用允许更大的灵活性。每个训练样本都与一个变量 ξ _i 相关联，当然，这个问题要求它们最小化。然而，通过这个技巧，一些边界样本也可以位于超平面的相对侧(足够靠近它)，即使它们继续被识别为内联。最后一个要考虑的元素是本文中最重要的元素之一，涉及超参数ν∞(0，1) 。在前述论文中，作者证明了每当 ρ ≠ 0 时，ν 可以解释为训练样本分数的上界，实际上是离群值。在这一章的开头，我们已经说明了在一个新奇的检测问题中，数据集必须是干净的。不幸的是，这并不总是正确的；因此， ν 和松弛变量的联合使用也允许我们处理包含一小部分异常值的数据集。在概率方面，如果 X 是从被噪声部分破坏的数据生成过程中提取的， ν 是在 X 中找到异常值的概率。

现在，让我们分析一个基于用元组(年龄、身高)标识的学生数据集的二维示例。我们将从二元高斯分布和 200 个均匀采样的测试点生成 2000 个内联:

import numpy as np

nb_samples = 2000
nb_test_samples = 200

X = np.empty(shape=(nb_samples + nb_test_samples, 2))

X[:nb_samples] = np.random.multivariate_normal([15, 160], np.diag([1.5, 10]), size=nb_samples)
X[nb_samples:, 0] = np.random.uniform(11, 19, size=nb_test_samples)
X[nb_samples:, 1] = np.random.uniform(120, 210, size=nb_test_samples)

由于比例不同，最好在训练模型之前对数据集进行归一化:

from sklearn.preprocessing import StandardScaler

ss = StandardScaler()
Xs = ss.fit_transform(X)

下图显示了标准化数据集的图表:

Dataset for the one-class SVM example

主斑点主要由内联体组成，部分测试样本位于同一高密度区域；因此，我们可以合理地假设在包含所有样本的数据集中约有 20%的异常值(因此， ν=0.2 )。当然，这样的选择是基于我们的一个假设，在任何现实场景中， ν 的值必须始终反映数据集中预期异常值的实际百分比。当这条信息不可用时，最好从一个较大的值开始(例如， ν=0.5 )，并通过减少它来继续，直到找到最佳配置(即，误分类的概率最小)。

同样重要的是要记住，训练过程有时会找到次优解；因此，一些内联可能被标记为异常值。在这些情况下，最好的策略是测试不同核的效果，例如，当使用多项式核时，增加它们的复杂性，直到找到最优解(不一定排除所有错误)。

现在，让我们使用一个径向基函数内核(特别适合高斯数据生成过程)初始化 scikit-learn OneClassSVM类的一个实例，并训练模型:

from sklearn.svm import OneClassSVM

ocsvm = OneClassSVM(kernel='rbf', gamma='scale', nu=0.2)
Ys = ocsvm.fit_predict(Xs)

我们选择了推荐值gamma='scale'，该值基于以下公式:

这种选择通常是最佳的起点，可以改变(根据结果是否不可接受而增加或减少)。在我们的例子中，由于数据集是二维的( n=2 )和归一化的 (std(X) = 1) ， γ = 0.5 ，这对应于单位方差高斯分布(因此，我们应该期望它是最合适的选择)。此时，我们可以通过突出显示异常值来绘制结果:

Classification result (left). Outliers from the test set (right)

从左图中可以看到，该模型已经成功识别了数据集的高密度部分，并且还将密集斑点外部区域的一些样本标记为异常值。它们对应于二元高斯下的低概率值，在我们的例子中，我们假设它们是应该被滤除的噪声样本。在右图中，只能看到异常区域，当然，这是高密度斑点的补充。我们可以总结说，一级 SVM，即使有点容易过度，也能够帮助我们以非常小的错误概率识别新奇事物。这也是由于数据集的结构(然而，这在许多情况下非常常见)，可以使用径向基函数内核轻松管理。不幸的是，对于高维数据，这种简单性往往会丧失，需要更彻底的超参数搜索来最小化错误率。

利用隔离森林进行异常检测

刘 F T，Ting K M，周 z 在文章隔离森林，，第八届 IEEE 国际数据挖掘会议， 2008)中提出了一种非常强大的异常检测方法，它基于集成学习的一般框架。由于这个主题非常广泛，并且主要涵盖在有监督的机器学习书籍中，我们邀请读者在必要时查看建议的资源之一。相反，在这种情况下，我们将描述模型，而不是非常强烈地引用所有潜在的理论。

先说一个森林是一组独立的模型叫做决策树。顾名思义，它们不仅仅是算法，而是划分数据集的一种非常实用的方法。从根开始，对于每个节点，选择一个特征和一个阈值，并将样本分成两个子集(对于非二叉树不是这样，但一般来说，所有涉及的树都是这些模型是二进制的)，如下图所示:

Generic structure of a binary decision tree

在有监督的任务中，选择元组(特征、阈值)是根据特定的标准来选择的，这些标准最小化了子对象的杂质。这意味着目标通常是分割一个节点，使得结果子集包含属于单个类的大多数样本。当然，这很容易理解，当所有的叶子都是纯的或者达到最大深度时，这个过程就结束了。相反，在这个特定的上下文中，我们从一个非常特殊的(但经验证明的)假设开始:如果属于隔离森林的树是在每次选择随机特征和随机阈值时生长的，则从根到包含任何内联的叶子的路径的平均长度比隔离离群点所需的路径长。这一假设的原因可以通过考虑一个二维例子很容易理解，如作者所示:

Bidimensional random partitioning. On the left, an inlier is isolated. On the right, an outlier belonging to a low-density region is detected

可以观察到，内联体通常属于高密度区域，需要更多的分区来隔离样本。相反，位于低密度区域的异常值可以用较少的划分步骤来检测，因为所需的粒度与斑点的密度成比例。因此，构建隔离林的目标是测量所有内联器的平均路径长度，并将其与新样本所需的路径长度进行比较。当这样的长度较短时，成为异常值的概率增加。作者提出的异常分数基于指数函数:

在前面的公式中， m 是属于训练集 X ，avg(h(X_I)是考虑所有树的样本 x _i 的平均路径长度， c(m) 是只依赖于 m 的归一化项。当 s(x _i ，m) → 1 时，样品x_IT19】被确定为异常。因此，由于s(**)被限制在 0 和 1 之间，如果我们考虑 0.5 的阈值，则正常样本与值 s(x _i ，m) < < 0.5 相关联。

现在让我们考虑一下葡萄酒数据集，它包含 178 个样本 x _i ∈ ℜ ¹³ ，其中每个特征都是特定的化学属性(例如，酒精、苹果酸、灰分等)，并训练一个隔离森林来检测一种新的葡萄酒是可以被认为是内联的(例如，现有品牌的变体)还是离群的，因为它的化学属性与每个现有样本都不兼容。第一步包括加载和规范化数据集:

import numpy as np

from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler

wine = load_wine()
X = wine['data'].astype(np.float64)

ss = StandardScaler()
X = ss.fit_transform(X)

我们现在可以实例化一个IsolationForest类并设置最重要的超参数。第一个是n_estimators=150，通知模型训练 150 棵树。另一个基本参数(类似于一类支持向量机中的 ν )称为contamination，其值表示训练集中异常值的预期百分比。因为我们信任数据集，所以我们选择了一个等于 0.01 (1%)的值来处理可以忽略不计的奇怪样本。出于兼容性原因，插入了behaviour='new'参数(查看官方文档了解更多信息)random_state=1000保证实验的可重复性。一旦类被初始化，就可以训练模型:

from sklearn.ensemble import IsolationForest

isf = IsolationForest(n_estimators=150, behaviour='new', contamination=0.01, random_state=1000)
Y_pred = isf.fit_predict(X)

print('Outliers in the training set: {}'.format(np.sum(Y_pred == -1)))

上一个片段的输出是:

因此，隔离林已经成功识别了 178 个内联器中的 176 个。我们可以接受这个结果，但是，像往常一样，我建议调整参数，以便获得一个与每个特定情况兼容的模型。此时，我们可以生成几个噪声样本:

import numpy as np

X_test_1 = np.mean(X) + np.random.normal(0.0, 1.0, size=(50, 13))
X_test_2 = np.mean(X) + np.random.normal(0.0, 2.0, size=(50, 13))
X_test = np.concatenate([X_test_1, X_test_2], axis=0)

测试集被分成两个块。第一个数组X_test_1包含噪声水平相对较低的样本( σ=1 )，而第二个数组X_test_2包含更多噪声样本( σ=2 )。因此，我们预计第一组中的异常值数量较少，第二组中的异常值数量较多。数组X_test是两个测试集的有序连接。现在让我们预测一下状态。由于这些值是双极性的，我们希望将它们与训练结果区分开来，因此我们将乘以预测时间2(即，-1表示训练集中的异常值，1表示训练集中的内联值，-2表示测试集中的异常值，2表示测试集中的内联值):

Y_test = isf.predict(X_test) * 2

Xf = np.concatenate([X, X_test], axis=0)
Yf = np.concatenate([Y_pred, Y_test], axis=0)

print(Yf[::-1])

输出如下:

[ 2 2 -2 -2 -2 -2 -2 2 2 2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 2 -2 2 2 -2 -2 -2 2 -2 -2 -2 -2 2 2 -2 -2 -2 -2 -2 -2 2 2 -2 2 -2 2 -2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 -2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 1 1 1 1 1 1 1 1 1 1 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

随着顺序的保留和反转，我们可以看到属于X_test_2(高方差)的大部分样本被归类为异常，而低方差的大部分样本被识别为内联。为了得到进一步的视觉确认，我们可以执行 t-SNE 降维，考虑到最终结果是二维分布，其与原始(13 维)分布的库尔巴克-莱布勒散度最低。这意味着所得维度的可解释性非常低，并且该图只能用于理解二维空间的哪些区域更有可能被内联体占据:

from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=5, n_iter=5000, random_state=1000)
X_tsne = tsne.fit_transform(Xf)

结果图如下图所示:

t-SNE plot for the novelty detection with the wine dataset

可以看到，许多接近训练内联器的样本本身就是内联器，一般来说，几乎所有远测试样本都是异常值。然而，由于强降维，很难得出更多的结论。然而，我们知道当噪声足够小时，找到内联的概率很大(这是一个合理的结果)。作为练习，我邀请读者检查单个化学属性(可在https://sci kit-learn . org/stable/datasets/index . html # wine-dataset上获得)，并针对每个属性或针对组，找出将内联者转化为外联者的阈值(例如，回答这个问题:与训练集兼容的最大酒精量是多少？).

摘要

在这一章中，我们讨论了概率密度函数的性质，以及如何利用它们来计算实际概率和相对可能性。我们已经看到了如何创建直方图，这是将值分组到预定义的箱中后表示值的频率的最简单方法。由于直方图有一些重要的限制(它们非常不连续，很难找到最佳的面元大小)，我们引入了核密度估计的概念，这是一种使用平滑函数估计密度的稍微复杂的方法。

我们分析了最常见的核(高斯核、Epanechnikov 核、指数核和均匀核)的特性，以及两种可用于找出每个数据集的最佳带宽的经验方法。使用这样的技术，我们试图解决一个基于合成数据集的非常简单的单变量问题。我们分析了 KDD 杯 99 数据集的 HTTP 子集，其中包含了几个正常和恶意网络连接的日志记录。我们使用 KDE 技术创建了一个基于两个阈值的简单异常检测系统，我们还解释了在处理这类问题时必须考虑哪些因素。

在最后一部分，我们分析了两种常见的方法，可以用来进行新颖性检测。单类支持向量机利用核的能力将复杂的数据集投影到特征空间，在那里它们可以线性分离。下一步是基于假设所有的训练集(除了一小部分)都是内联的，因此它们属于同一个类。该模型以最大化内联体和特征空间原点之间的分离为目标进行训练，结果基于样本相对于分离超平面的位置。相反，隔离森林是基于假设的集成模型，即在随机训练的决策树中，从根到样本的路径对于异常值来说平均较短。

因此，在训练森林之后，可以考虑给定新样本的平均路径长度来计算异常分数。当这样的分数接近 1 时，我们就可以断定出现异常的概率也是非常大的。相反，非常小的分值表明新奇事物反而是潜在的内联。

在下一章中，我们将讨论最常见的降维和字典学习技术，当需要管理具有大量特征的数据集时，这些技术非常有用。

问题

一个人身高 1.70 m 的概率是 p(高)= 0.75 ，而明天要下雨的概率是 P(雨)= 0.2 。 *P(高，雨)*的概率是多少？(即一个人身高 1.70 m，明天要下雨的概率)。
给定一个数据集 X ，我们构建一个包含 1000 个面元的直方图，我们发现其中许多面元是空的。为什么会这样？
直方图包含三个面元，分别包含 20、30 和 25 个样本。第一个仓的范围为 0 < x < 2 ，第二个 2 < x < 4 ，第三个 4 < x < 6 。 P(x) > 2 的大概概率是多少？
给定正态分布 N(0，1) ，样本 x 与 p(x) = 0.35 是否可以认为是异常？
有 500 个样本的数据集 X 有标准值(X) = 2.5 和 IQR 值(X) = 3.0 。最佳带宽是多少？
一位专家告诉我们，一个分布在两个值附近达到极大的峰值，密度突然从峰值的平均值下降 0.2 个标准差。哪种内核最合适？
给定一个样本 x (从 10，000 个样本的流动人群中收集)，我们不确定这是异常还是新奇，因为 p(x) = 0.0005 。在另外 10，000 次观察之后，我们重新训练模型并且 x 保持 p(x) < 0.001 。我们能断定 x 是异常吗？

进一步阅读

Epanechnikov V . A，多元概率密度的非参数估计，概率论及其应用，14， 1969
《数理统计年鉴》，1962 年
谢瑟·S·J，六种流行的带宽选择方法在一些真实数据集上的性能(附讨论)，计算统计学，7， 1992
施科尔科夫，普拉特·J·C，肖-泰勒·J·C，斯摩拉·A·J，威廉姆森·R·C，估计高维分布的支持，神经计算，13/7， 2001
刘 F T，丁 K M，周 z，隔离森林，ICDM 2008 ，第八届 IEEE 数据挖掘国际会议，2008
《理论神经科学》，麻省理工学院出版社，2005 年
机器学习算法第二版，博纳科索格，帕克特出版，2018

七、降维和成分分析

在这一章，我们将介绍和讨论一些非常重要的技术，可以用来进行降维和组件提取。在前一种情况下，目标是将高维数据集转换为低维数据集，以尽量减少信息损失。后者是一个过程，需要找到一个可以混合的原子字典，以建立样本。

特别是，我们将讨论以下主题:

主成分分析 ( 主成分分析
奇异值分解 ( 奇异值分解)和白化
核主成分分析
稀疏主成分分析和字典学习
要素分析
独立成分分析 ( ICA
非负矩阵分解 ( NNMF )
潜在狄利克雷分配 ( LDA )

技术要求

本章将介绍的代码需要以下内容:

Python 3.5+(强烈推荐蟒蛇分布(www.anaconda.com/distributio…)
以下库:
- SciPy 0.19+
- NumPy 1.10+
- 学习 0.20+
- 熊猫 0.22+
- Matplotlib 2.0+
- seaborn 0.9+

示例可在 GitHub 资源库中获得，网址为https://GitHub . com/PacktPublishing/HandsOn-Unsupervised-Learning-with-Python/tree/master/chapter 07。

主成分分析

降低数据集维数的最常见方法之一是基于样本协方差矩阵的分析。一般来说，我们知道随机变量的信息含量与其方差成正比。例如，给定一个多变量高斯，熵(我们用来测量信息的数学表达式)如下:

上式中，σ为协方差矩阵。如果我们(不失一般性)假设σ是对角的，那么很容易理解熵(按比例)大于每个单个分量的方差， σ _i ² 。这并不奇怪，因为一个低方差的随机变量集中在均值附近，惊喜的概率很低。另一方面，当 σ ² 变得越来越大时，潜在的结果随着不确定性而增加，这与信息量成正比。

当然，组件的影响一般是不同的；因此，主成分分析 ( PCA )的目标是找到样本的线性变换，能够将它们投影到更低维的子空间上，从而保留最大可能量的初始方差。实际上，让我们考虑一个数据集，x∈ℜ^m×nT7:

我们想要找到的线性变换是一个新的数据集，如下所示:

在应用这样的转换之后，我们期望得到以下结果:

让我们开始考虑样本协方差矩阵(出于我们的目的，我们也可以采用有偏估计)；为简单起见，我们还假设 X 的平均值为零:

这样的矩阵是对称的，是正半定的(不熟悉这些概念也没关系，但是对于证明后面的步骤很重要)，所以它的特征向量构成了正交基。简单回顾一下，如果 A 是正方形矩阵，向量 v _i 被称为与特征值相关联的特征向量 λ _i ，如果下列情况成立:

换句话说，特征向量被转换成其自身的扩展或收缩版本(不能发生旋转)。证明协方差矩阵的特征向量定义协方差分量的方向(即数据集具有特定协方差分量的方向)并不难(但所有数学细节将被省略)。然而，原因很简单；事实上，在变换之后，新的协方差矩阵(变换数据集的， Z )是不相关的(即它是对角的)，因为新的轴与协方差分量对齐。这意味着一个 versor(例如， v ₀ = (1，0，0，...，0) )转化为σ_I²v_I，所以是一个关联特征值与 i ^th 分量方差成正比的特征向量。

因此，为了找出哪些元素可以被丢弃，我们可以对特征值进行排序，从而得出以下结论:

对应的特征向量( v ₁ ，v ₂ ，...，v _n )分别确定最大方差对应的分量，以此类推，直到最后一个。形式上，我们定义这样的特征向量为主成分；因此，第一主成分是与最大方差相关的方向，第二主成分与第一主成分正交，并且它与第二大方差相关，以此类推。在二维数据集的情况下，这个概念显示在下面的屏幕截图中:

Principal components of a bidimensional dataset; the first principal component lies along the axis with the largest variance, while the second one is orthogonal, and it's proportional to the residual variance

至此，问题差不多解决了；事实上，如果我们只选择第一个 k 主成分( v _i ∈ ℜ ^{n × 1} )，我们就可以构建一个变换矩阵，a_k∈ℜ^n×k，从而将与第一个 k 特征值关联的特征向量作为行:

因此，我们可以通过使用以下矩阵乘法来转换整个数据集:

新数据集 Z 的维度等于 k <(或< < ) n ，并且它包含与组件数量成比例的原始方差量。例如，考虑到上一个截图中显示的例子，如果我们选择单个分量，所有的向量都被转换成沿着第一个主分量的点。当然，有一些信息的丢失，这必须逐案考虑；在接下来的几节中，我们将讨论如何评估这种损失并做出合理的决定。现在，我们将简要说明如何以有效的方式提取主成分。

奇异值分解的主成分分析

即使我们将采用完整的主成分分析实现，了解如何有效地执行这样的过程也是有帮助的。当然，最明显的方法是基于样本协方差矩阵的计算，它的特征分解(将输出特征值和相应的特征向量)，然后最终，有可能建立转换矩阵。这种方法很简单，但不幸的是，它的效率也很低。主要原因是我们需要计算样本协方差矩阵，这对于大型数据集来说可能是一项非常长的任务。

奇异值分解 ( SVD )提供了一种高效得多的方式，这是一个线性代数过程，具有一些重要的特性:它可以直接对数据集进行操作，当提取到所需数量的组件时可以停止，并且有可以小批量工作的增量版本，克服了内存不足的问题。特别是考虑到数据集 X ∈ ℜ ^{m × n} ，SVD 可以表示如下:

U 为酉矩阵(即UU^T= U^TU = I，故U^T= U^-1)包含作为行的左手奇异向量( XX ^T 的特征向量)； V (也是酉)包含右奇异向量作为行(对应于 X ^T X 的特征向量)，而λ是包含mσ_sT29】的奇异值的对角矩阵(它们是XX^TT33】和 X ^{T 的特征值的平方根特征值按降序排序，特征向量重新排列以匹配相应的位置。由于 1/m 因子是乘法常数，不影响特征值的相对大小；因此，排序顺序保持不变。因此，我们可以直接使用 V 或 U 并从}λ中选择第一个顶级 k 特征值。特别是我们可以观察到以下结果(作为变换矩阵， A ，等于 V ):**

因此，通过使用 U _k (仅包含顶部 k 特征向量)和λ_kT9】(仅包含顶部 k 特征值)的截断版本，我们可以直接获得低维变换数据集(具有 k 分量)，如下所示:

这种方法快速、有效，并且当数据集太大而无法放入内存时，可以很容易地进行缩放。即使我们在本书中没有处理这样的场景，提到 scikit-learn TruncatedSVD类(执行仅限于 k 顶特征值的 SVD)和IncrementalPCA类(对小批量执行 PCA)也是有帮助的。出于我们的目的，我们将使用标准的PCA类和一些重要的变体，它们要求整个数据集适合内存。

白粉

奇异值分解的一个重要应用是白化过程，该过程强制数据集 X 具有零均值(即 E[X] = 0 或零中心)，以具有单位协方差矩阵 C (它是实对称的)。这种方法非常有助于提高许多监督算法的性能，这些算法可以受益于所有组件共享的统一单一方差。

将分解应用于 C ，我们得到以下结果:

矩阵 V 的列是 C 的特征向量，而λ是包含特征值的对角矩阵(记住 SVD 输出奇异值，奇异值是特征向量的平方根)。因此，我们需要找到一个线性变换， z = Ax ，这样E【Z^TZ】= I。当使用前面的分解时，这很简单:

从前面的等式中，我们可以推导出变换矩阵的表达式 A :

现在，我们将使用一个小的测试数据集展示美白的效果，如下所示:

import numpy as np

from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=300, centers=1, cluster_std=2.5, random_state=1000)

print(np.cov(X.T))

显示数据集协方差矩阵的上一个块的输出如下:

[[6.37258226 0.40799363]
 [0.40799363 6.32083501]]

用于在通用数据集上执行白化的whiten()函数(零中心化是过程的一部分)如下所示(白化后correct参数强制进行比例校正):

import numpy as np

def zero_center(X):
    return X - np.mean(X, axis=0)

def whiten(X, correct=True):
    Xc = zero_center(X)
    _, L, V = np.linalg.svd(Xc)
    W = np.dot(V.T, np.diag(1.0 / L))
    return np.dot(Xc, W) * np.sqrt(X.shape[0]) if correct else 1.0

应用于X数组的白化结果如下截图所示:

Original dataset (left); whitened dataset (right)

我们现在可以检查新的协方差矩阵，如下所示:

import numpy as np

Xw = whiten(X)
print(np.cov(Xw.T))

输出如下:

[[1.00334448e+00 1.78229783e-17]
 [1.78229783e-17 1.00334448e+00]]

可以看到，矩阵现在是一个恒等式(误差最小)，数据集也有一个空平均值。

MNIST 数据集的主成分分析

现在，让我们应用主成分分析，以降低 MNIST 数据集的维度。我们将使用 scikit-learn 提供的压缩版本(1，797，8 × 8 图像)，但我们的任何考虑都不会受到这一选择的影响。让我们从加载和规范化数据集开始:

from sklearn.datasets import load_digits

digits = load_digits()
X = digits['data'] / np.max(digits['data'])

从理论讨论中，我们知道协方差矩阵特征值的大小与相应主成分的相对重要性(即解释的方差，因此也是信息量)成正比。因此，如果按降序排序，可以计算以下差异:

当组件数量 k → n 时，重要性趋于降低，我们可以通过选取第一个最大差异来选择最佳的 k ，这表明以下所有组件解释的差异量大幅下降。为了更好地理解这种机制，让我们计算特征值及其差(作为协方差矩阵， C ，是正半定的，我们确定λ_I≥0∀I∑(1，n) ):

import numpy as np

C = np.cov(X.T)
l, v = np.linalg.eig(C)
l = np.sort(l)[::-1]
d = l[:l.shape[0]-1] - l[1:]

展平图像(64 维阵列)的差异如下图所示:

Eigenvalue differences for each principal component

可以看出，第一主成分的差异非常大，对应于第四主成分(λ₄-λ₃T5)达到最大值；但是，接下来的差别还是很大的，而对应λ₆T9】则有一个急跌。在这一点上，趋势几乎是稳定的(除了一些残余振荡)，直到 λ ₁₁ ，然后开始非常迅速地下降，趋于零。由于我们还是要有正方形的图像，我们就要选择 k = 16 (相当于每边除以四)。在另一个任务中，你可以选择，例如 k = 15 ，甚至k = 8；然而，为了更好地理解降维引起的误差，分析解释的方差也将是有帮助的。因此，让我们从执行主成分分析开始:**

from sklearn.decomposition import PCA

pca = PCA(n_components=16, random_state=1000)
digits_pca = pca.fit_transform(X)

digits_pca阵列是在拟合模型并将所有样本投影到对应于前 16 个主成分的子空间之后获得的。如果我们想要将原始图像与其重建进行比较，我们需要调用inverse_transform()方法，该方法执行到原始空间的投影。所以，如果 PCA 在本例中是一个变换，f(x):ℜ⁶⁴→ℜ¹⁶T7】，逆变换是g(z):ℜ¹⁶→ℜ⁶⁴t13】。下面的截图显示了前 10 位数字及其重构之间的比较:**

Original samples (top row); reconstructions (bottom row)

重建很明显是有损耗的，但是数字仍然是可区分的。现在，让我们通过对explained_variance_ratio_数组的所有值求和来检查总的解释方差，该数组包含每个分量的解释方差的相对量(因此任何 k < n 分量的和总是小于 1):

print(np.sum(pca.explained_variance_ratio_))

上一个片段的输出如下:

0.8493974642542452

因此，在降维到 16 个分量的情况下，我们解释了大约 85%的原始方差，这是一个合理的值，考虑到我们为每个样本丢弃了 48 个分量。

显示所有单个贡献的图如下截图所示:

Explained variance ratio corresponding to each principal component

正如预期的那样，这种贡献往往会减少，因为在这种情况下，第一个主要组成部分是负责任的；例如，对于颜色，一条线(如黑色或白色)，而其余的线形成灰色阴影。这种行为非常普遍，几乎在每种情况下都能观察到。通过这个图表，也很容易找到额外的损失，以便进一步减少。例如，我们可以立即发现，对 3 个分量的严格限制将解释大约 40%的原始方差；所以，剩下的 45%被分成了剩下的 13 个部分。我邀请你重复这个例子，试图找到人类区分所有数字所需的最小组件数。

核主成分分析

有时，数据集不是线性可分的，标准的主成分分析不能提取正确的主成分。这个过程与第 3 章、高级聚类中讨论的过程没有什么不同，当时我们面对的是非凸聚类的问题。在这种情况下，由于几何形状的原因，一些算法无法成功分离。在这种情况下，目标是根据主成分的结构来区分不同的类(在纯无监督的场景中，我们考虑特定的分组)。因此，我们希望处理转换后的数据集， Z ，并检测可区分阈值的存在。例如，让我们考虑以下截图:

Original dataset (left); PCA projected version (right)

由于原始数据集是线性可分的，在 PCA 投影后，我们可以立即找到允许检测第一个分量的阈值(这是唯一真正需要的)，以便区分两个斑点。但是，如果数据集不是线性可分的，我们会得到一个不可接受的结果，如下图所示:

Original dataset (left); PCA projected version (right)

当几何形状更复杂时，找到可区分的阈值可能是不可能的。然而，我们知道将数据投影到更高维的空间可以使它们线性分离。特别是如果 x ∈ ℜ ⁿ ，我们可以选择一个合适的函数， f(x) ，这样y = f(x)∈ℜ^pT9】，用p>t23】n*。不幸的是，将这种转换应用于整个数据集可能非常昂贵；实际上，给定一个变换矩阵， A ，(带 n 个分量)，单个主分量，A^(t)，投影后，可以写成如下形式(记住它们是协方差矩阵的特征向量):*

因此，单个向量的变换如下:

可以看到，变换需要点积的计算，f(x_I)^Tf(x_I)。在这些情况下，我们可以使用所谓的内核技巧，它表示有特定的函数 K(，)，称为内核，具有有趣的属性，如下所示:

换句话说，我们可以通过简单地计算每两个点的核来计算在高维空间中主成分上的投影，而不是执行点积，这需要在计算*f()*之后进行 n 次乘法。

一些常见的内核如下:

径向基函数 ( 径向基函数，或高斯核:
多项式核，次数为 p :
乙状结肠核:

对于非常大的数据集，这个过程仍然相当昂贵(但是为了避免额外的时间，可以预先计算和存储内核值)，但是它比标准的投影要高效得多。此外，它的优点是允许在可能进行线性判别的空间中提取主成分。现在，让我们将径向基函数核主成分分析应用于上一张截图中显示的半月形数据集。gamma参数等于1/σ²T4。在这种特殊情况下，主要问题是存在双重重叠。考虑到原始标准差约为 1.0(也就是 σ ² = 1 ，我们至少需要三个标准差才能恰当区分；因此，我们将设置 γ = 10 :

from sklearn.datasets import make_moons
from sklearn.decomposition import KernelPCA

X, Y = make_moons(n_samples=800, noise=0.05, random_state=1000)

kpca = KernelPCA(n_components=2, kernel='rbf', gamma=10.0, random_state=1000)
X_pca = kpca.fit_transform(X)

投影结果如下图所示:

Original dataset (left); kernel PCA projected version (right)

可以看到，即使在这种情况下，第一个组件也足以做出决定(由于噪声，具有最小容差)，将阈值设置为零可以分离数据集。我邀请读者测试其他内核的效果并应用它们，以便区分包含全 0 和全 1 的 MNIST 子集。

通过因子分析增加异方差噪声的鲁棒性

标准主成分分析的主要问题之一是这种模型在异方差噪声方面的固有弱点。如果你不熟悉这个术语，引入两个定义会很有帮助。多元去相关噪声项的特征在于对角协方差矩阵 C ，其可以具有两种不同的配置，如下所示:

C = diag(σ ² ，σ ² ，...，σ ² ) :这种情况下，噪声定义为同态(所有分量方差相同)。
C = diag(σ₁T3】2，σ₂T7】2，...，σ_n²**)，同σ₁²≠σ₂²...≠ σ _n ² :这种情况下，噪声定义为异方差(每个分量都有自己的方差)。

有可能证明，当噪声是同态时，主成分分析可以很容易地管理它，因为单个成分的解释方差以同样的方式受到噪声项的影响(也就是说，这相当于没有噪声)。相反，当噪声为异方差时，主成分分析的性能下降，结果可能是绝对不可接受的。为此，Rubin 和 Thayer(在 ML 因子分析的 EM 算法中，Rubin D .和 Thayer D .，mentorimetrika，47，1982)提出了一种替代的降维方法，称为因子分析，可以解决这类问题。

假设我们有一个以零为中心的数据集, X ,包含 m 样本x_I∈ℜⁿT9】。我们的目标是找到一组潜在变量，z_I∈ℜ^p(带 p < n )和一个矩阵， A (称为因子加载矩阵)，这样每个样本都可以重写，如下所示:

因此，我们现在假设一个样本， x _i 是一组高斯潜变量加上一个额外的异方差噪声项的组合。由于潜在变量的维数较低，这个问题与标准主成分分析非常相似，主要区别在于，现在我们考虑了异方差噪声(当然，术语 n 也可以为空，或同方差)。因此，当确定分量(即潜在变量)时，不同噪声方差的影响被包括在模型中，最终效果是部分滤波(去噪)。在前面提到的论文中，作者提出了一种优化算法，它在形式上并不复杂，但需要许多数学运算(为此，我们省略了任何证明)。该方法基于期望最大化 ( EM )算法，便于找到最大化对数似然的参数集。在本书中，我们不讨论所有的数学细节(可以在原始论文中找到)，而是检查这种方法的属性，并将结果与标准主成分分析进行比较。

让我们从加载 Olivetti faces 数据集开始，将其置零，并创建一个异方差噪声版本，如下所示:

import numpy as np

from sklearn.datasets import fetch_olivetti_faces

faces = fetch_olivetti_faces(shuffle=True, random_state=1000)
X = faces['data']
Xz = X - np.mean(X, axis=0)

C = np.diag(np.random.uniform(0.0, 0.1, size=Xz.shape[1]))
Xnz = Xz + np.random.multivariate_normal(np.zeros(shape=Xz.shape[1]), C, size=Xz.shape[0])

下面的截图显示了一些原始且有噪声的图像:

Original images (upper line); noisy versions (lower line)

现在，让我们评估平均对数似然性(通过score()方法，在PCA和FactorAnalysis类中都可用)如下:

主成分分析，带有原始数据集和128组件
主成分分析，带噪声数据集和128组件
因子分析，含噪声数据集和128成分(潜在变量)

在下面的代码片段中，所有 3 个模型都被实例化和训练:

from sklearn.decomposition import PCA, FactorAnalysis

pca = PCA(n_components=128, random_state=1000)
pca.fit(Xz)
print('PCA log-likelihood(Xz): {}'.format(pca.score(Xz)))

pcan = PCA(n_components=128, random_state=1000)
pcan.fit(Xnz)
print('PCA log-likelihood(Xnz): {}'.format(pcan.score(Xnz)))

fa = FactorAnalysis(n_components=128, random_state=1000)
fa.fit(Xnz)
print('Factor Analysis log-likelihood(Xnz): {}'.format(fa.score(Xnz)))

上一个片段的输出如下:

PCA log-likelihood(Xz): 4657.3828125
PCA log-likelihood(Xnz): -2426.302304948351
Factor Analysis log-likelihood(Xnz): 1459.2912218162423

这些结果表明了在异方差噪声存在时因子分析的有效性。主成分分析获得的最大平均对数似然约为4657，在有噪声的情况下降至-2426。相反，因子分析获得的平均对数似然值约为 1，460，比使用主成分分析获得的对数似然值大得多(即使噪声的影响尚未完全滤除)。因此，每当数据集包含(或数据科学家怀疑它包含)异方差噪声时(例如，样本是通过不同仪器捕获的源的叠加获得的)，我强烈建议将因子分析作为主要的降维方法。当然，如果需要其他条件(例如，非线性、稀疏性等)，可以在做出最终决定之前评估本章中讨论的其他方法。

稀疏主成分分析和字典学习

标准 PCA 一般是密集分解；也就是说，向量一旦被变换，就是所有具有非零系数的分量的线性组合:

在前面的表达式中，系数 α _i 几乎总是不同于零，因此所有的组件都参与了重建过程。出于降维的目的，这不是一个问题，因为我们对每个组件解释的方差更感兴趣，以便限制它们。然而，对于一些任务，分析更大的一组建筑原子是有帮助的，假设每个向量可以表示为它们的稀疏组合。最经典的例子是文本语料库，其中词典包含的术语比每个文档中包含的术语多得多。这种模型通常被称为字典学习算法，因为原子集合定义了一种字典，包含了所有可以用来创建新样本的单词。当原子的数量 k 大于样本的维度 n 时，字典称为过完备，表示往往比较稀疏。反之，当 k < n 时，字典称为**欠完备，**向量需要更密集。

通过对解的 L ₁ 范数施加惩罚，这样的学习问题可以很容易地通过函数的最小化来解决。这种约束导致稀疏的原因超出了本书的范围，但感兴趣的人可以在掌握机器学习算法， 博纳科尔索 g .**Packt Publications2018 中找到更长的讨论。

字典学习(以及稀疏主成分分析)的问题可以正式表达如下:

这是一种算法的特殊情况，其中组件 U _k 被强制具有单位长度(除非有normalize_components=False参数)，系数 V 被惩罚，以便增加它们的稀疏度(与系数成比例， α )。

让我们考虑 MNIST 数据集，执行 30 个分量的稀疏主成分分析(产生不完全字典)和中高稀疏度水平(例如， α = 2.0 )。数组X应该包含归一化样本，如下面的主成分分析示例所示:

from sklearn.decomposition import SparsePCA

spca = SparsePCA(n_components=30, alpha=2.0, normalize_components=True, random_state=1000)
spca.fit(X)

在训练过程结束时，components_数组包含原子，如下图所示:

Components extracted by the sparse PCA algorithm

不难理解，每个数字都可以用这些原子组成；然而，考虑到原子的数量，稀疏的数量不可能非常大。我们来考虑一下，比如数字X[0]的变换:

y = spca.transform(X[0].reshape(1, -1)).squeeze()

系数的绝对值如下图所示:

Absolute coefficients for the sparse transformation of the digit X[0]

显然有一些主导成分(例如 2、7、13、17、21、24、26、27 和 30)，一些次要成分(例如 5、8 等)和一些无效或可忽略的成分(例如 1、3、6 等)。如果以相同的码长(30 个分量)增加稀疏度，则空分量对应的系数将降至零，而如果同样增加码长(例如 k = 100 ，则字典将变得过完备，空系数的数量也将增加。

非负矩阵分解

当数据集 X 为非负时，有可能应用因子分解技术，当任务的目标是提取与样本的结构部分相对应的原子时，该技术已被证明(例如在通过非负矩阵因子分解学习物体的部分，Lee D. D .，和 Seung，S. H .，Nature， 401，10/1999)更可靠。例如，在图像的情况下，它们应该是几何元素或者甚至更复杂的部分。非负矩阵分解 ( NNMF )强加的主要条件是所有涉及的矩阵必须是非负的并且 X = UV 。因此，一旦定义了一个规范 N (例如，弗罗贝尼乌斯)，简单的目标就变成了以下内容:

当还需要稀疏性时，这通常是不可接受的(此外，为了允许更大的灵活性来改变解决方案以满足特定的要求)，该问题通常通过在 Frobenius(对 L ₂ 的矩阵扩展)和 L ₁ 规范(例如，在 ElasticNet 中)上增加惩罚来表达(例如在 scikit-learn 中):

双重正则化通过避免类似于监督模型的过拟合的效果(由于解是次优的，它也更灵活地适应从相同的数据生成过程中提取的新样本)，允许您获得样本的稀疏性和部分之间的更好匹配；这增加了通常可实现的可能性)。

现在，让我们考虑 MNIST 数据集，让我们将其分解为 50 个原子，最初设置 α = 2.0 和 β = 0.1 (在 scikit-learn 中称为l1_ratio)。这种配置将强制中等稀疏度和强L2/弗罗贝纽斯正则化。该过程简单明了，类似于稀疏主成分分析:

from sklearn.decomposition import NMF

nmf = NMF(n_components=50, alpha=2.0, l1_ratio=0.1, random_state=1000)
nmf.fit(X)

在培训过程结束时，组件(原子)如下图所示:

Atoms extracted by the NNMF algorithm

与我们在标准词典学习中观察到的相反，原子现在的结构更加结构化，它们再现了数字的特定部分(例如，垂直或水平笔画、圆圈、点等)；因此，我们可以期待更多的稀疏表示，因为更少的组件足以构建一个数字。考虑到上一节所示的例子(数字X[0]，所有组件的绝对贡献如下图所示:

Absolute coefficients for the NNMF of the digit X[0]

三个成分占优势(3、24 和 45)；因此，我们可以尝试将样本表示为它们的组合。系数分别为 0.19、0.18 和 0.16。结果如下截图所示(数字X[0]表示零):

Deconstruction of the digit X[0], based on three main components

有趣的是该算法是如何选择原子的。即使这个过程受到 α 和 β 参数的强烈影响，通过范数，我们可以观察到，例如第三个原子(截图中的第一个)可以被许多零、三、八共享；最后一个原子对 0 和 9 都有帮助。每当原子的粒度太粗时，带有较弱的 L ₁ 惩罚的过完备字典可能会有所帮助。当然，每个问题都需要具体的解决方案；因此，我强烈建议用领域专家来检查原子的结构。作为练习，我邀请您将 NNMF 应用于另一个小图像数据集(例如，Olivetti、Cifar-10 或 STL-10)，并尝试找到隔离固定数量的结构部分所需的正确参数(例如，对于面部，它们可以是眼睛、鼻子和嘴巴)。

独立成分分析

当使用标准主成分分析(或其他技术，如因子分析)时，这些成分是不相关的，但不能保证它们在统计上是独立的。换句话说，假设我们有一个数据集， X ，从联合概率分布中得出，p(X)；如果有 n 组件，我们不能总是确定以下等式成立:

然而，有许多重要的任务，基于一个叫做鸡尾酒会的常见模式。在这种情况下，我们可以假设(或者我们知道)许多不同且独立的来源(例如，声音和音乐)重叠并产生单个信号。此时，我们的目标是通过对每个样本应用线性变换来尝试分离源。让我们考虑一个白化的数据集， X (因此所有的成分都具有相同的信息内容)，我们可以假设它是从高斯分布 N(0，I) 中采样的(这不是限制性条件，因为许多不同来源的重叠很容易收敛到正态分布)。因此，目标可以表达如下:

换句话说，我们将每个样本表示为多个独立因素的乘积，具有基于指数函数的先验分布。唯一必须绝对执行的条件是非高斯性(否则，组件变得不可区分)；因此，函数 f _k (z) 不能是二次多项式。在实践中，我们也希望包括适度的稀疏性，因此我们期望峰值和重尾分布(也就是说，概率仅在非常短的范围内很高，然后突然下降到几乎为零)。这种情况可以通过检查归一化的第四矩来验证，称为峰度:

对于高斯分布，峰度为 3。由于这通常是一个参考值，所有具有峰度(X) > 3 的分布称为超高斯或细峰度，而具有峰度(X) < 3 的分布称为亚高斯或扁峰度。前一类分布的一个例子是拉普拉斯分布，如下图所示:

Probability density functions of Gaussian distribution (left) and Laplace distribution (right)

不幸的是，峰度的使用因其相对于异常值缺乏鲁棒性而受到阻碍(也就是说，由于它涉及四次幂，即使很小的值也可以被放大并改变最终结果；例如，尾部有异常值的噪声高斯可以表现为超高斯)。为此，作者 Hyvarinen 和 Oja(在独立成分分析:算法与应用，Hyvarinen A .和 Oja，e .神经网络 13，2000)基于负熵的概念，提出了一种称为快速独立成分分析(FT6】astICAT8】)的算法。我们不打算在这本书里描述整个模型；然而，理解基本思想是有帮助的。可以证明，在所有方差相同的分布中，高斯分布的熵最大。因此，如果数据集 X (零中心)是从具有协方差的分布σ中提取的，则可以将 X 的负熵定义为高斯N(0；σ**)**X的熵:

因此，我们的目标是通过减少 H(X) 来减少 H _N (X) (总是大于等于零)。FastICA 算法基于通过特定函数的组合对 H _N (X) 的近似。最常见的称为 logcosh (也是 scikit-learn 中的默认值)，如下所示:

有了这个技巧，负熵可以更容易地优化，最终的分解肯定包含独立的成分。现在，让我们将 FastICA 算法应用于 MNIST 数据集(为了强制更好的精度，我们正在设置max_iter=10000和tol=1e-5):

from sklearn.decomposition import FastICA

ica = FastICA(n_components=50, max_iter=10000, tol=1e-5, random_state=1000)
ica.fit(X)

算法找到的 50 个独立组件(始终可通过components_实例变量获得)如下截图所示:

Independent components extracted by FastICA

在这种情况下，组件可以立即识别为数字的一部分(考虑到数据集的维度，我邀请读者通过将组件的数量减少和增加到 64 个来重复该示例，这是最大数量)。分量趋向于达到相应分布的平均位置；因此，用更少的数量，可以区分更结构化的模式(可以被认为是不同的重叠信号)，而更大数量的组件导致更多以特征为中心的元素。然而，与 NNMF 相反，FastICA 并不保证提取样本的实际部分，而是提取更完整的区域。换句话说，虽然 NNMF 可以很容易地检测到，例如，一些单笔画，但快速独立分量分析倾向于将样本视为不同信号的总和，在图像的情况下，这通常涉及样本的整个维度，除非分量的数量急剧增加。为了更好地理解这个概念，让我们考虑奥利维蒂人脸数据集，它包含 400 幅 64 × 64 灰度肖像:

from sklearn.datasets import fetch_olivetti_faces

faces = fetch_olivetti_faces(shuffle=True, random_state=1000)

下面的截图显示了前 10 张脸:

Sample faces extracted from the Olivetti faces dataset

现在，让我们提取 100 个独立的组件:

ica = FastICA(n_components=100, max_iter=10000, tol=1e-5, random_state=1000)
ica.fit(faces['data'])

前 50 个组件绘制在下面的截图中:

50 (out of 100) independent components extracted by FastICA

正如你所看到的，每个组件都类似于一个元面(有时称为特征面)，由于所有剩余的特征(即使它们在精确的样本集中不能立即识别)，它们具有一些特定的、主要的特征以及次要的贡献。当组件数量增加到 350 时，效果变得更加明显，如下所示:

50 (out of 350) independent components extracted by FastICA

在这种情况下，次要特征不太占优势，因为有更多的重叠分布，并且它们中的每一个都集中在更原子的特征面上。当然，没有完整的领域知识，无法定义组件的最佳数量。例如，在 Olivetti 人脸数据集的情况下，识别特定的子元素(例如，眼镜的位置)或更完整的面部表情可能会有所帮助。在前一种情况下，较大数量的组件产生更集中的解决方案(即使它们在全局上不太容易区分)，而在后一种情况下，较小数量的组件(如前一个示例)产生更完整的结果，在这种情况下可以评估不同的影响因素。就信号而言，分量的数量应该等于预期重叠因子的数量(假设它们是独立的)。例如，音频信号可以包含一个人在机场说话的录音，背景声音宣布航班。在这种情况下，场景可以由三个部分组成:两个声音和噪音。由于噪声将被部分分解为主要成分，最终数量将等于 2。

基于潜在狄利克雷分配的主题建模

我们现在将考虑另一种在处理文本文档时非常有用的分解(也就是 NLP)。理论部分不是很容易，因为需要很深的概率论和统计学习知识(可以在原论文潜伏狄利克雷分配，机器学习研究杂志， Blei D. ， Ng A. ，和 Jordan M. ，3，(2003)993-1022)；因此，我们只讨论主要元素，没有任何数学参考(更紧凑的描述也出现在机器学习算法第二版，Bonaccorso，g .， Packt Publications ，2018)。让我们考虑一组文本文档， d _j (称为语料库)，其原子(或成分)是单词， w _i :

收集所有单词后，我们可以建立一个字典:

我们还可以陈述以下不等式(*N()*计算一个集合的元素数):

这意味着文档之间的单词分布是稀疏的，因为单个文档中只使用了几个单词，而前者的选择是对称狄利克雷分布(模型以其命名)，这是极其稀疏的(而且，它是分类分布的共轭先验，是一个一次多项式分布，所以很容易纳入模型)。概率密度函数(由于分布是对称的，参数 α _i = α ∀ i )如下:

现在，让我们考虑将文档语义分组为主题， t _k ，并且假设每个主题都由少量特殊单词表征:

这意味着话题之间的词语分布也是稀疏的。所以我们有全联合概率(单词、话题)，我们要确定条件概率 p(w _i |t _k ) 和p(t_k| w_I)。换句话说，给定一个文档，它是术语的集合(每个术语都有一个边际概率 p(w _i ) ，我们要计算这样的文档属于特定主题的概率。由于文档被柔和地分配给所有的主题(也就是说，它可以不同程度地属于多个主题)，我们需要考虑一个稀疏的主题-文档分布，从中可以得出主题-混合( θ _i ):

类似地，我们需要考虑话题词分布(因为一个词可以不同程度地被更多话题分享)，从中我们可以得出话题词混合样本， β _j :

潜在狄利克雷分配 ( LDA )是一个生成模型(训练目标，以一种简单的方式，包括找到最优参数 α 和 γ )，该模型能够从语料库中提取固定数量的主题，并用一组词来表征它们。给定一个样本文档，它能够通过提供主题混合概率向量( θ _i = (p(t ₁ )、p(t ₂ )将其分配给一个主题，...，p(t_k))；它还可以处理看不见的文档(使用相同的字典)。

现在，让我们将 LDA 应用于 20 个新闻组数据集的一个子集，该子集包含数千条公开发布用于自然语言处理研究的消息。特别是，我们想要对rec.autos和comp.sys.mac.hardware子组进行建模。我们可以使用内置的 scikit-learn fetch_20newsgroups()功能，要求去掉所有不必要的页眉、页脚和引号(其他帖子附在答案上):

from sklearn.datasets import fetch_20newsgroups

news = fetch_20newsgroups(subset='all', categories=('rec.autos', 'comp.sys.mac.hardware'), remove=('headers', 'footers', 'quotes'), random_state=1000)

corpus = news['data']
labels = news['target']

此时，我们需要对语料库进行矢量化。换句话说，我们需要将每个文档转换成包含词汇表中每个单词的频率(计数)的稀疏向量:

我们将使用CountVectorizer类来执行这一步，要求去掉重音并删除相对频率很高但不具有代表性的终止词(例如，英语中的 and、the 等)。此外，我们正在强制标记器排除所有非纯文本的标记(通过设置token_pattern='[a-z]+')。在其他情况下，这种模式可能会有所不同，但在这种情况下，我们不想依赖数字和符号:

from sklearn.feature_extraction.text import CountVectorizer

cv = CountVectorizer(strip_accents='unicode', stop_words='english', analyzer='word', token_pattern='[a-z]+')
Xc = cv.fit_transform(corpus)

print(len(cv.vocabulary_))

上一个片段的输出如下:

因此，每个文档都是一个 14，182 维的稀疏向量(很明显，大多数值都是空的)。我们现在可以通过强加n_components=2来执行 LDA，因为我们期望提取两个主题:

from sklearn.decomposition import LatentDirichletAllocation

lda = LatentDirichletAllocation(n_components=2, learning_method='online', max_iter=100, random_state=1000)
Xl = lda.fit_transform(Xc)

训练过程结束后，components_实例变量包含每对(单词和话题)的相对频率(以计数表示)。因此，在我们的例子中，它的形状是(2，14，182)，而components_[i, j]元素，用 i ∈ (0，1) 和 j ∈ (0，14，181) ，可以解释为这个词的重要性， j ，为了定义这个话题， i 。因此，我们将有兴趣检查，例如，这两个主题的前 10 个单词:

import numpy as np

Mwts_lda = np.argsort(lda.components_, axis=1)[::-1]

for t in range(2):
    print('\nTopic ' + str(t))
    for i in range(10):
        print(cv.get_feature_names()[Mwts_lda[t, i]])

输出如下:

Topic 0
compresion
progress
deliberate
dependency
preemptive
wv
nmsu
bpp
coexist
logically

Topic 1
argues
compromising
overtorque
moly
forbid
cautioned
sauber
explosion
eventual
agressive

很容易(考虑到一些非常特殊的术语)理解Topic 0被分配给comp.sys.mac.hardware，另一个被分配给rec.autos(不幸的是，这个过程不能基于自动检测，因为语义必须由人类来解释)。为了评估模型，让我们考虑两个示例消息，如下所示:

print(corpus[100])
print(corpus[200])

输出限于几行，如下所示:

I'm trying to find some information on accelerator boards for the SE. Has
anyone used any in the past, especially those from Extreme Systems, Novy or
MacProducts? I'm looking for a board that will support extended video,
especially Radius's two-page monitor. Has anyone used Connectix Virtual in
conjunction with their board? Any software snafus? Are there any stats
anywhere on the speed difference between a board with an FPU and one
without? Please send mail directly to me. Thanks.

...

The new Cruisers DO NOT have independent suspension in the front.  They
still
run a straight axle, but with coils.  The 4Runner is the one with
independent
front.  The Cruisers have incredible wheel travel with this system. 

The 91-up Cruiser does have full time 4WD, but the center diff locks in
low range.  My brother has a 91 and is an incredibly sturdy vehicle which
has done all the 4+ trails in Moab without a tow.  The 93 and later is even
better with the bigger engine and locking diffs.

所以第一个帖子明显和图形有关，第二个是政治信息。让我们计算两者的主题混合，如下所示:

print(Xl[100])
print(Xl[200])

输出如下:

[0.98512538 0.01487462]
[0.01528335 0.98471665]

因此，第一条消息属于Topic 0的概率约为 98%，而第二条消息几乎不属于Topic 1。这证实了分解工作正常。为了更好地了解总体分布情况，将属于每个类别的消息的混合可视化将会很有帮助，如下图所示:

Topic-mixtures for comp.sys.mac.hardware (left) and rec.autos (right)

如你所见，主题几乎是正交的。属于rec.autos的大部分消息都有 p(t ₀ ) < 0.5 和 p(t ₁ ) > 0.5 ，而comp.sys.mac.hardware有轻微的重叠，其中没有p(t₀)>0.5和 p(t _{1】的消息组这可能是因为存在可以将两个主题限定为同等重要性的词语(例如，术语}讨论或辩论可以同等地出现在两个新闻组中)。作为练习，我邀请您使用更多子集，并尝试证明主题的正交性，并检测可能导致不正确作业的单词。

摘要

在这一章中，我们介绍了可用于降维和字典学习的不同技术。主成分分析是一种非常众所周知的方法，它涉及找到与方差较大的方向相关联的数据集的最重要的组成部分。这种方法具有对角化协方差矩阵和提供每个特征重要性的直接度量的双重效果，从而简化选择并最大化残差解释方差(可以用较少数量的分量解释的方差量)。由于主成分分析本质上是一种线性方法，它不能经常用于非线性数据集。为此，开发了一个基于内核的变体。在我们的例子中，您看到了径向基函数核如何能够将非线性可分数据集投影到子空间上，其中主成分分析可以确定判别分量。

稀疏主成分分析和字典学习是广泛使用的技术，当需要提取可以混合(线性组合)的建筑原子以产生样本时使用。在许多情况下，目标是找到一个所谓的过完备字典，这相当于说我们期望比实际用于构建每个样本的原子多得多的原子(这就是为什么表示是稀疏的)。虽然主成分分析可以提取不相关的成分，但它很少找不到统计独立的成分。为此，我们引入了独立分量分析的概念，这是一种为了从样本中提取重叠源而开发的技术，这些样本可以被认为是独立原因(例如，声音或视觉元素)的总和。另一种具有特殊特征的方法是 NNMF，它可以产生稀疏表示和一组类似于样本特定部分的成分(例如，对于一张脸，它们可以代表眼睛、鼻子等)。最后一节介绍了 LDA 的概念，LDA 是一种主题建模技术，在给定一个文档语料库(即一个文档属于每个特定主题的概率)的情况下，可以用来找到主题混合。

在下一章中，我们将介绍一些基于无监督范式的神经模型。具体来说，将讨论深度信念网络、自动编码器和无需协方差矩阵的特征分解(nor SVD)就能提取数据集主成分的模型。

问题

数据集 X 有一个协方差矩阵 C=diag(2，1) 。你对 PCA 有什么期待？
考虑到前面的问题，如果 X 以零为中心，球 B _0.5 (0，0) 为空，我们能假设 x = 0 (第一主成分)的阈值允许水平判别吗？
主成分分析提取的成分在统计上是独立的。这是正确的吗？
一个 Kurt(X) = 5 的分布适合 ICA。这是正确的吗？
包含样本( 1，2 )和( 0，-3 )的数据集 X 的 NNMF 是什么？
一个包含 10 个文档的语料库与一个包含 10 个术语的词典相关联。我们知道每份文件的固定长度是 30 个字。字典是不是太全了？
核主成分分析采用二次核。如果原始维数为 2，那么进行 PCA 的新空间的维数是多少？

进一步阅读

稀疏编码在线词典学习，J. Mairal ， F .巴赫， J. Ponce ，和 G. Sapiro ，2009
通过非负矩阵分解学习物体的各个部分，李德德，承宪，自然，401，10/1999
用于 ML 因子分析的 EM 算法、鲁宾 D. 、和泰尔 D. 、心理测量学卡，47，1982
独立成分分析:算法和应用，海瓦里宁和欧嘉，神经网络 13，2000
信息论的数学基础，钦钦人工智能，多佛出版物
潜在狄利克雷分配，机器学习研究杂志，布莱 D. ， Ng A. ，和乔丹 M. ，3，(2003) 993-1022
机器学习算法第二版，博纳科索格，帕克特出版，2018
掌握机器学习算法，博纳科索格，帕克特出版，2018

八、无监督神经网络模型

在本章中，我们将讨论一些可用于无监督任务的神经模型。神经网络(通常是深度神经网络)的选择允许您解决具有需要复杂处理单元(例如图像)的特定特征的高维数据集的复杂性。

特别是，我们将涵盖以下内容:

自动编码器
去噪自动编码器
稀疏自动编码器
可变自动编码器
主成分分析神经网络；
桑格网络
Rubner-Tavan 的网络
无监督的深度信念网络 ( DBN )

技术要求

本章中的代码要求如下:

Python 3.5+(强烈推荐 Anaconda 发行版
以下库:
- SciPy 0.19+
- NumPy 1.10+
- 学习 0.20+
- 熊猫 0.22+
- Matplotlib 2.0+
- seaborn 0.9+
- 张量流 1.5+
- 深度信仰网络(github.com/albertbup/d…)

这些例子可以在 GitHub 资源库中找到，网址为https://GitHub . com/packktpublishing/HandsOn-Unsupervised-Learning-with-Python/chapter 08。

自动编码器

在第 7 章、降维和分量分析中，我们讨论了一些常见的方法，这些方法可以用来降低数据集的维数，给定其特有的统计特性(例如，协方差矩阵)。然而，当复杂度增加时，即使核主成分分析 ( k ernel PCA )也可能无法找到合适的低维表示。换句话说，信息的丢失可以克服一个阈值，这个阈值保证了有效重建样本的可能性。自动编码器是利用神经网络的极端非线性的模型，以便找到给定数据集的低维表示。具体来说，我们假设 X 是从数据生成过程中抽取的一组样本， p _data (x) 。为简单起见，我们将考虑x_I∈ℜⁿ，但对支架的结构没有限制(例如，对于 RGB 图像，x_I∈ℜ^n×m×3)。自动编码器在形式上分为两个部分:一个是编码器，将高维输入转换为较短的代码，另一个是解码器，执行相反的操作(如下图所示):

Structural schema of a generic autoencoder

如果代码是一个 p 维向量，编码器可以定义为一个参数化函数，e():

以类似的方式，解码器是另一个参数化函数，d():

因此，一个完整的自动编码器是一个复合函数，给定一个输入样本， x _i ，它提供一个最优重构作为输出:

由于自动编码器通常通过神经网络实现，因此使用反向传播算法训练自动编码器，通常基于均方误差成本函数:

或者，考虑到数据生成过程，我们可以考虑参数化条件分布*q()*来重新表达目标:

因此，成本函数现在可以变成例如 *p _数据()和q()*之间的库尔巴克-莱布勒散度:

由于 p _数据 的熵为常数，可以通过优化过程排除；因此，散度的最小化相当于最小化 p _数据 和 q 之间的交叉熵。如果假设 p _数据 和 q 为高斯型，则库尔巴克-莱布勒成本函数相当于均方误差，因此两种方法可以互换。在某些情况下，当数据在(0，1)范围内归一化时，可以对 p _数据 和 q 采用伯努利分布。形式上，这并不完全正确，因为伯努利分布是二元的并且x_I∑{ 0，1 }^d；然而，sigmoid 输出单元的使用也保证了连续样本的成功优化，x_I∑(0，1) ^d 。在这种情况下，成本函数如下:

深度卷积自动编码器示例

让我们实现一个基于 TensorFlow 和 Olivetti faces 数据集的深度卷积自动编码器(该数据集相对较小，但提供了良好的表达能力)。让我们从加载图像和准备训练集开始:

from sklearn.datasets import fetch_olivetti_faces

faces = fetch_olivetti_faces(shuffle=True, random_state=1000)
X_train = faces['images']

这些样本是 400，64 × 64 灰度图像，我们将把它们的大小调整到 32 × 32，以加快计算速度并避免内存问题(此操作会导致视觉精度略有下降，如果您有足够的计算资源，可以将其删除)。我们现在可以定义主要常数(纪元数(nb_epochs)、batch_size、code_length)和graph:

import tensorflow as tf

nb_epochs = 600
batch_size = 50
code_length = 256 
width = 32
height = 32

graph = tf.Graph()

因此，我们将为 600 个时代训练模型，每批 50 个样本。由于每幅图像为 64 × 64 = 4，096 ，压缩比为 4，096/256 = 16 倍。当然，这种选择并不是一个规则，我邀请您始终检查不同的配置，以便最大限度地提高收敛速度和最终精度。在我们的案例中，我们使用以下层对编码器进行建模:

2D 卷积与 16 (3 × 3)滤波器，(2 × 2)步距，ReLU 激活和相同的填充
2D 卷积与 32 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充
2D 卷积与 64 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充
2D 卷积与 128 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充

解码器利用转置卷积序列(也称为去卷积):

2D 转置卷积具有 128 (3 × 3)个滤波器、(2 × 2)个步长、ReLU 激活和相同的填充
2D 转置卷积具有 64 (3 × 3)个滤波器，(1 × 1)步距，ReLU 激活和相同的填充
2D 转置卷积具有 32 (3 × 3)个滤波器、(1 × 1)个步长、ReLU 激活和相同的填充
具有 1 (3 × 3)滤波器、(1 × 1)步距、Sigmoid 激活和相同填充的 2D 转置卷积

损失函数基于重建图像和原始图像之间差异的 L ₂ 范数。优化器是亚当，学习率 η=0.001 。TensorFlow DAG 的编码器部分如下:

import tensorflow as tf

with graph.as_default():
    input_images_xl = tf.placeholder(tf.float32, 
                                     shape=(None, X_train.shape[1], X_train.shape[2], 1))
    input_images = tf.image.resize_images(input_images_xl, (width, height), 
                                          method=tf.image.ResizeMethod.BICUBIC)

    # Encoder
    conv_0 = tf.layers.conv2d(inputs=input_images,
                              filters=16,
                              kernel_size=(3, 3),
                              strides=(2, 2),
                              activation=tf.nn.relu,
                              padding='same')

    conv_1 = tf.layers.conv2d(inputs=conv_0,
                              filters=32,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

    conv_2 = tf.layers.conv2d(inputs=conv_1,
                              filters=64,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

    conv_3 = tf.layers.conv2d(inputs=conv_2,
                              filters=128,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

DAG 的代码部分如下:

import tensorflow as tf

with graph.as_default():   
    # Code layer
    code_input = tf.layers.flatten(inputs=conv_3)

    code_layer = tf.layers.dense(inputs=code_input,
                                 units=code_length,
                                 activation=tf.nn.sigmoid)

    code_mean = tf.reduce_mean(code_layer, axis=1)

DAG 的解码器部分如下:

import tensorflow as tf

with graph.as_default(): 
    # Decoder
    decoder_input = tf.reshape(code_layer, (-1, int(width / 2), int(height / 2), 1))

    convt_0 = tf.layers.conv2d_transpose(inputs=decoder_input,
                                         filters=128,
                                         kernel_size=(3, 3),
                                         strides=(2, 2),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_1 = tf.layers.conv2d_transpose(inputs=convt_0,
                                         filters=64,
                                         kernel_size=(3, 3),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_2 = tf.layers.conv2d_transpose(inputs=convt_1,
                                         filters=32,
                                         kernel_size=(3, 3),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_3 = tf.layers.conv2d_transpose(inputs=convt_2,
                                         filters=1,
                                         kernel_size=(3, 3),
                                         activation=tf.sigmoid,
                                         padding='same')

    output_images = tf.image.resize_images(convt_3, (X_train.shape[1], X_train.shape[2]), 
                                           method=tf.image.ResizeMethod.BICUBIC)

loss函数和 Adam 优化器在以下代码片段中定义:

import tensorflow as tf

with graph.as_default():
    # Loss
    loss = tf.nn.l2_loss(convt_3 - input_images)

    # Training step
    training_step = tf.train.AdamOptimizer(0.001).minimize(loss)

一旦定义了完整的 DAG，我们就可以初始化会话和所有变量:

import tensorflow as tf

session = tf.InteractiveSession(graph=graph)
tf.global_variables_initializer().run()

一旦初始化了 TensorFlow，就可以开始训练过程，如下所示:

import numpy as np

for e in range(nb_epochs):
    np.random.shuffle(X_train)

    total_loss = 0.0
    code_means = []

    for i in range(0, X_train.shape[0] - batch_size, batch_size):
        X = np.expand_dims(X_train[i:i + batch_size, :, :], axis=3).astype(np.float32)

        _, n_loss, c_mean = session.run([training_step, loss, code_mean],
                                        feed_dict={
                                            input_images_xl: X
                                        })
        total_loss += n_loss
        code_means.append(c_mean)

    print('Epoch {}) Average loss per sample: {} (Code mean: {})'.
          format(e + 1, total_loss / float(X_train.shape[0]), np.mean(code_means)))

上一个片段的输出如下:

Epoch 1) Average loss per sample: 11.933397521972656 (Code mean: 0.5420681238174438)
Epoch 2) Average loss per sample: 10.294102325439454 (Code mean: 0.4132006764411926)
Epoch 3) Average loss per sample: 9.917563934326171 (Code mean: 0.38105469942092896)
...
Epoch 600) Average loss per sample: 0.4635812330245972 (Code mean: 0.42368677258491516)

在训练过程结束时，每个样本的平均损失约为 0.46(考虑 32 × 32 个图像)，代码的平均值为 0.42。该值表示编码相当密集，因为单个值预期在范围(0，1)内均匀分布；因此，平均值为 0.5。在这种情况下，我们对这个数据不感兴趣，但是我们也将在寻找稀疏性时比较结果。

下图显示了一些示例图像的自动编码器输出:

Sample output of the deep convolutional autoencoder

放大到 64 × 64 部分影响重建质量；然而，通过降低压缩率和增加代码长度可以获得更好的结果。

去噪自动编码器

自动编码器的一个非常有用的应用并不严格与它们寻找低维表示的能力相关，而是依赖于从输入到输出的转换过程。特别是，让我们假设一个以零为中心的数据集， X ，以及一个噪声版本，其样本具有以下结构:

在这种情况下，自动编码器的目标是去除噪声项并恢复原始样本，x_IT3。从数学角度来看，标准和去噪自动编码器没有特别的区别；但是，考虑此类型号的容量需求非常重要。由于它们必须恢复原始样本，给定一个损坏的输入(其特征占据更大的样本空间)，图层的数量和尺寸可能比标准自动编码器更大。当然，考虑到复杂性，不经过几次测试是不可能有清晰的洞察力的；因此，我强烈建议从较小的模型开始，增加容量，直到最佳成本函数达到合适的值。对于噪声的添加，有几种可能的策略:

腐蚀每批样品(在整个时代)。
使用噪声层作为编码器的输入 1。
使用缺失层作为编码器的输入 1(例如，带有椒盐噪声)。在这种情况下，丢失的概率可以是固定的，或者可以在预定义的间隔内随机采样(例如，(0.1，0.5))。

如果假设噪声是高斯噪声(这是最常见的选择)，就有可能同时产生同异方差噪声和异方差噪声。在第一种情况下，所有分量的方差保持不变(即N(I)∞N(0，σ ² I) )，而在后一种情况下，每个分量都有自己的方差。根据问题的性质，另一种解决方案可能更合适；然而，当没有限制时，为了提高系统的整体鲁棒性，最好使用异方差噪声。

向深度卷积自动编码器添加噪声

在本例中，我们将修改之前开发的深度卷积自动编码器，以便管理有噪声的输入样本。DAG 几乎是等效的，不同的是，现在，我们需要同时提供噪声图像和原始图像:

import tensorflow as tf

with graph.as_default():
    input_images_xl = tf.placeholder(tf.float32, 
                                     shape=(None, X_train.shape[1], X_train.shape[2], 1))
    input_noisy_images_xl = tf.placeholder(tf.float32, 
                                           shape=(None, X_train.shape[1], X_train.shape[2], 1))

    input_images = tf.image.resize_images(input_images_xl, (width, height), 
                                          method=tf.image.ResizeMethod.BICUBIC)
    input_noisy_images = tf.image.resize_images(input_noisy_images_xl, (width, height), 
                                                method=tf.image.ResizeMethod.BICUBIC)

    # Encoder
    conv_0 = tf.layers.conv2d(inputs=input_noisy_images,
                              filters=16,
                              kernel_size=(3, 3),
                              strides=(2, 2),
                              activation=tf.nn.relu,
                              padding='same')
...

loss函数当然是通过考虑原始图像来计算的:

...

# Loss
loss = tf.nn.l2_loss(convt_3 - input_images)

# Training step
training_step = tf.train.AdamOptimizer(0.001).minimize(loss)

变量标准初始化后，我们可以开始训练过程，考虑一个加性噪声，N_I∞N(0，0.45) (即 σ ≈ 0.2 ):

import numpy as np

for e in range(nb_epochs):
    np.random.shuffle(X_train)

    total_loss = 0.0
    code_means = []

    for i in range(0, X_train.shape[0] - batch_size, batch_size):
        X = np.expand_dims(X_train[i:i + batch_size, :, :], axis=3).astype(np.float32)
        Xn = np.clip(X + np.random.normal(0.0, 0.2, size=(batch_size, X_train.shape[1], X_train.shape[2], 1)), 0.0, 1.0)

        _, n_loss, c_mean = session.run([training_step, loss, code_mean],
                                        feed_dict={
                                            input_images_xl: X,
                                            input_noisy_images_xl: Xn
                                        })
        total_loss += n_loss
        code_means.append(c_mean)

    print('Epoch {}) Average loss per sample: {} (Code mean: {})'.
          format(e + 1, total_loss / float(X_train.shape[0]), np.mean(code_means)))

一旦模型经过训练，就有可能用一些有噪声的样本来测试它。结果显示在下面的截图中:

Noisy samples (upper row); denoised images (lower row)

正如您所看到的，自动编码器已经成功地学会了如何对输入图像进行降噪，即使它们已经非常损坏。我邀请您用其他数据集测试模型，寻找允许合理良好重建的最大噪声方差。

稀疏自动编码器

标准自动编码器生成的代码通常很密集；然而，正如第 7 章、降维和成分分析中所讨论的，有时，最好使用过完备的字典和稀疏编码。实现这一目标的主要策略是简单地在成本函数中添加一个 L ₁ 惩罚(在代码层上):

α 常数决定了将要达到的稀疏程度。当然，由于C_sT5 的最优值与原值不对应，为了达到同样的精度，往往需要更多的纪元和更长的码层。吴恩达提出的另一种方法(在斯坦福大学的稀疏自动编码器、* CS294A 中)是基于一种稍微不同的方法。代码层被认为是一组独立的伯努利随机变量；因此，给定另一组具有小平均值的伯努利变量(例如，p_r∞B(0.05))，就有可能尝试找到最佳代码，该代码也使 z _i 和这样的参考分布之间的 kulback-Leibler 散度最小化:*

因此，新的成本函数如下:

最终的效果与使用L₁T3】点球取得的效果没有太大区别。事实上，在这两种情况下，模型都被迫学习次优表示，还试图最小化一个目标，如果单独考虑，该目标将导致输出代码始终为空。因此，全代价函数将达到保证重构能力和稀疏性的最小值(稀疏性必须总是与代码长度相平衡)。因此，一般来说，代码越长，可能实现的稀疏性就越大。

向深度卷积自动编码器添加稀疏性约束

在这个例子中，我们希望通过使用L₁T4】惩罚来增加代码的稀疏性。DAG 和训练过程与主例完全相同，唯一不同的是loss函数，现在变成如下:

...
sparsity_constraint = 0.01 * tf.reduce_sum(tf.norm(code_layer, ord=1, axis=1))
loss = tf.nn.l2_loss(convt_3 - input_images) + sparsity_constraint
...

我们添加了稀疏约束α= 0.01；因此，我们可以通过检查平均代码长度来重新训练模型。该过程的输出如下:

Epoch 1) Average loss per sample: 12.785746307373048 (Code mean: 0.30300647020339966)
Epoch 2) Average loss per sample: 10.576686706542969 (Code mean: 0.16661183536052704)
Epoch 3) Average loss per sample: 10.204148864746093 (Code mean: 0.15442773699760437)
...
Epoch 600) Average loss per sample: 0.8058895015716553 (Code mean: 0.028538944199681282)

如您所见，代码现在变得非常稀疏，最终平均值约为 0.03。这条信息表明，大多数代码值接近于零，在解码图像时只能考虑其中的少数几个。作为练习，我邀请您分析一组选定图像的代码，试图根据它们的激活/失活来理解它们的值的语义。

可变自动编码器

让我们考虑一个数据集， X ，来自一个数据生成过程， p _数据 。变分自动编码器是一种生成模型(基于标准自动编码器的主要概念)，由 Kingma 和 Welling(在*自动编码变分贝叶斯中，Kingma D. P .和 Welling m .*ArXiv:1312.6114【stat。ML] )，旨在再现数据生成过程。为了实现这个目标，我们需要从一个基于一组潜在变量 z 和一组可学习参数 θ 的通用模型开始。给定一个样本， x _i ∈ X ，模型的概率为 p(x，z； θ) 。因此，训练过程的目标是找到最大化可能性的最佳参数，p(x；θ) ，可通过全联合概率边缘化得到:

前面的表达式很简单，但不幸的是，它很少能以封闭的形式处理。主要原因是我们没有关于先验的有效信息，p(z； θ) 。而且，即使假设例如z∨N(0，σ)(例如 N(0，I) )，找到有效样本的概率也极其稀疏。换句话说，给定一个值 z ，我们也不太可能生成一个实际上属于 p _数据 的样本。为了解决这个问题，作者提出了一种变分方法，我们将简要地看一下(在前面提到的论文中给出了完整的解释)。假设标准自动编码器的结构，我们可以通过将编码器建模为q(z | x；θ _q ) 。此时，我们可以计算q()和实际条件概率p(z | x；θ) :

当期望值运算符在 z 上工作时，可以提取最后一个术语并将其移动到表达式的左侧，如下所示:

经过另一个简单的操作，前面的等式变成如下:

左边是模型下样本的对数似然，右边是一个非负项(KL-散度)和另一个名为证据下界 ( ELBO )的项之和:

正如我们将要讨论的，使用 ELBO 比处理公式的剩余部分更容易，并且，由于 KL-散度不能引入负贡献，如果我们最大化 ELBO，我们也最大化对数似然。

在上式中， p 是码长，所以是均值和对角协方差向量的维数。右边的表达式非常容易计算，因为σ是对角的(即迹是元素的和，行列式是积)。然而，当采用随机梯度下降 ( SGD )算法时，该公式的最大化虽然正确，但不是可微的操作。为了克服这个问题，作者建议重新参数化分布。

当出现批次时，对正态分布进行采样，得到α∨N(0，I) 。使用该值，可以通过使用概率编码器的输出来构建所需的样本:μ(z | x；θ_q)+ασ(z | x；θ_q**)²。这个表达式是可微的，因为 α 在每批期间是恒定的(当然，作为μ(z | x；θ _q ) 和σ(z | x；θ _q ) 用神经网络参数化，它们是可微的)。

ELBO 右侧的第二项是 log p(x|z)的期望值； θ) 。很容易看出，这样的表达式对应于原始分布和重构之间的交叉熵:

这是标准自动编码器的成本函数，假设使用伯努利分布，我们将把它最小化。所以，公式变成如下:

深度卷积变分自动编码器示例

在这个例子中，我们想要建立和训练一个基于 Olivetti 人脸数据集的深度卷积变分自动编码器。该结构与我们第一个示例中使用的结构非常相似。编码器有以下几层:

2D 卷积与 16 (3 × 3)滤波器，(2 × 2)步距，ReLU 激活和相同的填充
2D 卷积与 32 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充
2D 卷积与 64 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充
2D 卷积与 128 (3 × 3)滤波器、(1 × 1)步距、ReLU 激活和相同的填充

解码器具有以下转置卷积:

2D 转置卷积具有 128 (3 × 3)个滤波器、(2 × 2)个步长、ReLU 激活和相同的填充
2D 转置卷积具有 128 (3 × 3)个滤波器、(2 × 2)个步长、ReLU 激活和相同的填充
2D 转置卷积具有 32 (3 × 3)个滤波器、(1 × 1)个步长、ReLU 激活和相同的填充
具有 1 (3 × 3)滤波器、(1 × 1)步距、Sigmoid 激活和相同填充的 2D 转置卷积

噪音的产生完全由 TensorFlow 管理，它基于理论部分解释的技巧。包含图形定义和编码器的 DAG 的第一部分显示在下面的代码片段中:

import tensorflow as tf

nb_epochs = 800
batch_size = 100
code_length = 512
width = 32
height = 32

graph = tf.Graph()

with graph.as_default():
    input_images_xl = tf.placeholder(tf.float32, 
                                     shape=(batch_size, X_train.shape[1], X_train.shape[2], 1))
    input_images = tf.image.resize_images(input_images_xl, (width, height), 
                                          method=tf.image.ResizeMethod.BICUBIC)

    # Encoder
    conv_0 = tf.layers.conv2d(inputs=input_images,
                              filters=16,
                              kernel_size=(3, 3),
                              strides=(2, 2),
                              activation=tf.nn.relu,
                              padding='same')

    conv_1 = tf.layers.conv2d(inputs=conv_0,
                              filters=32,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

    conv_2 = tf.layers.conv2d(inputs=conv_1,
                              filters=64,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

    conv_3 = tf.layers.conv2d(inputs=conv_2,
                              filters=128,
                              kernel_size=(3, 3),
                              activation=tf.nn.relu,
                              padding='same')

DAG 中定义代码层的部分如下:

import tensorflow as tf

with graph.as_default():
    # Code layer
    code_input = tf.layers.flatten(inputs=conv_3)

    code_mean = tf.layers.dense(inputs=code_input,
                                units=width * height)

    code_log_variance = tf.layers.dense(inputs=code_input,
                                        units=width * height)

    code_std = tf.sqrt(tf.exp(code_log_variance))

DAG 的解码器部分如下:

import tensorflow as tf

with graph.as_default():  
    # Decoder
    decoder_input = tf.reshape(sampled_code, (-1, int(width / 4), int(height / 4), 16))

    convt_0 = tf.layers.conv2d_transpose(inputs=decoder_input,
                                         filters=128,
                                         kernel_size=(3, 3),
                                         strides=(2, 2),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_1 = tf.layers.conv2d_transpose(inputs=convt_0,
                                         filters=128,
                                         kernel_size=(3, 3),
                                         strides=(2, 2),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_2 = tf.layers.conv2d_transpose(inputs=convt_1,
                                         filters=32,
                                         kernel_size=(3, 3),
                                         activation=tf.nn.relu,
                                         padding='same')

    convt_3 = tf.layers.conv2d_transpose(inputs=convt_2,
                                         filters=1,
                                         kernel_size=(3, 3),
                                         padding='same')

    convt_output = tf.nn.sigmoid(convt_3)

    output_images = tf.image.resize_images(convt_output, (X_train.shape[1], X_train.shape[2]), 
                                           method=tf.image.ResizeMethod.BICUBIC)

DAG 的最后一部分包含损失函数和 Adam 优化器，如下所示:

import tensorflow as tf

with graph.as_default():
    # Loss
    reconstruction = tf.nn.sigmoid_cross_entropy_with_logits(logits=convt_3, labels=input_images)
    kl_divergence = 0.5 * tf.reduce_sum(
            tf.square(code_mean) + tf.square(code_std) - tf.log(1e-8 + tf.square(code_std)) - 1, axis=1)

    loss = tf.reduce_sum(tf.reduce_sum(reconstruction) + kl_divergence)

    # Training step
    training_step = tf.train.AdamOptimizer(0.001).minimize(loss)

损耗函数由两部分组成:

基于交叉熵的重构损失
码分布与参考正态分布之间的库尔巴克-莱布勒散度

此时，像往常一样，我们可以初始化会话和所有变量，并开始 800 个时期和每批 100 个样本的训练过程:

import tensorflow as tf
import numpy as np

session = tf.InteractiveSession(graph=graph)
tf.global_variables_initializer().run()

for e in range(nb_epochs):
    np.random.shuffle(X_train)

    total_loss = 0.0

    for i in range(0, X_train.shape[0] - batch_size, batch_size):
        X = np.zeros((batch_size, 64, 64, 1), dtype=np.float32)
        X[:, :, :, 0] = X_train[i:i + batch_size, :, :]

        _, n_loss = session.run([training_step, loss],
                                feed_dict={
                                    input_images_xl: X
                                })
        total_loss += n_loss

    print('Epoch {}) Average loss per sample: {}'.format(e + 1, total_loss / float(batch_size)))

在训练过程的最后，我们可以测试几个样本的重建。结果如下图所示:

Sample reconstructions produced by the variational autoencoder

作为练习，我邀请读者修改 DAG，以便接受通用输入代码并评估模型的生成属性。或者，可以获取训练样本的代码并应用一些噪声，以便观察对输出重建的影响。

基于 Hebbian 的主成分分析

在本节中，我们将分析两个神经模型(桑格和鲁布纳-塔万网络)，它们可以执行主成分分析 ( 主成分分析)而不需要特征分解协方差矩阵或执行截断的奇异值分解。它们都是基于赫布边学习的概念(更多细节请参考达扬，p .和雅培，L. F .，理论神经科学，麻省理工学院出版社， 2005 或博纳科索，g .，掌握机器学习算法， Packt，2018)，这是最早关于非常简单神经元动力学的数学理论之一。然而，这样的概念有非常有趣的含义，特别是在组件分析领域。为了更好地理解网络的动力学，快速概述神经元的基本模型将是有帮助的。让我们考虑一个输入， x ∈ ℜ ⁿ ，以及一个权重向量， w ∈ ℜ ⁿ 。神经元执行点积(无偏差)，以产生标量输出， y :

现在，如果我们想象两个神经元，第一个叫做突触前单位，另一个叫做突触后单位。 Hebb 法则规定突触强度在突触前和突触后单位输出相同符号的值(特别是两者都为正)时必须增加，而在符号不同时必须减弱。这样一个概念的数学表达如下:

常数 η 是学习速率。完整的分析超出了本书的范围，但是有可能证明一个赫伯神经元(通过一些非常简单的修改，需要控制 w 的生长)修改了突触权重，使得在足够多的迭代之后，它沿着数据集的第一个主成分 X 对齐。从这个结果出发(我们就不证明了)，我们可以引入桑格的网络。

桑格网络

桑格的网络模型是由桑格提出的(在桑格，T. D .，单层线性前馈神经网络中的最优无监督学习，神经网络， 2，1989)，目的是通过在线过程(相反，标准 PCA 是需要整个数据集的批处理)以降序提取数据集的第一个 k 主成分， X 。即使有一个基于特定版本的奇异值分解的增量算法，这些神经模型的主要优势是它们处理单个样本而不损失性能的内在能力。在展示网络结构之前，有必要先介绍一下对 Hebb 规则的修改，称为 Oja 规则:

引入这一规则是为了解决标准 Hebbian 神经元无限增长的问题。其实很容易理解，如果点积 w ^T x 为正，δw 会越来越多地通过增加 w 的量级来更新权重。因此，在大量迭代之后，模型可能会遇到溢出。Oja 规则通过引入一个自动限制来克服这个问题，该限制迫使权重大小饱和，而不影响神经元找到第一个主成分方向的能力。实际上，用 w _k 表示 k ^次迭代后的权重向量，可以证明如下:

桑格的网络是基于欧嘉规则的修改版本，定义为广义赫布边学习 ( GHL )。假设我们有一个数据集， *X，*包含 m 向量，x_I∈ℜⁿ。网络结构如下图所示:

Structure of a generic Sanger's network

权重组织成矩阵，W = { W_ij}(W_ijT7】是连接突触前单位 i 与突触后单位 j 的权重)；因此，输出的激活可以通过使用以下公式来计算:

然而，在这种网络中，我们更感兴趣的是最终的权重，因为它们必须等于第一个 n 个主成分。不幸的是，如果我们不加任何修改地应用 Oja 规则，所有的神经元都会找到相同的成分(第一个)；因此，必须采用不同的策略。从理论上，我们知道主成分必须是正交的；因此，如果w₁T5】是第一个分量方向的向量，我们可以强制w₂T9】与w₁T13】正交，依此类推。该方法基于克-施密特正交化程序。让我们考虑两个向量——已经收敛的 w ₁ 和w₂₀，在没有任何干预的情况下，也将收敛到 w ₁ 。通过考虑这个向量在 w ₁ 上的投影，我们可以找到 w ₂₀ 的正交分量:**

此时，正交分量 w ₂ 等于以下值:

第三个分量必须正交于 w ₁ 和 w ₂ ，因此必须对所有 n 单元重复该过程，直到最终收敛。此外，我们现在使用的是已经融合的组件，而是一个并行更新的动态系统；因此，有必要将此过程纳入学习规则，如下所示:

之前的更新是指单个权重， w _ij ，给定一个输入， x 。很容易理解，第一部分是标准的赫布规则，而剩下的部分是正交项，它扩展到 y _i 之前的所有单位。

在矩阵形式中，更新如下:

*Tril()*函数计算正方形矩阵的下三角形部分。收敛性证明不是无足轻重的，但是在 η 单调递减的温和条件下，可以看到这个模型是如何按降序收敛到第一个 n 主成分的:

这样的约束不难实现；但是，一般来说，算法在 η < 1 时也能达到收敛，并且在迭代过程中保持不变。

桑格网络的一个例子

让我们考虑一个样本二维零中心数据集，它是用 scikit-learn make_blobs()效用函数获得的:

import numpy as np

def zero_center(Xd):
    return Xd - np.mean(Xd, axis=0)

X, _ = make_blobs(n_samples=500, centers=3, cluster_std=[5.0, 1.0, 2.5], random_state=1000)
Xs = zero_center(X)

Q = np.cov(Xs.T)
eigu, eigv = np.linalg.eig(Q)

print('Covariance matrix: {}'.format(Q))
print('Eigenvalues: {}'.format(eigu))
print('Eigenvectors: {}'.format(eigv.T))

上一个片段的输出如下:

Covariance matrix: [[18.14296606  8.15571356]
 [ 8.15571356 22.87011239]]
Eigenvalues: [12.01524122 28.99783723]
Eigenvectors: [[-0.79948496  0.60068611]
 [-0.60068611 -0.79948496]]

特征值分别约为 12 和 29，说明第一主成分(对应转置特征向量矩阵的第一行，所以 (-0.799，0.6) )比第二主成分短很多。当然，在这种情况下，我们已经通过特征分解协方差矩阵计算了主成分，但这只是为了教学目的。桑格网络将按降序提取组件；因此，我们期望找到第二列作为权重矩阵的第一列，第一列作为第二列。让我们从初始化权重和训练常数开始:

import numpy as np

n_components = 2
learning_rate = 0.01
nb_iterations = 5000
t = 0.0

W_sanger = np.random.normal(scale=0.5, size=(n_components, Xs.shape[1]))
W_sanger /= np.linalg.norm(W_sanger, axis=1).reshape((n_components, 1))

In order to reproduce the example, it's necessary to set the random seed equal to 1,000; that is, np.random.seed(1000).

在这种情况下，我们执行固定次数的迭代(5000 次)；但是，我邀请您修改该示例，以便使用基于两个后续时间步骤计算的权重之差的范数(例如，Frobenius)的容差和停止标准(该方法可以通过避免无用的迭代来加快训练速度)。

下图显示了初始配置:

Initial configuration of the Sanger's network

此时，我们可以开始训练循环，如下所示:

import numpy as np

for i in range(nb_iterations):
    dw = np.zeros((n_components, Xs.shape[1]))
    t += 1.0

    for j in range(Xs.shape[0]):
        Ysj = np.dot(W_sanger, Xs[j]).reshape((n_components, 1))
        QYd = np.tril(np.dot(Ysj, Ysj.T))
        dw += np.dot(Ysj, Xs[j].reshape((1, X.shape[1]))) - np.dot(QYd, W_sanger)

    W_sanger += (learning_rate / t) * dw
    W_sanger /= np.linalg.norm(W_sanger, axis=1).reshape((n_components, 1))

print('Final weights: {}'.format(W_sanger))
print('Final covariance matrix: {}'.format(np.cov(np.dot(Xs, W_sanger.T).T)))

上一个片段的输出如下:

Final weights: [[-0.60068611 -0.79948496]
 [-0.79948496  0.60068611]]
Final covariance matrix: [[ 2.89978372e+01 -2.31873305e-13]  
 [-2.31873305e-13 1.20152412e+01]]

如您所见，最终协方差矩阵如预期的那样去相关，权重已经收敛到 C 的特征向量。重量(主要部件)的最终配置如下图所示:

Final configuration of the Sanger's network

第一主成分对应权重， w ₀ ，最大， w ₁ 为第二主成分。我邀请您使用高维数据集测试网络，并基于协方差矩阵的奇异值分解或特征分解，将性能与标准算法进行比较。

Rubner-Tavan 的网络

另一个可以进行主成分分析的神经网络是由 Rubner 和 Tavan 提出的(在 Rubner，j .和 Tavan，p .，一个用于主成分分析的自组织网络，欧洲物理学快报， 10(7)，1989)。然而，他们的方法是基于协方差矩阵的去相关，协方差矩阵是主成分分析的最终结果(也就是说，它就像用自下而上的策略操作，而标准过程是自上而下的)。让我们考虑一个以零为中心的数据集， X ，以及一个输出为 y ∈ ℜ ^m 向量的网络。因此，输出分布的协方差矩阵如下:

Structure of a generic Rubner-Tavan's network

如您所见，与桑格网络的主要区别在于每个输出单元之前都有求和节点(第一个除外)。这种方法被称为分层横向连接，因为每个节点 y _i (i > 0) 都是由一个直接分量 n _i 组成的，加起来就是所有先前的加权输出。因此，假设符号 v ⁽ⁱ⁾ ，表示向量的 i ^th 分量， i ，网络输出如下:

已经证明，这个模型，有一个特定的权重更新规则(我们将讨论)，收敛到一个单一的，稳定的，不动点，输出被迫成为相互去相关。从模型的结构来看，操作顺序如下:

第一个输出保持不变
第二输出被迫与第一输出去相关
第三输出被迫与第一和第二输出去相关，以此类推
最后一个输出被迫与所有先前的输出去相关

经过多次迭代，每一次生产，y_Iy_jT5【同】 i ≠ j 变为空， C 变为对角协方差矩阵。此外，在上述论文中，作者证明了特征值(对应于方差)是按降序排序的；因此，可以通过取包含第一个 p 行和列的子矩阵来选择顶部的 p 组件。

Rubner-Tavan 网络通过使用两个不同的规则进行更新，每个权重层一个规则。内部权重 w _ij 通过使用 Oja 规则进行更新:

该规则保证了主成分的提取不会出现w_ijT3】的无限增长。相反，外部权重 v _jk 通过使用反希伯来人规则来更新:

前一个公式的第一项*-ηy^(j)y^(k)负责去相关，而第二项类似于 Oja 规则，充当防止权重溢出的自限制正则化器。特别是-ηy^(I)y^(k)项可以解释为更新规则的反馈信号，δw_ij，受δv_JK*项修正的实际输出影响。考虑到桑格网络的行为，不容易理解的是，一旦输出被去相关，内部权重 w _ij 就变成正交的，代表 X 的第一主成分。

在矩阵形式中，权重 w _ij 可以立即排列成 W = {w _ij } ，这样在训练过程结束时，每一列都是 C 的特征向量(降序)。相反，对于外部砝码 v _jk ，我们需要再次使用*Tril()*运算符:

因此，迭代 t+1 时的输出变为如下:

有趣的是，这样的网络有一个循环输出；因此，一旦应用了输入，需要几次迭代才能使 y 稳定下来(理想情况下，更新必须持续到*| | y^(t+1)-y^(t)| |→0*)。

鲁布纳-塔万网络的一个例子

在本例中，我们将使用桑格网络示例中定义的数据集，以便使用鲁布纳-塔万网络执行主成分提取。为了方便起见，让我们重新计算特征分解:

import numpy as np

Q = np.cov(Xs.T)
eigu, eigv = np.linalg.eig(Q)

print('Eigenvalues: {}'.format(eigu))
print('Eigenvectors: {}'.format(eigv.T))

上一个片段的输出如下:

Eigenvalues: [12.01524122 28.99783723]
Eigenvectors: [[-0.79948496 0.60068611]
 [-0.60068611 -0.79948496]]

我们现在可以初始化超参数，如下所示:

n_components = 2
learning_rate = 0.0001
max_iterations = 1000
stabilization_cycles = 5
threshold = 0.00001

W = np.random.normal(0.0, 0.5, size=(Xs.shape[1], n_components))
V = np.tril(np.random.normal(0.0, 0.01, size=(n_components, n_components)))
np.fill_diagonal(V, 0.0)

prev_W = np.zeros((Xs.shape[1], n_components))
t = 0

因此，我们选择使用等于 0.00001 的停止阈值(比较基于权重矩阵的两次连续计算的弗罗贝纽斯范数)和最大 1000 次迭代。我们还设置了五个稳定周期和一个固定的学习速率， η=0.0001 。我们可以开始学习过程，如下所示:

import numpy as np

while np.linalg.norm(W - prev_W, ord='fro') > threshold and t < max_iterations:
    prev_W = W.copy()
    t += 1

    for i in range(Xs.shape[0]):
        y_p = np.zeros((n_components, 1))
        xi = np.expand_dims(Xs[i], 1)
        y = None

        for _ in range(stabilization_cycles):
            y = np.dot(W.T, xi) + np.dot(V, y_p)
            y_p = y.copy()

        dW = np.zeros((Xs.shape[1], n_components))
        dV = np.zeros((n_components, n_components))

        for t in range(n_components):
            y2 = np.power(y[t], 2)
            dW[:, t] = np.squeeze((y[t] * xi) + (y2 * np.expand_dims(W[:, t], 1)))
            dV[t, :] = -np.squeeze((y[t] * y) + (y2 * np.expand_dims(V[t, :], 1)))

        W += (learning_rate * dW)
        V += (learning_rate * dV)

        V = np.tril(V)
        np.fill_diagonal(V, 0.0)

        W /= np.linalg.norm(W, axis=0).reshape((1, n_components))

print('Final weights: {}'.format(W))

上一个块的输出如下:

Final weights: [[-0.60814345 -0.80365858]
 [-0.79382715 0.59509065]]

正如预期的那样，权重收敛到协方差矩阵的特征向量。让我们也计算最终的协方差矩阵，以便检查它的值:

import numpy as np

Y_comp = np.zeros((Xs.shape[0], n_components))

for i in range(Xs.shape[0]):
    y_p = np.zeros((n_components, 1))
    xi = np.expand_dims(Xs[i], 1)

    for _ in range(stabilization_cycles):
        Y_comp[i] = np.squeeze(np.dot(W.T, xi) + np.dot(V.T, y_p))
        y_p = y.copy()

print('Final covariance matrix: {}'.format(np.cov(Y_comp.T)))

输出如下:

Final covariance matrix: [[28.9963492 0.31487817]
 [ 0.31487817 12.01606874]]

最后的协方差矩阵也是去相关的(误差可以忽略不计)。Rubner-Tavan 的网络通常比 Sanger 的网络更快，因为反 Hebbian 反馈加快了收敛速度；因此，当采用这种模型时，它们应该是首选。然而，为了避免振荡，调整学习速率非常重要。我建议从一个小值开始，稍微增加它，直到迭代次数达到最小值。或者，可以从更高的学习速率开始，允许更快的初始校正，并通过使用线性(如桑格网络)或指数衰减来逐步降低它。

无监督深度信念网络

在本节中，我们将讨论一个非常著名的生成模型，在无监督的情况下，可以使用该模型对从预定义的数据生成过程中提取的输入数据集 X 进行降维。由于这本书没有特别的先决条件，而且数学复杂性相当高，我们将简要介绍概念，没有证明，也没有对算法的结构进行深入分析。在讨论深度信念网络 ( DBN )之前，有必要介绍另一个模型，受限玻尔兹曼机 ( RBM )，它可以被认为是一个 DBN 的积木。

受限玻尔兹曼机器

这个网络，也被称为调和，在动力系统中的信息处理:调和理论的基础 ，并行分布式处理，第 1 卷，麻省理工学院出版社， 1986)中被提出作为一个概率生成模型。换句话说，RBM 的目标是学习未知的分布(即数据生成过程)，以便生成所有可能的样本。通用结构如下图所示:

Structure of a generic Restricted Boltzmann Machine

神经元 x _i 是可观察的(也就是说，它们代表 RBM 必须学习的过程产生的向量)，而h_jT7】是潜在的(也就是说，它们是隐藏的，对 x _i 假设的值有贡献)。没有任何进一步的细节，我们需要说这个模型具有一个马尔可夫随机场 ( MRF )的结构，这得益于同一层的神经元之间没有连接(即描述网络的图是二分图)。MRFs 的一个重要性质是对全联合概率建模的可能性， p(x，h；θ) ，吉布斯分布:

指数 *E(x，h，**【θ】*起着物理系统能量的作用，在我们的例子中，它等于:

这个公式的主要假设是所有的神经元都是伯努利分布的(即 x _i ，h_j∞B(0，1) )，术语 b _i 和 c _j 是对可观察和潜在单位的偏差。给定一个数据生成过程， p _数据 ，一个 RBM 必须被优化，使得可能性，p(x；**【θ)，被最大化。跳过所有的中间步骤(可以在前面提到的论文中找到)，可以证明以下几点:

在前面的公式中，*σ()*是 sigmoid 函数。给定这两个表达式，就有可能导出(操作被省略)对数似然相对于所有可学习变量的梯度:

很容易理解，所有梯度的第一项非常容易计算，而所有第二项需要所有可能的可观察值的总和。这显然是一个无法用封闭形式解决的棘手问题。为此，Hinton(在 Hinton，G. *《训练受限玻尔兹曼机器的实用指南》，*多伦多大学计算机科学系，2010)提出了一种称为对比发散的算法，可用于寻找近似解。这种方法的解释需要马氏链的知识(这不是先决条件)；然而，我们可以总结这种策略，说它通过有限(和少量)的采样步骤来计算梯度的近似值(通常，一个步骤就足以获得好的结果)。这种方法可以非常有效地训练 RBM，并使深度信念网络易于使用且极其有效。

深层信念网络

DBN 是一个基于成果管理制的堆叠模型。通用结构如下图所示:

Structure of a generic DBN

第一层包含可见单元，而其余所有单元都是潜在的。在无监督的情况下，目标是学习未知的分布，找出样本的内部表示。事实上，当潜在单元的数量小于输入单元的数量时，该模型学习如何利用低维子空间对分布进行编码。辛顿和奥森德罗(在辛顿和奥森德罗的《深度信念网的快速学习算法》，神经计算， 18/7，2005)提出了一个逐步贪婪训练过程(通常实现的)。每对层被认为是 RBM，并使用对比发散算法进行训练。一旦训练了一个 RBM，隐藏的层就变成了后续 RBM 的可观察层，这个过程一直持续到最后一层。因此，DBN 开发了一系列内部表示(这就是为什么它被定义为深度网络)，其中每个层次都是在较低层次的特征上训练的。这个过程与变化的自动编码器没有什么不同；然而，在这种情况下，模型的结构更加严格(例如，不可能使用卷积单元)。此外，输出不是输入的重构，而是内部表示。因此，考虑到上一节中讨论的公式，如果有必要反转流程(即，给定内部表示，获取输入)，则有必要从最顶层进行采样，应用以下公式:

当然，这个过程必须向后重复，直到到达实际的输入层。数据库网络非常强大(例如，在天体物理学领域有几个科学应用)，即使它们的结构不像其他更近的模型那样灵活；然而，复杂性通常更高，因此，我总是建议从较小的模型开始，只有在最终精度不足以满足特定目的的情况下，才增加层和/或神经元的数量。

无人监管的 DBN 的例子

在本例中，我们希望使用 DBN 来寻找 MNIST 数据集的低维表示。由于这些模型的复杂性很容易增加，我们将把过程限制在 500 个随机样本。该实现基于深信念网络包(github.com/albertbup/d…，支持 NumPy 和 TensorFlow。在前一种情况下，类(名称保持不变)必须从dbn包导入，而在后一种情况下，包是dbn.tensorflow。在本例中，我们将使用需求较少的 NumPy 版本，但也邀请读者查看 TensorFlow 版本。

让我们从加载和规范化数据集开始，如下所示:

import numpy as np

from sklearn.datasets import load_digits
from sklearn.utils import shuffle

nb_samples = 500

digits = load_digits()

X_train = digits['data'] / np.max(digits['data'])
Y_train = digits['target']

X_train, Y_train = shuffle(X_train, Y_train, random_state=1000)
X_train = X_train[0:nb_samples]
Y_train = Y_train[0:nb_samples]

我们现在可以实例化UnsupervisedDBN类，其结构如下:

64 个输入神经元(从数据集中隐含检测到)
32 个乙状结肠神经元
32 个乙状结肠神经元
16 个乙状结肠神经元

因此，最后一个表示由 16 个值组成(原始维度的四分之一)。我们正在设置 η=0.025 的学习率，每批 16 个样本(当然，为了尽量减少重构误差，请你查看其他配置)。以下代码片段初始化并训练了模型:

from dbn import UnsupervisedDBN

unsupervised_dbn = UnsupervisedDBN(hidden_layers_structure=[32, 32, 16],
                                   learning_rate_rbm=0.025,
                                   n_epochs_rbm=500,
                                   batch_size=16,
                                   activation_function='sigmoid')

X_dbn = unsupervised_dbn.fit_transform(X_train)

在训练过程的最后，我们可以在将分布投影到二维空间后，对其进行分析。像往常一样，我们将使用 t-SNE 算法，它保证找到最相似的低维分布:

from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=10, random_state=1000)
X_tsne = tsne.fit_transform(X_dbn)

下图显示了投影样本的图表:

t-SNE plot of the unsupervised DBN output representations

正如你所看到的，大部分的块是相当内聚的，这表明一个数字的特殊属性已经在低维空间中被成功地表现出来了。在某些情况下，同一个数字组被分成更多的簇，但一般来说，有噪声的(孤立的)点的数量极低。例如，包含数字2的组用符号x表示。大部分样品在0<x₀T32】30，x₁T33-40范围内；但是，一个子组也位于范围*-10<x₁<10*内。如果我们检查这个小簇的邻居，他们是由代表数字8(由正方形代表)的样本组成的。很容易理解，一些格式错误的二进制与格式错误的八进制非常相似，这证明了原始簇的拆分是正确的。从统计学的角度来看，解释的差异可能会有不同的影响。在某些情况下，几个组件足以确定一个类的特殊特征，但这通常是不正确的。当属于不同类别的样本显示出相似性时，只能通过次要成分的差异来进行区分。当处理包含几乎(甚至部分)重叠样本的数据集时，这一考虑非常重要。当执行降维时，数据科学家的主要任务不是检查总体解释的方差，而是理解是否存在受降维负面影响的区域。在这种情况下，可以定义多个检测规则(例如，当一个样本,*x_I∈R₁或x_I∈R₄→x_I*具有 y _k 标签)或尽量避免创建这种分割的模型(例如

摘要

在这一章中，我们讨论了一些非常常见的用于解决无监督任务的神经模型。自动编码器允许您找到数据集的低维表示，而不会对其复杂性有特定的限制。特别地，深度卷积网络的使用有助于检测和学习高级和低级几何特征，当内部代码也比原始维度短得多时，这可以导致非常精确的重构。我们还讨论了如何为自动编码器增加稀疏性，以及如何使用这些模型对样本进行降噪。标准自动编码器的一个稍有不同的变体是变分自动编码器，它是一个生成模型，可以提高学习数据生成过程的能力，数据集应该从该过程中绘制。

Sanger 和 Rubner-Tavan 的网络是神经模型，能够在没有任何统计预处理的情况下提取数据集的第一个 k 主成分。它们还具有以在线方式自然工作的优势(尽管标准主成分分析通常需要整个数据集，即使存在性能略差于离线算法的增量变体)，以及以降序提取组件的优势。我们讨论的最后一个模型是无监督环境下的数据库网络。我们描述了它们的构造块 RBM 的生成属性，然后我们分析了这样的模型如何学习数据生成过程的内部(通常是低维的)表示。

在下一章中，我们将讨论其他神经模型— 生成性对抗网络 ( GANs )和自组织映射 ( SOMs )。前者可以学习输入分布，并从中产生新的样本，而后者基于大脑某些特定区域的功能，并训练它们的单位，以便接受特定的输入模式。

问题

在自动编码器中，编码器和解码器必须结构对称。这是正确的吗？
给定数据集 X 及其变换 Y ，基于自动编码器产生的代码， X 中包含的所有信息都可以在 Y 中找到。这是正确的吗？
一码，z_I∑(0，1)¹²⁸T5】有和(z _i ) = 36* 。稀疏吗？*
如果 std(z _i ) = 0.03 ，是不是代码稀疏？
桑格网络需要协方差矩阵列作为输入向量。这是正确的吗？
我们如何确定由 Rubner-Tavan 网络提取的每个组件的重要性？
给定一个随机向量，h_I∈ℜ^mt5】(m是一个 DBN 的输出维数)，有可能确定最可能对应的输入样本吗？

进一步阅读

堆叠去噪自动编码器:利用局部去噪标准在深度网络中学习有用的表示，文森特，P. ，拉罗彻尔，h，拉约伊，I，本吉奥，y，和曼扎戈尔，P. ，机器学习研究杂志 11 ，2010
稀疏自动编码器，CS294A， Ng，A. ，斯坦福大学
自动编码变分贝叶斯，金马。D. P】和威林* *，m .ArXiv:1312.6114【stat。ML]
理论神经科学，达扬和阿博特，洛杉矶，麻省理工学院出版社，2005 年
单层线性前馈神经网络中的最优无监督学习，神经网络，桑格，T. D. ， 2，1989
用于主成分分析的自组织网络，欧洲物理学快报，鲁布纳，j .和塔万，P. ， 10(7)，1989
动力系统中的信息处理:和声理论的基础，并行分布式处理，斯摩棱斯基，保罗，第一卷，麻省理工学院出版社，1986 年
训练受限玻尔兹曼机器实用指南，辛顿，G. ，多伦多大学计算机科学系，2010
深度信念网络的快速学习算法，辛顿，奥辛德罗，和特怀威，神经计算，18/7，2005
机器学习算法，第二版，博纳科尔索，g，Packt，2018
掌握机器学习算法，博纳科索，g，Packt，2018