R-应用无监督学习-二-R 应用无监督学习（二）第六章：第五章数据比较方法学习目标到本章结束时，你将能够：创建

R 应用无监督学习（二）

原文：annas-archive.org/md5/4c0eaccb71c3cc25e04cea28e06cad11

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章：第五章

数据比较方法

学习目标

到本章结束时，你将能够：

创建数据哈希
创建图像签名
比较图像数据集
执行因子分析以隔离潜在变量
使用因子分析比较调查和其他数据集

在本章中，我们将探讨不同的数据比较方法。

简介

无监督学习关注于分析数据的结构以得出有用的结论。在本章中，我们将探讨使我们能够利用数据结构来比较数据集的方法。我们将重点研究的方法包括哈希函数、分析签名和潜在变量模型。

哈希函数

假设你想把一个 R 脚本发送给你的朋友。然而，你和你的朋友在文件上遇到了技术问题——也许你们的电脑被恶意软件感染了，或者也许有黑客正在篡改你的文件。所以，你需要一种方法来确保你的脚本在发送给朋友时是完整的，没有被损坏或更改。检查文件是否完整的一种方法就是使用哈希函数。

哈希函数可以为数据创建类似指纹的东西。我们所说的指纹是指一种小而易于检查的东西，使我们能够验证数据是否具有我们认为是其身份的东西。因此，在你创建想要发送的脚本之后，你将对脚本应用哈希函数并获取其指纹。然后，你的朋友可以在收到文件后使用相同的哈希函数对文件进行处理，并确保指纹匹配。如果发送的文件指纹与接收到的文件指纹匹配，那么这两个文件应该是相同的，这意味着文件是完整发送的。以下练习展示了如何创建和使用一个简单的哈希函数。

练习 29：创建和使用哈希函数

在这个练习中，我们将创建和使用一个哈希函数：

指定需要使用哈希函数的数据。我们一直在探讨你想发送的 R 脚本场景。以下是一个简单的 R 脚本示例：
```
string_to_hash<-"print('Take the cake')"
```
在这里，我们有一个打印字符串Take the cake的脚本。我们将其保存为名为string_to_hash的变量。
指定可能的哈希值总数。我们的目标是为我们脚本创建一个指纹。我们需要指定我们将允许存在的指纹总数。我们希望指定一个足够低以便于操作，但又足够高以至于不同脚本偶然具有相同指纹的可能性不大的数字。在这里，我们将使用 10,000：
```
total_possible_hashes<-10000
```
将脚本（目前是一个字符串）转换为数值。哈希函数通常是算术性的，因此我们需要用数字而不是字符串来工作。幸运的是，R 有一个内置函数可以为我们完成这项工作：
```
numeric<-utf8ToInt(string_to_hash)
```
这已经将我们的脚本中的每个字符转换为整数，基于每个字符在 UTF-8 编码方案中的编码。我们可以通过打印到控制台来查看结果：
```
print(numeric)
```
输出如下：
```
  [1] 112 114 105 110 116  40  39  84  97 107 101  32 116 104 101  32  99  97 107[20] 101  39  41
```
应用我们的散列函数。我们将使用以下函数来生成我们的最终散列，或者说，脚本的指纹：
```
hash<-sum((numeric+123)²) %% total_possible_hashes
```
在 R 中运行这一行后，我们发现hash的最终值是 2702。由于我们使用了模运算符（在 R 中为%%），hash的值将始终在 0 到 10,000 之间。

我们用来将数值向量转换为最终散列值的简单函数并不是唯一的散列函数。专家们设计了许多具有不同复杂程度的此类函数。好的散列函数应该具有许多属性，其中之一是最重要的属性是抗碰撞性，这意味着很难找到两个产生相同散列值的数据集。

练习 30：验证我们的散列函数

在这个练习中，我们将验证我们的散列函数是否使我们能够有效地比较不同的数据，通过检查不同的信息产生不同的散列值。这个练习的结果将是不同信息的散列函数，我们将比较以验证不同的信息产生不同的散列值：

创建一个执行散列的函数。我们将把前面练习中引入的代码组合到一个函数中：
```
get_hash<-function(string_to_hash, total_possible_hashes){
numeric<-utf8ToInt(string_to_hash)
hash<-sum((numeric+123)²) %% total_possible_hashes
return(hash)
}
```
这个函数接收我们想要散列的字符串和可能的散列总数，然后应用我们在上一个练习中使用的相同的散列计算来计算散列值，并返回该散列值。
比较不同输入的散列值。我们可以如下比较来自不同输入的散列值：
```
script_1<-"print('Take the cake')"
script_2<-"print('Make the cake')"
script_3<-"print('Take the rake')"
script_4<-"print('Take the towel')"
```
在这里，我们有四个不同的字符串，表达不同的信息。我们可以看到它们的散列值如下：
```
print(get_hash(script_1,10000))
print(get_hash(script_2,10000))
print(get_hash(script_3,10000))
print(get_hash(script_4,10000))
```
第一个脚本返回散列值 2702，正如我们在前面的练习中找到的那样。第二个脚本，尽管它只与第一个脚本有一个字符不同，返回的散列值是 9853，第三个脚本，也只与第一个脚本有一个字符不同，返回的散列值是 9587。最后一个脚本返回的散列值是 5920。尽管这四个脚本有相当大的相似性，但它们有不同的指纹，我们可以使用这些指纹来比较和区分它们。

这些散列对于您和您的消息接收者来说很有用，可以用来验证您的脚本在发送过程中未被篡改。当您发送脚本时，您可以告诉您的朋友确保脚本的散列值是 2702。如果您的朋友收到的脚本散列值不是 2702，那么您的朋友可以得出结论，脚本在发送和接收过程中被篡改了。如果您的朋友能够可靠地检测文件是否损坏，那么您可以避免传播恶意软件或与您的朋友因误解而争吵。

在线分发的软件有时会附带一个用户可以用来检查文件损坏的哈希值。为此，专业人士使用比前面练习中简单函数更高级的哈希函数。专业人士使用的其中一个哈希函数称为 MD5，在 R 中使用digest包可以非常容易地应用它：

install.packages('digest')
library(digest)
print(digest(string_to_hash,algo='md5'))

输出如下：

[1] "a3d9d1d7037a02d01526bfe25d1b7126"

在这里，我们可以看到我们简单 R 脚本的 MD5 哈希值。您可以自由尝试其他数据的 MD5 哈希值，以比较结果。

分析签名

在第四章，降维中，我们讨论了降维——这些方法使我们能够以给我们数据洞察力的方式简洁地表达数据。之前讨论的哈希函数是另一种实现降维的方法。哈希函数在许多用途中都很有用，包括我们讨论的文件验证用例。在那个场景中，我们感兴趣的是确定两个脚本是否完全相同。即使数据有细微的差异，例如将“take”一词改为“make”，也可能完全破坏预期的信息，因此需要精确性。

在其他情况下，我们可能希望在不需要比较的两个数据集具有完全相同性的情况下，进行有意义的比较。考虑检测版权侵犯的情况。假设一个网站托管了来自其用户的图像。它想要确保用户没有提交受版权保护的图像。因此，每次它收到上传时，它都希望检查该上传是否与大型版权图像数据库中的任何图像相同。

仅检查图像是否完全相同是不够的，因为一些不择手段的上传者可能会进行细微的修改并仍然尝试上传。例如，他们可能会改变一个像素的颜色，或者非常轻微地裁剪图像，或者将其压缩到比原始图像大或小。即使有这些细微的修改，图像仍然会违反版权法。一个检查完全相同性的哈希函数将无法识别这些版权侵犯。因此，图像托管网站将想要检查代表两张图片的数据中是否存在任何实质上相似的底层结构。

我们提到，哈希函数创建了一种类似于数据指纹的东西。指纹应该在每次观察时都完全相同，即使两个指纹在大多数方面相似，除非它们完全匹配，否则我们不会认为它们彼此匹配。

在这种情况下，我们需要的更像是一个签名，而不是指纹。每次您签名时，您的签名应该看起来或多或少相同。但是，即使是同一个人在尝试匹配之前的签名时，每个签名之间也会有细微的差异。为了验证签名是否匹配，我们需要检查实质性的相似性，而不是完美的身份。我们在这里提供的代码将展示如何将任何大小图像编码为一个小巧且健壮的编码，这允许在数据集之间进行快速且准确的近似比较。这种编码图像数据的方法可以被称为创建分析签名。

我们创建图像签名的步骤如下。首先，我们将图像分成一个 10x10 的网格。然后，我们将测量网格每个部分的亮度。之后，我们将比较每个部分的亮度与其相邻部分的亮度。最终的签名将包含一个向量，该向量包含每个网格部分与其每个相邻部分的比较。这种方法是由 Wong、Bern 和 Goldberg 发明的，并在一篇名为《任何类型图像的图像签名》的论文中发表。

在我们创建分析签名之前，我们需要进行一些数据准备，如下一练习所述。

练习 31：为创建图像分析签名进行数据准备

在这个练习中，我们将为阿拉莫的照片创建分析签名进行数据准备。

注意

对于所有需要导入外部 CSV 文件或图像的练习和活动，请转到RStudio-> 会话-> 设置工作目录-> 到源文件位置。您可以在控制台中看到路径已自动设置。

首先，我们需要配置 R 以能够读取并处理我们的图像数据。我们需要安装imager包。您可以在 R 控制台中执行install.packages('imager')来安装此包，然后您可以通过在控制台中运行library('imager')来加载它。
接下来，我们需要读取数据。在我们的例子中，我们将使用这张阿拉莫的照片：

图 5.1：阿拉莫图像

首先，从github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/Lesson05/Exercise31/alamo.jpg 下载到您的计算机上，并将其保存为alamo.jpg。确保它保存在 R 的工作目录中。如果它不在 R 的工作目录中，那么请使用setwd()函数更改 R 的工作目录。然后，您可以将此图像加载到名为im（代表图像）的变量中，如下所示：
```
filepath<-'alamo.jpg'
im <- imager::load.image(file =filepath) 
```
我们将要探索的其余代码将使用这个名为im的图像。在这里，我们已经将阿拉莫的照片加载到im中。然而，你可以通过将图像保存到你的工作目录并在filepath变量中指定其路径来运行其余的代码。
我们正在开发的签名是为了用于灰度图像。因此，我们将使用imager包中的函数将此图像转换为灰度：
```
im<-imager::rm.alpha(im)
im<-imager::grayscale(im)
im<-imager::imsplit(im,axis = "x", nb = 10)   
```
这段代码的第二行是将图像转换为灰度。最后一行将图像分割成 10 等份。
以下代码创建了一个空矩阵，我们将用有关我们 10x10 网格每个部分的详细信息来填充它：
```
matrix <- matrix(nrow = 10, ncol = 10)
```
接下来，我们将运行以下循环。这个循环的第一行使用了imsplit命令。这个命令之前也被用来将 x 轴分成 10 等份。这次，对于 x 轴的每个 10 等份，我们将在 y 轴上进行分割，也将它分成 10 等份：
```
for (i in 1:10) {
  is <- imager::imsplit(im = im[[i]], axis = "y", nb = 10)
  for (j in 1:10) {
    matrix[j,i] <- mean(is[[j]])
  }
}
```
在沿 y 轴分割后，矩阵通过mean(is[[j]])更新。这是所选部分的平均亮度的度量。

这段代码的结果是一个 10x10 的矩阵，其中i-j元素包含原始照片i-j部分的平均亮度。

如果你打印这个矩阵，你可以看到照片每个部分的亮度数字：
```
print(matrix)
```
输出应该看起来像以下这样：

![Figure 5.2: Screenshot of the output matrix img/C12628_05_02.jpg

图 5.2：输出矩阵的截图

你可以将这些亮度数字与原始照片的外观进行比较。

我们可以在这里停止，因为我们已经生成了一个复杂数据集的压缩编码。然而，我们可以采取一些进一步的步骤来使这个编码更有用。

我们可以做的事情之一是创建一个brightnesscomparison函数。这个函数的目的是比较图像中两个不同部分的相对亮度。最终，我们将比较我们分析的每张图像的所有不同部分。我们的最终指纹将包含许多这样的亮度比较。这个练习的目的是创建一个亮度比较函数，这将最终使我们能够创建最终的指纹。

请注意，这个练习是在上一个练习的基础上构建的，这意味着你应该在运行这个练习的代码之前运行上一个练习中的所有代码。

练习 32：创建亮度比较函数

在这个函数中，我们传递两个参数，x和y：每个参数代表图片中特定部分的亮度。如果x和y相当相似（小于 10%的差异），那么我们可以说它们基本上是相同的，我们返回 0，表示亮度差异大约为 0。如果x比y大 10%以上，我们返回 1，表示x比y亮，如果x比y小 10%以上，我们返回-1，表示x比y暗。

创建亮度比较函数。brightnesscomparison 函数的代码如下：

brightnesscomparison<-function(x,y){
compared<-0
if(abs(x/y-1)>0.1){
if(x>y){
compared<-1
}
if(x<y){
compared<-(-1)
}
}
return(compared)
}

我们可以使用这个函数来比较我们为图片形成的 10x10 网格的两个部分。例如，为了找到与直接左侧部分的亮度比较，我们可以执行以下代码：
```
i<-5
j<-5
left<-brightnesscomparison(matrix[i,j-1],matrix[i,j])
```
在这里，我们查看矩阵的第 5 行和第 5 列。我们将这个部分与它左侧直接的部分进行比较——第 5 行和第 4 列，我们通过指定 j-1 来访问这个部分。
使用亮度比较函数来比较图像部分与其上面的邻居。我们可以进行类似的操作来比较这个部分与其上面的部分：
```
i<-5
j<-5
top<-brightnesscomparison(matrix[i-1,j],matrix[i,j])
```

在这里，top 是第五部分与它上面紧邻的节点的亮度比较，我们通过指定 i-1 来访问这个节点。

这个练习的重要输出是 top 和 left 的值，它们都是图像部分与其他相邻部分的比较。在这种情况下，left 等于零，意味着我们选择的图像部分的亮度与左侧的图像部分大致相同。同样，top 等于 1，意味着我们选择的节点的直接上方部分比我们选择的节点亮度更高。

在下一个练习中，我们将创建一个 neighborcomparison 函数。这个函数将比较我们 10x10 网格中的每个节点的亮度，与它的邻居进行比较。这些邻居包括我们刚才比较过的左侧邻居，以及上面的邻居。总的来说，我们图片的每个部分（顶部、底部、左侧、右侧、左上、右上、左下和右下）都有八个邻居。我们想要这个邻居比较函数的原因是，它将使我们很容易得到最终的解析特征。

请注意，这个练习建立在之前的练习之上，你应该在运行这段代码之前运行所有之前的代码。

练习 33：创建一个函数来比较图像部分与所有相邻部分

在这个练习中，我们将创建一个 neighborcomparison 函数来比较图像部分与其他所有相邻部分。为此，执行以下步骤：

创建一个函数，比较图像部分与其左侧的邻居。我们在之前的练习中做过这个操作。对于任何图像部分，我们可以比较其亮度与其左侧邻居的亮度，如下所示：
```
i<-5
j<-5
left<-brightnesscomparison(matrix[i,j-1],matrix[i,j])
```
创建一个函数，比较图像部分与其上面的邻居。我们在之前的练习中做过这个操作。对于任何图像部分，我们可以比较其亮度与其上方邻居的亮度，如下所示：
```
i<-5
j<-5
top<-brightnesscomparison(matrix[i-1,j],matrix[i,j])
```
如果你查看 Step 1和 Step 2，你可以开始注意到这些邻居比较中的模式。要比较图像部分与其左侧的部分，我们需要访问矩阵的j-1索引部分。要比较图像部分与其右侧的部分，我们需要访问矩阵的j+1索引部分。要比较图像部分与其上方的部分，我们需要访问矩阵的i-1索引部分。要比较图像部分与其下方的部分，我们需要访问矩阵的i+1索引部分。

因此，我们将有每个图像部分与其上方、下方、左侧和右侧的每个邻居的比较。以下代码显示了我们将进行的比较，除了在 Step 1和 Step 2中进行的顶部和左侧比较：
```
i<-5
j<-5
top_left<-brightnesscomparison(matrix[i-1,j-1], matrix[i,j])
bottom_left<-brightnesscomparison(matrix[i+1,j-1],matrix[i,j])
top_right<-brightnesscomparison(matrix[i-1,j+1],matrix[i,j])
right<-brightnesscomparison(matrix[i,j+1],matrix[i,j])
bottom_right<-brightnesscomparison(matrix[i+1,j+1],matrix[i,j])
bottom<-brightnesscomparison(matrix[i+1,j],matrix[i,j])
```
初始化一个向量，该向量将包含部分与其每个邻居的最终比较：
```
comparison<-NULL
```
我们将使用这个comparison向量来存储我们最终生成的所有邻居比较。它将包含图像部分与其每个邻居的比较。

在步骤 1-4 中，我们展示了邻居比较函数的各个部分。在这个步骤中，我们将它们组合起来。你可以看到的邻居比较函数接受一个图像矩阵作为参数，并且还有i和j值，指定我们正在关注的图像矩阵的部分。该函数使用我们为top和left比较编写的代码，并为其他邻居添加了其他比较，例如top_left，它比较图像亮度级别与上方左侧部分的图像亮度。总的来说，每个图像部分应该有八个邻居：左上、上、右上、左、右、左下、下和右下。在这个步骤中，我们将进行这八个比较并将它们存储在comparison向量中。最后，有一个return语句，它返回所有比较。

这里是我们可以使用来获取所有邻居比较的函数：

neighborcomparison<-function(mat,i,j){
comparison<-NULL
top_left<-0
if(i>1 & j>1){
top_left<-brightnesscomparison(mat[i-1,j-1],mat[i,j])
}
left<-0
if(j>1){
left<-brightnesscomparison(mat[i,j-1],mat[i,j])
}
bottom_left<-0
if(j>1 & i<nrow(mat)){
bottom_left<-brightnesscomparison(mat[i+1,j-1],mat[i,j])
}
top_right<-0
if(i>1 & j<nrow(mat)){
top_right<-brightnesscomparison(mat[i-1,j+1],mat[i,j])
}
right<-0
if(j<ncol(mat)){
right<-brightnesscomparison(mat[i,j+1],mat[i,j])
}
bottom_right<-0
if(i<nrow(mat) & j<ncol(mat)){
bottom_right<-brightnesscomparison(mat[i+1,j+1],mat[i,j])
}
top<-0
if(i>1){
top<-brightnesscomparison(mat[i-1,j],mat[i,j])
}
bottom<-0
if(i<nrow(mat)){
bottom<-brightnesscomparison(mat[i+1,j],mat[i,j])
}
comparison<-c(top_left,left,bottom_left,bottom,bottom_right,right,top_right,top)
return(comparison)
}

此函数返回一个包含八个元素的向量：每个元素对应于我们网格中特定部分的邻居。你可能已经注意到，10x10 网格的一些部分似乎没有八个邻居。例如，10x10 网格的 1-1 元素下面有一个邻居，但没有上面的邻居。对于没有特定邻居的网格位置，我们说它们对该邻居的亮度比较为 0。这降低了创建和解释亮度比较的复杂性水平。

最终输出是一个名为comparison的向量，它包含图像部分与其八个邻居之间的亮度级别比较。

在下一个练习中，我们将完成创建我们的分析签名。每个图像的分析签名将包括每个图像部分与其八个邻居的比较。我们将使用两个嵌套的 for 循环来遍历我们 10x10 网格的每个部分。本练习的预期输出将是一个生成图像分析签名的函数。

练习 34：创建一个生成图像分析签名的函数

在这个练习中，我们将创建一个函数，用于为图像生成一个分析签名。为此，请执行以下步骤：

我们首先创建一个名为 signature 的变量，并用 NULL 值初始化它：
```
signature<-NULL
```
当我们完成时，这个 signature 变量将存储完整的签名。
现在，我们可以遍历我们的网格。对于网格的每个部分，我们向签名中添加八个新元素。我们添加的元素是我们之前介绍的 neighborcomparison 函数的输出：
```
for (i in 1:nrow(matrix)){
for (j in 1:ncol(matrix)){
signature<-c(signature,neighborcomparison(matrix,i,j))
}
}
```
我们可以通过在控制台中运行 print(signature) 来查看我们的指纹是什么样的。它是一个包含 800 个值的向量，所有这些值都等于 0（表示相似的亮度或没有邻居）、1（表示某个区域比其邻居更亮）或 -1（表示某个区域比其邻居更暗）。
将 Step 1 和 Step 2 结合在一个函数中，该函数可以生成任何图像矩阵的签名：
```
get_signature<-function(matrix){
signature<-NULL
for (i in 1:nrow(matrix)){
for (j in 1:ncol(matrix)){
signature<-c(signature,neighborcomparison(matrix,i,j))
}
}
return(signature)
}
```
此代码定义了一个名为 get_signature 的函数，它使用来自 Step 1 和 Step 2 的代码来获取该签名。我们可以使用之前创建的图像矩阵来调用此函数。
由于我们稍后还将创建更多签名，我们将把这个签名保存到一个变量中，该变量指明了它代表的是什么。在这种情况下，我们将称之为 building_signature，因为它是一个建筑图像的签名。我们可以这样做：
```
building_signature<-get_signature(matrix)
building_signature
```
输出如下：

图：5.3：building_signature 矩阵

存储在 building_signature 中的向量是这个练习的最终输出，也是我们在本章中一直试图开发的图像签名。

这个签名旨在类似于人类的亲笔签名：小巧且表面上与其他签名相似，但足够独特，使我们能够将其与数百万其他现有签名区分开来。

我们可以通过读取一个完全不同的图像并比较生成的签名来检查我们找到的签名解决方案的鲁棒性。这就是以下活动的场景。

活动 11：为人物照片创建图像签名

让我们尝试为这张图像创建一个图像指纹，这是一张伟大的豪尔赫·路易斯·博尔赫斯的照片。

要完成这个活动，您可以遵循本章中我们已经遵循的所有步骤。以下步骤将为您概述这个过程。请记住，在我们之前进行的图像签名练习中，我们使用了一个 10x10 的亮度测量矩阵。然而，10x10 的矩阵可能不适合某些情况，例如，如果我们处理的图像特别小，或者我们有数据存储限制，或者我们期望使用不同的矩阵大小可以获得更高的精度。因此，在以下活动中，我们也将使用 9x9 矩阵计算一个签名：

注意

这可以在任何给定的矩阵大小上执行。它可能是一个 5x5 的矩阵用于数据存储，或者是一个 20x20 的矩阵用于精确的签名。

图 5.4：豪尔赫·路易斯·博尔赫斯图像

这些步骤将帮助我们完成活动：

将图像加载到您的 R 工作目录中。将其保存到名为 im 的变量中。
将您的图像转换为灰度并分成 100 个部分。
创建一个亮度值矩阵。
使用我们之前创建的 get_signature 函数创建一个签名。将 Borges 图像的签名保存为 borges_signature。

注意

这个活动的解决方案可以在第 227 页找到。

这个活动的最终输出是 borges_signature 变量，它是 Borges 照片的分析签名。此外，我们还创建了 borges_signature_ninebynine 变量，它也是一个分析签名，但基于 9x9 而不是 10x10 的矩阵。我们可以在分析中使用它们中的任何一个，但我们将使用 borges_signature 变量。如果您已经完成了迄今为止的所有练习和活动，那么您应该有两个分析签名：一个名为 building_signature，另一个名为 borges_signature。

签名比较

接下来，我们可以比较这两个签名，看看它们是否将我们的不同图像映射到不同的签名值。

您可以使用以下一行 R 代码比较签名：

comparison<-mean(abs(borges_signature-building_signature))

这种比较计算了两个签名中每个元素之间的差的绝对值，然后计算这些值的平均值。如果两个签名完全相同，那么这个差值将为 0。comparison 的值越大，两个图像的差异就越大。

在这种情况下，comparison 的值为 0.644，这表明平均而言，相应的签名条目之间大约相差 0.644。对于值仅在 1 和 -1 之间变化的数据库来说，这种差异是显著的。因此，我们看到我们的签名创建方法为非常不同的图像创建了非常不同的签名，正如我们所期望的那样。

现在，我们可以计算一个与我们的原始图像非常相似但又不完全相同的图像的签名：

图 5.5：标记的阿拉莫图像

为了获得这张图像，我从一个原始的阿拉莫图像开始，并在四个地方添加了单词水印，模拟了有人可能对图像进行的修改。由于图像现在与原始图像不同，简单的版权检测软件可能会被这个水印欺骗。我们的分析签名方法不应该如此天真。我们将在以下活动中完成这项任务。

活动十二：为水印图像创建图像签名

在这个活动中，我们将为水印图像创建一个图像签名：

将图像加载到您的 R 工作目录中。将其保存到名为im的变量中。
将您的图像转换为灰度并分成 100 个部分。
创建一个亮度值矩阵。
使用我们之前创建的get_signature函数创建一个签名。将阿拉莫图像的签名保存为watermarked_signature。
将水印图像的签名与原始图像的签名进行比较，以确定签名方法是否能够区分图像。

输出将如下所示：

![图 5.6：预期水印图像的签名 img/C12628_05_06.jpg

图 5.6：预期水印图像的签名

注意

本活动的解决方案可以在第 230 页找到。

为了检测版权侵权，我们可以计算数据库中每个受版权保护图像的签名。然后，对于每个新上传的图像，我们将新上传图像的签名与数据库中的签名进行比较。如果任何受版权保护的图像的签名与新上传的图像的签名相同或非常接近，我们将它们标记为潜在的匹配项，需要进一步调查。比较签名可能比比较原始图像快得多，并且它具有对水印等小变化具有鲁棒性的优点。

我们刚才执行签名方法是一种将数据编码以实现不同数据集之间比较的方法。有许多其他编码方法，包括一些使用神经网络的方法。每种编码方法都将具有其独特的特征，但它们都将共享一些共同特征：它们都将尝试返回压缩数据，以便于在不同数据集之间进行简单和准确的比较。

将其他无监督学习方法应用于分析签名

到目前为止，我们只使用了哈希和分析签名来比较两张图像或四个简短字符串。然而，可以应用于哈希或分析签名的无监督学习方法没有限制。为一系列图像创建分析签名可能是深入分析的第一步，而不是最后一步。在为一系列图像创建分析签名之后，我们可以尝试以下无监督学习方法：

聚类：我们可以将第一、二章节中讨论的任何聚类方法应用于由分析特征组成的数据库。这可能使我们能够找到所有倾向于彼此相似的一组图像，可能是因为它们是同一类型物体的照片。有关聚类方法的更多信息，请参阅第一、二章节。
异常检测：我们可以将第六章中描述的异常检测方法应用于由分析特征组成的数据库。这将使我们能够找到与数据集中其他图像非常不同的图像。

潜在变量模型 – 因子分析

本节将介绍潜在变量模型。潜在变量模型试图用少量隐藏或潜在的变量来表示数据。通过找到与数据集相对应的潜在变量，我们可以更好地理解数据，甚至可能理解数据来源或生成方式。

考虑到学生在各种不同课程中获得的分数，从数学到音乐、外语到化学。心理学家或教育工作者可能对使用这些数据更好地理解人类智力感兴趣。研究人员可能想要在数据中测试几种不同的智力理论，例如：

理论 1：有两种不同类型的智力，拥有一种类型的人将在一组课程中表现出色，而拥有另一种类型的人将在其他课程中表现出色。
理论 2：只有一种类型的智力，拥有它的人将在所有类型的课程中表现出色，而没有它的人则不会。
理论 3：一个人可能在一个或几个他们所学的课程标准下非常聪明，但在其他标准下并不聪明，每个人都会有一套不同的课程，他们在这些课程中表现出色。

这些理论中的每一个都表达了一个潜在变量的概念。数据只包含学生成绩，这可能是智力的表现，但不是智力本身的直接衡量标准。这些理论表达了不同类型的智力如何以潜在的方式影响成绩。即使我们不知道哪个理论是正确的，我们也可以使用无监督学习工具来理解数据的结构以及哪个理论最适合数据。

任何记得自己曾是学生的人可能都厌倦了别人评价他们的智力。因此，我们将使用一个稍微不同的例子。我们不会评价智力，而是评价个性。我们不会寻找左右脑，而是寻找人们个性的不同特征。但我们将采取相同的方法——使用潜在变量来识别复杂数据是否可以由少数几个因素解释。

为了为我们的分析做准备，我们需要安装和加载正确的包。下一个练习将涵盖如何加载因子分析所需的包。

练习 35：准备因子分析

在这个练习中，我们将为因子分析准备数据：

安装必要的包：

install.packages('psych')
install.packages('GPArotation')
install.packages('qgraph')

然后，你可以运行以下行将它们加载到你的 R 工作空间中：
```
library(psych)
library(GPArotation)
library(qgraph)
```
接下来，加载数据。我们将使用的数据是 500 份对名为“修订版 NEO 人格问卷”的人格测试的记录。该数据集包含在名为 qgraph 的 R 包中。首先，将数据读入你的 R 工作空间。以下代码将允许你将其作为名为 big5 的变量访问：
```
data(big5)
```
你可以通过在 R 控制台中运行以下行来查看数据的顶部部分：
```
print(head(big5))
```
输出如下：

![图 5.7：数据顶部部分]

图 5.7：数据顶部部分
你可以通过以下方式查看你的数据行和列的数量：
```
print(nrow(big5))
```
输出如下：
```
500
```
要检查列，执行以下操作：
```
print(ncol(big5))
```
输出如下：
```
240
```
这份数据包含 500 行和 240 列。每一行是与一个调查受访者相关的完整记录。每一列记录了一个问题的答案。因此，有 500 名调查受访者每人回答了 240 个问题。每个答案都是数值型的，你可以看到以下回答的范围：
```
print(range(big5))
```
输出如下：
```
[1] 1 5
```
答案的范围是 1 到 5。

本练习的预期输出是一个已加载 big5 数据的 R 工作空间。你可以确信数据已加载，因为当你控制台中运行 print(range(big5)) 时，你会得到 1 5 作为输出。

在这个调查中，问题是通过向调查受访者展示一个问题并询问他们同意该答案的程度来提出的，其中 5 代表强烈同意，1 代表强烈不同意。这些陈述是关于一个人特定特征的说法，例如：

“我喜欢结识新朋友。”

“我有时会犯错误。”

“我喜欢修理东西。”

你可以想象，其中一些问题可能会测量相同的潜在“潜在”人格特质。例如，我们经常描述人们为好奇，从他们喜欢学习和尝试新事物的意义上来说。那些强烈同意“我喜欢结识新朋友”这一说法的人，如果这两个说法都在测量潜在的好奇心，那么他们很可能也会同意“我喜欢修理东西”这一说法。或者，可能是有些人具有对人们感兴趣的人格特质，而其他人则具有对事物感兴趣的人格特质。如果是这样，那些同意“我喜欢结识新朋友”的人就不太可能也同意“我喜欢修理东西”。在这里，因子分析的目的在于发现哪些问题对应于一个潜在的想法，并了解这些潜在想法是什么。

你可以看到，大多数列都以字母表中的一个字母以及一个数字开头。目前，你可以忽略这些标签，只需记住每个问题都有细微的差别，但其中许多是相似的。

我们现在准备好开始我们的潜在变量模型。在这种情况下，我们将执行因子分析。因子分析是一种强大且非常常见的方法，可以执行我们之前讨论过的潜在变量分析。因子分析假设存在某些潜在因子控制数据集，然后向我们展示这些因子是什么以及它们如何与我们的数据相关。

要开始我们的因子分析，我们需要指定我们要找多少个因子。有几种方法可以决定我们应该寻找多少个因子。目前，我们将从查看五个因子开始，因为这项数据被称为大五因子，然后我们稍后会检查是否其他数量的因子更好。

我们需要为我们新的数据创建一个相关矩阵。相关矩阵是一个矩阵，其中i-j项是存储在列i中的变量和存储在列j中的变量之间的相关系数。这与我们在第四章，降维中讨论的协方差矩阵类似。你可以按照以下方式创建这个数据的相关矩阵：

big_cor <- cor(big5)

现在，因子分析本身相当简单。我们可以通过使用fa命令来完成它，这是psych R 包的一部分：

solution <- fa(r = big_cor, nfactors = 5, rotate = "oblimin", fm = "pa")

你会注意到这个命令中有几个参数。首先，我们指定r=big_cor。在这种情况下，R 的psych包的作者决定使用小写r来指代因子分析中使用的协方差矩阵。下一个参数是nfactors=5。这指定了我们在因子分析中要寻找的因子数量。这次我们选择了五个因子，但稍后我们会看看使用更多或更少的因子。

最后两个参数对我们这里的目的是不太重要的。第一个参数说rotate="oblimin"。因子分析在向我们展示结果之前，在幕后对我们的数据进行旋转，这就是所谓的旋转。有许多技术可以用来完成这种幕后旋转，而oblimin是fa函数的作者选择作为默认的旋转方法。你可以自由地尝试其他旋转方法，但它们通常会产生实质上相似的结果。最后一个参数，就像旋转参数一样，指定了一个在幕后使用的方法。在这里，fm代表因子分解方法，pa代表主成分。你也可以尝试其他因子分解方法，但再次强调，它们应该会产生实质上相似的结果。

注意

你可以使用 R 中的?命令查找与任何其他命令相关的文档。在这种情况下，如果你对我们刚刚使用的fa命令感到好奇，你可以运行?fa来从你的 R 控制台加载文档。

现在，我们可以查看我们的因子分析输出：

print(solution)

输出如下：

图 5.8：输出的一部分

当我们这样做时，大部分输出是一个包含 240 行，标签为Standardized loadings的 DataFrame。这个数据框的每一行对应于我们原始数据框中的一列，例如N1。请记住，这些是我们数据来源的个性测试问题。这个数据框的前五行标签为PA1至PA5。这些对应于我们要找的五个因素。特定个性问题和特定因素的数量称为负载。例如，我们在数据框中有 0.54 的条目，对应于个性问题N1和因素PA1。我们说问题N1在因素PA1上的负载为 0.54。

我们可以将负载解释为对总分有贡献的因素。为了得到最终分数，你可以将每个特定负载乘以调查受访者对每个问题的回答，并将结果相加。用方程式表示，我们可以写成以下形式：

受访者 1 的因素 1 得分 =

（问题 1 上的因素 1 负载）*（受访者 1 对问题 1 的回答）+

（问题 2 上的因素 1 负载）*（受访者 1 对问题 2 的回答）+

...

（问题 240 上的因素 1 负载）*（受访者 1 对问题 240 的回答）

因此，如果一个特定问题的因素 1 负载很高，那么受访者对该问题的回答将对总因素 1 得分有较大的贡献。如果负载较低，那么该问题的回答对因素 1 得分的贡献就不那么大，换句话说，它对因素 1 的影响就不那么重要。这意味着每个负载都是衡量每个特定问题在因素测量中的重要程度的度量。

每个因素都有 240 个总负载——对应于个性调查中的每个问题。如果你查看负载矩阵，你可以看到许多问题在一个因素上有很大的负载，而在所有其他因素上的负载都很小（接近零）。研究人员经常试图根据每个因素的最高负载问题来假设每个因素的解释。

在我们的案例中，我们可以看到第一个因子，称为PA1，对第一个问题（标记为N1）有较高的负荷（0.54）。它还对第六个问题（N6）有相对较高的负荷（0.45），以及第十一题（N11）有较高的负荷（0.62）。在这里很容易看出一个模式——标记为N的问题往往对这个第一个因子有较高的负荷。结果证明，原始测试上的这些N问题都是为了衡量心理学家所说的“神经质”。神经质是一种基本的人格特质，使人们在面对困难时产生强烈的负面反应。调查中的N问题各不相同，但每个问题都是旨在衡量这种人格特质。我们在因子分析中发现，这个第一个因子往往对神经质问题有较高的负荷——因此我们可以称这个因子为“神经质因子”。

我们可以在其他负荷中看到类似的模式：

因素 2 似乎对“A”问题有较高的负荷，这些问题是用来衡量宜人性的。
因素 3 似乎对“E”问题有较高的负荷，这些问题是用来衡量外向性的。
因素 4 似乎对“C”问题有较高的负荷，这些问题是用来衡量尽责性的。
因素 5 似乎对“O”问题有较高的负荷，这些问题是用来衡量开放性的。

这些标签对应于“大五”人格理论，该理论认为这五种人格特质是人格最重要的最基本方面。在这种情况下，我们的五因素分析产生了与数据集中预先存在的标签相匹配的因子负荷模式。然而，我们不需要标记的问题来从因子分析中学习。如果我们获得了未标记的问题数据，我们仍然可以运行因子分析，并仍然可以在特定问题的特定因子负荷中找到模式。在找到这些模式后，我们必须仔细查看在相同因子上有高负荷的问题，并试图找出这些问题的共同之处。

因子分析背后的线性代数

因子分析是一种强大且灵活的方法，可以用多种方式使用。在以下练习中，我们将更改我们使用的因子分析命令的一些细节，以便更好地了解因子分析是如何工作的。

请注意：以下练习基于之前的因子分析代码。您必须运行之前突出显示的代码，特别是big_cor <- cor(big5)，才能成功运行以下练习中的代码。

练习 36：使用因子分析的进一步探索

在这个练习中，我们将详细探讨因子分析：

我们可以改变fa函数中的几个参数以获得不同的结果。记住上次我们运行以下代码来创建我们的解决方案：
```
solution <- fa(r = big_cor, nfactors = 5, rotate = "oblimin", fm = "pa") 
```
这次，我们可以更改fa函数的一些参数，如下所示：
```
solution <- fa(r = big_cor, nfactors = 3, rotate = "varimax", fm = "minres")
```
在这种情况下，我们将旋转方法更改为varimax，将因素化方法更改为minres。这些都是对函数背后使用的方法的改变。对我们来说最重要的是，我们将因素数量（nfactors）从 5 更改为 3。检查此模型中的因素载荷：
```
print(solution)
```
输出如下：

图 5.9：输出的一部分

你也可以尝试在这些载荷中寻找模式。如果你发现有三个、四个或其他一些数量的载荷存在一个显著的载荷模式，你甚至可能拥有一个新的个性心理学理论。
确定在因子分析中要使用的因素数量。在这个时候，一个自然的问题是要问：我们应该如何选择我们寻找的因素数量？最简单的答案是，我们可以使用 R 的psych包中的另一个命令，称为fa.parallel。我们可以用以下方式运行此命令：
```
parallel <- fa.parallel(big5, fm = 'minres', fa = 'fa')
```
再次，我们对函数背后的行为做出了选择。你可以尝试不同的fm和fa参数的选择，但你应该每次都看到实质上相似的结果。

fa.parallel命令的一个输出是以下斯克里普图：

Figure 5.10: Parallel analysis scree plot

img/C12628_05_11.jpg

图 5.10：平行分析斯克里普图

我们在第四章，降维中讨论了斯克里普图。标记为主成分特征值的 y 轴显示了每个因素在解释我们模型方差中的重要性。在大多数因子分析中，前几个因素具有相对较高的特征值，然后特征值急剧下降，并出现一个长长的平台期。这种急剧下降和长长的平台期共同形成了一个类似肘部的图像。因子分析中的常见做法是选择一个接近这种模式形成的肘部的因素数量。心理学家们共同选择了五个作为通常用于个性调查表的因素数量。你可以检查这个斯克里普图，看看你是否认为这是正确的数量，并且可以自由尝试使用不同数量的因素进行因子分析。

我们刚刚进行的练习的最终重要结果是斯克里普图。

你可能已经注意到，因子分析似乎与主成分分析有一些共同之处。特别是，两者都依赖于绘制特征值作为衡量不同向量重要性的方法之一的斯克里普图。

本书中对因子分析和主成分分析之间的全面比较超出了本书的范围。只需说明的是，它们具有相似的线性代数基础：两者都是为了创建协方差矩阵的近似，但每种方法实现这一点的途径略有不同。我们在因子分析中找到的因子与我们在主成分分析中找到的主成分类似。我们鼓励您尝试在相同的数据集上使用因子分析和主成分分析，并比较结果。结果通常非常相似，但并非完全相同。在大多数实际应用中，可以使用任何一种方法——您使用的具体方法将取决于您自己的偏好以及您对哪种方法最适合当前问题的判断。

活动 13：执行因子分析

在这个活动中，我们将使用因子分析对新的数据集进行分析。您可以在github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson05/Activity13找到数据集。

注意

此数据集来自 UCI 机器学习仓库。您可以在archive.ics.uci.edu/ml/machine-learning-databases/00484/tripadvisor_review.csv找到数据集。我们已经下载了文件并将其保存到github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson05/Activity13/factor.csv。

此数据集由 Renjith、Sreekumar 和 Jathavedan 编制，可在 UCI 机器学习仓库中找到。此数据集包含有关个人对旅游目的地所写评论的信息。每一行对应一个特定的唯一用户，总共有 980 个用户。列对应旅游景点的类别。例如，第二列记录每个用户的艺术画廊平均评分，第三列记录每个用户的舞厅平均评分。共有 10 个旅游景点类别。您可以在archive.ics.uci.edu/ml/datasets/Travel+Reviews找到关于每列记录的文档。

通过因子分析，我们将寻求确定不同类别用户评分之间的关系。例如，果汁吧（第 4 列）和餐厅（第 5 列）的用户评分可能相似，因为它们都由相同的潜在因素决定——用户对食物的兴趣。因子分析将尝试找到这些控制用户评分的潜在因素。我们可以遵循以下步骤进行因子分析：

下载数据并将其读入 R。
加载 psych 包。
选择记录用户评分的列子集。
为数据创建一个相关矩阵。
确定应使用的因子数量。
使用fa命令进行因子分析。
检查并解释因子分析的结果。

输出应类似于以下内容：

![图 5.11：因子分析预期结果

![img/C12628_05_12.jpg]

图 5.11：因子分析预期结果

注意

本活动的解决方案可以在第 231 页找到。

摘要

在本章中，我们讨论了与数据比较相关的话题。我们首先讨论了数据指纹的概念。为了说明数据指纹，我们引入了一个可以将任何字符串转换为固定范围内数字的哈希函数。这种哈希函数很有用，因为它使我们能够确保数据具有正确的身份——就像现实生活中的指纹一样。在介绍哈希函数之后，我们讨论了数据签名的重要性。数据签名或分析签名很有用，因为它使我们能够看到两个数据集是否是近似匹配的——指纹需要精确匹配。我们通过图像数据说明了分析签名的作用。我们通过介绍潜在变量模型来结束本章。我们讨论的潜在变量模型是因子分析。我们讨论了因子分析的一个应用，该应用使用心理调查数据来确定人们性格之间的差异。在下一章中，我们将讨论异常检测，这将使我们能够找到与数据集其他部分不匹配的观测值。

第七章：第六章

异常检测

学习目标

到本章结束时，你将能够：

使用参数和非参数方法在单变量和多变量数据中寻找异常值
使用数据转换来识别单变量和多变量数据中的异常值
使用马氏距离工作
通过结合季节性模型来提高异常检测性能

在本章中，我们将探讨不同的异常检测技术。

简介

数据分析通常隐含地假设所有观测值都是有效、准确和可靠的。但这种情况并不总是合理的假设。考虑信用卡公司的例子，他们收集的数据包括个人信用卡的消费记录。如果他们假设所有消费都是有效的，那么就会给盗贼和欺诈者可乘之机。相反，他们检查他们的交易数据集，寻找异常情况——即偏离一般观察模式的交易。由于欺诈交易没有标记，他们必须使用无监督学习来寻找这些异常情况，以防止犯罪活动。

在许多其他情况下，异常检测都是有用的。例如，制造商可能使用异常检测方法来寻找产品中的缺陷。医学研究人员可能在正常的心跳模式中寻找异常，以诊断疾病。IT 安全专业人员试图在服务器或计算机上找到异常活动，以识别恶意软件。在每种情况下，无监督学习方法都可以帮助区分有效观测值和异常观测值。

本章将介绍几种异常检测技术。我们将首先使用参数和非参数方法来寻找单变量和多变量数据中的异常值。然后，我们将讨论使用数据转换来识别单变量和多变量数据中的异常值。接下来，我们将探讨马氏距离，这是一种用于异常检测的多变量工具。我们将通过讨论使用回归模型来提高异常检测性能，通过结合季节性模型和检测上下文和集体异常来完成本章。

单变量异常值检测

在单变量情况下，异常检测是最简单的，也就是说，每个观测值只是一个数字。在这种情况下，我们可能首先通过检查观测值是否缺失、NULL、NA，或者记录为无穷大或其他与观测值类型不匹配的内容来进行常识性检查。完成此检查后，我们可以应用真正的无监督学习。

对于单变量数据，异常检测包括寻找异常值。R 的内置boxplot函数使得对异常值进行初步探索性检查变得非常容易，如下面的练习所示。

练习 37：使用 R 的 boxplot 函数进行异常值探索性视觉检查

对于单变量数据，异常检测包括寻找异常值。R 的内置箱线图函数使得在本练习中演示的初始探索性异常值检查变得非常简单。我们将使用一个名为mtcars的数据集，它是 R 内置的。

在这个练习中，我们将创建一个可以在图 6.3 中看到的箱线图。箱线图是一种重要的单变量图。在这里，围绕 3 的粗水平线表示数据中的中位数。围绕这个中位线的箱子在第一四分位数（25 百分位数）处有下限，在第三四分位数（75 百分位数）处有上限。垂直虚线延伸到所有非异常数据的上下端。这些虚线被称为图的触须，因此这种类型的图有时被称为“箱线和触须”图。最后，在图顶部附近表示为圆形点的两个观察值（至少根据 R 的结果）是异常值。

百分位数也可以称为分位数，在本例中我们将使用“分位数”来指代箱线图。分位数是数据中的一个点，它大于所有数据点中的一些固定比例。例如，0.1 分位数是数据中的一个观察值，它大于所有观察值的 10%，而小于其余部分。0.25 分位数也称为 25 百分位数或第一四分位数，它大于数据的 25%，而 75 百分位数或 0.75 分位数大于数据的 75%。当我们取一个观察值并找出它大于多少比例的观察值时，这被称为找到它的分位数。当我们取一个如 0.25 的分位数并试图找到一个与该分位数相对应的观察值时，我们可以称之为取逆分位数。

要使用 R 的boxplot函数进行异常值探索性视觉检查，请执行以下步骤：

要加载数据，打开 R 控制台并输入以下命令：
```
data(mtcars)
```
执行以下命令以查看数据集的前六行：
```
head(mtcars)
```
输出如下：

图 6.1：mtcars 数据集的前六行
您可以在 R 控制台中执行以下操作来查找有关此数据集的详细信息：
```
?mtcars
```
文档如下：

图 6.2：输出部分
按以下方式创建汽车重量的箱线图：
```
boxplot(mtcars$wt) 
```
输出将如下所示：

图 6.3：表示汽车重量的箱线图
看起来有两个观察值被 R 归类为异常值，它们的值似乎高于 5。请注意，这些重量是以千磅为单位的，所以实际上这些重量高于 5,000 磅。我们可以通过在我们的数据上运行一个简单的过滤器来确定这些观察值：
```
highest<-mtcars[which(mtcars$wt>5),]
print(highest)
```
输出如下所示：

图 6.4：重量超过 5,000 磅的汽车

我们可以看到有三个车型的重量超过 5,000 磅。由于我们只观察到两个异常值，我们得出结论，重量第三高的凯迪拉克 Fleetwood 不是异常值。这留下了其他两个：林肯大陆和克莱斯勒帝国，作为显然具有异常重量的车型。与其他车型相比，这些车型构成了异常。研究人员的一个潜在下一步是调查为什么这些车型似乎具有异常高的汽车重量。

在第三章，概率分布中，我们讨论了数据集倾向于遵循的不同分布。许多数据集都有所谓的长尾或胖尾，这意味着与平均值相差很大的观测值所占的比例不均衡——不一定是因为它们是异常的异常值，只是因为它们分布的性质。如果我们碰巧正在处理一个遵循胖尾分布的数据集，我们定义异常值的标准应该改变。

在以下练习中，我们将转换一个遵循胖尾分布的数据集，并观察哪些观测值被报告为异常值。

练习 38：将胖尾数据集转换为改进异常值分类

在以下练习中，我们将转换一个遵循胖尾分布的数据集，并观察哪些观测值被报告为异常值。我们将使用预加载到 R 中的rivers数据集：

按照以下方式加载数据集：
```
data(rivers)
```
执行以下命令以查看数据集的前六个观测值：
```
head(rivers)
```
输出如下：
```
[1] 735 320 325 392 524 450
```
你可以看到rivers数据集是一个向量。你可以通过输入以下内容来了解更多关于它的信息：
```
?rivers
```
文档说明如下：

图 6.5：河流数据集的信息
观察异常值的分布。首先，尝试通过运行以下命令来绘制河流数据的箱线图：
```
boxplot(rivers)
```
箱线图看起来是这样的：

图 6.6：河流数据集的箱线图

你可以看到这个箱线图与之前我们查看的mtcars箱线图看起来不同。特别是，箱线和须占据了较小的绘图部分，而大量的绘图由 R 分类为异常值的观测值组成。然而，异常值在任何数据集中都不应该非常多——根据定义，它们应该是稀少的。当我们观察到这样的箱线图，其中包含许多占据绘图大部分的异常值时，我们可以合理地得出结论，我们的分布是胖尾分布。
为了得到更好的异常值分类，对数据进行转换可能有所帮助。许多与自然和自然现象相关的数据集已知遵循对数正态分布。如果我们对每个观测值取对数，得到的分布将是正态分布（因此不是厚尾分布）。要转换数据，可以使用以下命令：
```
log_rivers<-log(rivers)
```
观察转换数据的箱线图。最后，我们可以执行以下命令来查看转换数据的箱线图：
```
boxplot(log_rivers)
```
箱线图将如下所示：

图 6.7：转换数据集的箱线图

如预期，箱线和须线占据了更大的绘图比例，被分类为异常值的观测值更少。在这种情况下，我们可以使用 log_rivers 图而不是 rivers 图来分类异常值。

之前的练习展示了数据准备在异常检测中的重要性。我们也可以尝试在原始数据集上进行异常检测。但有时，就像我们之前看到的 rivers 示例一样，通过进行一些简单的数据准备，我们可以得到不同且更好的结果。需要注意的是，数据有多种可能的转换方式。我们已经使用了对数转换，但还有许多其他可能有效的方法。另一件需要注意的事情是，数据准备既有利也有弊：一些数据准备方法，如平均和数据平滑，可能会导致我们丢弃有价值的信息，从而使我们的异常检测效果降低。

到目前为止，我们一直依赖于 R 的内置异常检测方法，并且我们已经通过简单的箱线图视觉检查来确定哪些观测值是异常值。在下一个练习中，我们将自己确定哪些观测值是异常值，而不使用 R 的内置异常分类。我们将找到数据的分位数——具体来说是 0.25 和 0.75 分位数。

练习 39：不使用 R 的内置 boxplot 函数找出异常值

在这个练习中，我们将自己确定哪些观测值是异常值，而不使用 R 的内置异常值分类。我们将找到数据的分位数——具体来说是 .25 和 .75 分位数。我们将再次使用 rivers 数据，这是我们之前练习中使用的数据：

通过执行以下命令加载数据：
```
data(rivers)
```
四分位距是第一四分位数（25 百分位数）和第三四分位数（75 百分位数）之间的差值。因此，我们可以通过运行以下命令将四分位距存储在名为 interquartile_range 的变量中：
```
interquartile_range<-unname(quantile(rivers,.75)-quantile(rivers,.25))
```
在这种情况下，我们使用 unname 确保将 interquartile_range 只是一个数字，而不是 DataFrame、列表或其他数据类型。这将使后续的工作更容易且更可靠。
使用以下命令检查四分位距：
```
print(interquartile_range)
```
输出如下：
```
[1] 370
```
R 中 boxplot 函数的标准方法是使用 1.5 倍的四分位距作为非异常观察值可以分散的范围的上限。然后非异常值的上限是第三四分位数加 1.5 倍的四分位距，非异常值的下限是第一四分位数减去 1.5 倍的四分位距。我们可以这样计算：
```
upper_limit<-unname(quantile(rivers,.75)+1.5*interquartile_range)
lower_limit<-unname(quantile(rivers,.25)-1.5*interquartile_range)
```
我们识别出的异常值是那些高于我们的upper_limit或低于我们的lower_limit的观察值。我们可以这样确定：
```
rivers[which(rivers>upper_limit | rivers<lower_limit)]
```
这将输出一个我们将其分类为异常值的观察值列表：
```
[1] 1459 1450 1243 2348 3710 2315 2533 1306 1270 1885 1770
```
注意

这个练习使用了 R 中 boxplot 函数所使用的方法。但在无监督学习中，关于方法细节的灵活性总是存在的。
另一种寻找异常值的方法是使用前一个练习中的方法，但使用不同的上限和下限值。我们可以通过改变乘以四分位距的系数来实现这一点，如下所示：
```
upper_limit<-unname(quantile(rivers,.75)+3*interquartile_range)
lower_limit<-unname(quantile(rivers,.25)-3*interquartile_range)
```
我们已经将系数从 1.5 改为 3。这种改变使得我们的方法不太可能将任何特定观察值分类为异常值，因为它增加了上限并降低了下限。

通常，你可以尝试以你认为合理且能导致良好结果的方式创造性地改变无监督学习方法。

前一个练习中的方法是所谓的非参数方法。在统计学中，有些方法被称为参数方法，而有些被称为非参数方法。参数方法对数据的潜在分布做出假设，例如，假设数据遵循正态分布。非参数方法旨在摆脱这些约束性假设。前一个练习中的方法仅依赖于分位数，因此它不对数据的分布或与之相关的参数（如均值和方差）做出任何假设。正因为如此，我们称之为非参数方法。相比之下，参数异常检测方法是一种对数据的分布或其参数（如均值和方差）做出假设的方法。请注意，非参数方法和参数方法总是在寻找相同的异常：不存在参数异常或非参数异常，只有参数方法和非参数方法。我们将在下一个练习中讨论参数异常检测方法。

我们将计算一个称为z 分数的东西。z 分数是衡量观察值与平均值之间距离的标准测量。每个 z 分数都是以标准差为单位测量的。因此，z 分数为 0 表示观察值与平均值相差 0 个标准差；换句话说，它等于平均值。z 分数为 1 表示观察值比平均值高 1 个标准差。z 分数为-2.5 表示观察值比平均值低 2.5 个标准差。在某些情况下，将距离平均值两个标准差以上的观察值视为异常或异常是惯例。

练习 40：使用参数方法检测异常值

在这个练习中，我们将通过寻找 z 分数大于 2 或小于-2 的观察值来调查异常。我们将使用前一个练习中使用的河流数据集：

按如下方式加载数据并确定标准差：
```
data(rivers)
standard_deviation<-sd(rivers)
```
通过计算每个观察值与平均值相差多少个标准差来确定每个观察值的 z 分数：
```
z_scores<-(rivers-mean(rivers))/ standard_deviation
```
通过选择 z 分数大于 2 或小于-2 的观察值来确定哪些观察值是异常值：
```
outliers<-rivers[which(z_scores>2 | z_scores<(-2))]
outliers
```
输出如下：
```
[1] 2348 3710 2315 2533 1885 1770
```
在这种情况下，我们可以看到有六条河流被归类为异常值——它们的 z 分数都高于 2。在这里，“异常值”这个词可能对这些异常来说过于强烈，因为 2 的 z 分数并不特别巨大。没有严格的规则来定义异常值，是否使用这个术语将取决于具体情况和背景。

参数异常检测，如本练习所示，是一种常见做法。然而，其适用性取决于其背后的参数假设是否有效。在这种情况下，我们计算了标准差，并寻找了距离平均值两个标准差以上的异常值。如果数据来自高斯（正态）分布，这种做法是合理的。在这种情况下，一些证据表明河流数据并不来自高斯分布。对于对基础数据分布有疑问的情况，可以使用非参数方法或转换方法（如前几项练习中使用的）可能更合适。

多变量异常检测

到目前为止的所有方法都集中在单变量异常值检测。然而，在实践中，这样的数据实际上是罕见的。大多数数据集包含关于数据多个属性的调查。在这些情况下，不清楚如何计算一个点是否是异常值。我们可以对每个单独的维度计算 z 分数，但对于那些在一个维度上 z 分数高而在另一个维度上低，或者在一个维度上相对较大而在另一个维度上平均的观察值，我们应该怎么办？没有简单的答案。在这些情况下，我们可以使用称为马氏距离的东西来计算 z 分数的多维类似物。马氏距离是 z 分数的多维类似物。这意味着它测量两点之间的距离，但使用特殊的单位，就像 z 分数一样，依赖于数据的方差。在完成以下练习后，我们将更好地理解马氏距离的含义。

练习 41：计算马氏距离

在这个练习中，我们将学习一种新的测量方法，它最终将帮助我们找到在考虑身高和体重时，与一般人群相比哪些个体是异常值。这个练习的预期输出将是一个单一数据点的马氏距离，测量它距离数据中心的距离。

注意

在后续的练习中，我们将计算数据集中所有点的马氏距离。然后，我们可以使用这些测量结果来分类哪些个体在考虑身高和体重时与一般人群相比是异常值：

注意

此数据集来自在线统计计算资源。您可以在wiki.stat.ucla.edu/socr/index.php/SOCR_Data_Dinov_020108_HeightsWeights找到数据集。我们已经下载了文件并将其保存在github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson06/Exercise41-42/heightsweights.csv。我们使用了数据的前 200 行。

首先，在 R 中加载此数据。首先，从我们的 GitHub 存储库github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson06/Exercise41-42/heightsweights.csv下载它。然后，将其保存在 R 的工作目录中。然后，按照以下方式将其加载到 R 中：
```
filename<-'heightsweights.csv'
raw<-read.csv(filename, stringsAsFactors=FALSE)
```
我们已经给数据集命名为raw。你可以确保列有正确的名称如下：
```
names(raw)<-c('index','height','weight')
```
绘制数据并观察模式：
```
plot(raw$height,raw$weight)
```
图形如下所示：

图 6.8：身高和体重的图

我们可以看到，在这个样本中，身高和体重都有很大的变化。请注意，我们的单变量方法在这个数据中不会直接起作用。我们可以计算身高的标准差或四分位数范围，并找到身高的异常值。但身高的异常值不一定是体重的异常值——它们可能是正常体重，或者正好符合体重的预期。同样，体重的异常值可能有完全平均或预期的身高。我们并不立即清楚如何计算“完全异常值”，或者说在某些整体意义上是异常的观测值，考虑到身高和体重两个方面。
接下来，我们将计算多维平均值的等价物，称为centroid：
```
centroid<-c(mean(raw$height),mean(raw$weight))
```
计算任何给定点和质心之间的距离。作为一个例子，我们将选择数据集中的第一个观测值：
```
example_distance<-raw[1,c('height','weight')]-centroid
```
计算我们数据中身高和体重协方差矩阵的逆。首先，我们计算身高和体重数据的协方差矩阵：
```
cov_mat<-cov(raw[,c('height','weight')])
```
使用 R 中的solve函数计算其逆，该函数用于计算矩阵逆：
```
inv_cov_mat<-solve(cov_mat)
```
计算我们的点和数据集质心之间的马氏距离：
```
mahalanobis_dist<-t(matrix(as.numeric(example_distance)))%*% matrix(inv_cov_mat,nrow=2) %*% matrix(as.numeric(example_distance))
```
在这个情况下，我们使用 %*% 因为它表示矩阵乘法，这正是我们想要执行的，我们同样需要将每个参数转换为一个数值矩阵。

这个练习的输出是一个马氏距离，它是多维 z 分数的推广：也就是说，它是衡量每个点距离均值多少个标准差的广义度量。在这种情况下，我们找到的马氏距离是 1.71672。马氏距离类似于任何类型的距离测量——0 是最低可能的测量值，数值越高，距离越远。只有质心会有测量的马氏距离为 0。马氏距离的优势在于它们是以一种考虑了每个变量的方差的方式进行标准化的，这使得它们在异常值检测中非常有效。在本章的后面部分，我们将看到如何使用这种类型的测量来找到这个多维数据集中的异常值。

在簇中检测异常

在前两章中，我们讨论了不同的聚类方法。我们现在讨论的方法，马氏距离，如果想象我们正在查看的数据是某个特定聚类的数据，就可以在聚类应用中非常有用。例如，在划分聚类中，具有最高马氏距离的点可以被选为从聚类中移除的点。此外，马氏距离的范围可以用来表示任何给定聚类的分散程度。

多变量异常值检测的其他方法

对于多变量异常值检测，还有其他方法，包括一些被称为非参数的方法。与非参数方法一样，一些先前的练习可能依赖于分位数，换句话说，从大到小排列每个观测值的排名，以分类异常值。一些非参数方法使用此类排名的总和来理解数据的分布。然而，这些方法并不常见，并且在一般情况下并不比马氏距离更有效，因此我们建议在多变量异常值检测中依赖马氏距离。

练习 42：基于马氏距离比较分类异常值

在这个练习中，我们将使用马氏距离的比较来分类异常值。这个练习是之前练习的延续，它将依赖于相同的数据库和一些相同的变量。在那个练习中，我们找到了一个数据点的马氏距离；现在我们将为所有数据点找到它。在执行以下练习之前，您应该运行之前练习中的所有代码，并确保您熟悉那里提出的思想：

注意

此数据集来自 UCI 机器学习仓库。您可以在wiki.stat.ucla.edu/socr/index.php/SOCR_Data_Dinov_020108_HeightsWeights找到数据集。我们已经下载了文件，并将其保存在github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson06/Exercise41-42/heightsweights.csv。

创建一个NULL变量，用于存储我们计算出的每个距离：
```
all_distances<-NULL
```

遍历每个观测值，并计算它与数据质心的马氏距离。此循环中的代码来自之前的练习，我们在那里学习了如何计算马氏距离：

k<-1
while(k<=nrow(raw)){
the_distance<-raw[k,c('height','weight')]-centroid
mahalanobis_dist<-t(matrix(as.numeric(the_distance)))%*% matrix(inv_cov_mat,nrow=2) %*% matrix(as.numeric(the_distance))
all_distances<-c(all_distances,mahalanobis_dist)
k<-k+1
}

运行此循环后，我们为数据中的每个点测量了马氏距离。

绘制所有具有特别高 Mahalanobis 距离的观测值。在这种情况下，我们将特别高定义为 Mahalanobis 距离的最高 10%。这意味着所有高于 Mahalanobis 距离的 0.9 分位数（我们将在以下代码中选择）的 Mahalanobis 距离：
```
plot(raw$height,raw$weight)
points(raw$height[which(all_distances>quantile(all_distances,.9))],raw$weight[which(all_distances>quantile(all_distances,.9))],col='red',pch=19)
```
输出如下：

图 6.9：具有高 Mahalanobis 距离的观测值

此图显示了 Mahalanobis 距离最高的 10% 的每个点的实心红色点。我们可以看到，其中一些在身高和体重方面似乎都是异常值，并且可以通过执行我们的单变量异常检测方法来观察到。然而，图表上许多红色点既不是身高的单变量异常，也不是体重的单变量异常。它们只是相对于整个点云的异常，而 Mahalanobis 距离使我们能够量化并检测这一点。

在季节性数据中检测异常

到目前为止，我们只讨论了异常检测作为一种检测异常的方法。然而，异常检测不仅仅包括异常检测。一些异常不容易被检测为原始异常。接下来，我们将查看季节性、趋势性数据。在这类数据中，我们希望找到在季节性趋势或周期背景下发生的异常。

我们将使用 R 中的 expsmooth 包的数据。

我们将使用的数据记录了 1985 年至 2005 年间澳大利亚每月游客数量，以千人计。

在以下练习中，我们将处理具有时间趋势的数据。通过时间趋势，我们指的是随着时间的推移，观测值倾向于增加或减少（在这种情况下，它们倾向于增加）。为了检测异常，我们想要做的是去趋势。去趋势数据意味着尽可能去除其时间趋势，以便我们可以找到每个观测值与预期值的偏差程度。

练习 43：执行季节性建模

在这个练习中，我们将尝试对数据进行建模，以确定我们应该将哪些视为数据的预期值，哪些视为与预期值的偏差。预期的输出是一组误差值，我们将在未来的练习中使用这些误差值来分类异常——误差最大的观测值将被分类为数据集的异常：

首先，通过执行以下命令将此数据加载到 R 中：
```
install.packages("expsmooth")
library(expsmooth)
data(visitors)
plot(visitors)
```
输出如下：

图 6.10：1985 年至 2005 年间澳大利亚每月游客数量的图表

按以下方式检查数据的第一个六个观测值：

head(visitors)

输出如下：

       May   Jun   Jul   Aug   Sep   Oct
1985  75.7  75.4  83.1  82.9  77.3 105.7

按以下方式检查最后六个观测值：

tail(visitors)

输出如下：

       Jan   Feb   Mar   Apr May Jun Jul Aug Sep Oct   Nov   Dec
2004                                                 479.9 593.1
2005 462.4 501.6 504.7 409.5

由于日期可能难以处理，我们可以分配一个数值变量来跟踪日期的顺序。我们将把这个变量称为 period 并定义如下：
```
period<-1:length(visitors)
```
将visitors数据与刚刚创建的period变量结合，将它们都放入一个名为raw的 DataFrame 中：
```
raw<-data.frame(cbind(visitors,period))
```
这一步在技术上属于监督学习，而不是无监督学习。在这种情况下，我们将使用监督学习，但仅作为进行无监督学习过程中的一个中间步骤。为了找到数据中的时间趋势，我们可以运行一个线性回归，将销售额与数值时间周期相关联，如下所示：
```
timetrend<-lm(visitors~period+I(log(period)),data=raw)
```
接下来，我们可以获得这个时间趋势的拟合值，并将它们作为raw DataFrame 的一部分存储：
```
raw$timetrend<-predict(timetrend,raw)
```
去趋势的过程意味着我们从第 6 步中找到的预测趋势中减去。我们这样做的原因是我们想找到数据中的异常值。如果我们保留数据中的趋势，看起来像异常的东西实际上可能是趋势的预期结果。通过从数据中去除趋势，我们可以确保观察到的异常不是趋势的结果。我们可以按照以下方式完成去趋势：
```
raw$withouttimetrend<-raw$visitors-raw$timetrend
```
我们可以绘制去趋势数据的一个简单图表如下：
```
plot(raw$withouttimetrend,type='o')
```
绘图如下：

图 6.11：去趋势数据的绘图

在这个图表中，你应该注意到数据中没有明显的从左到右的趋势，这表明我们已经成功地去除了趋势。

我们的数据记录了澳大利亚每月的访客数量。澳大利亚的季节温度和天气有变化，因此合理地假设游客访问量可能会根据特定月份的天气是否宜人而增加或减少。甚至可能存在与季节和全年天气变化相关的商务或外交访问澳大利亚的变化。因此，合理地假设澳大利亚的访问量存在年度模式。我们的下一步将是将这些年度模式从我们的数据中去除。
首先，我们创建一个矩阵，其中每一列包含关于一年中不同月份的去趋势访客数据：
```
seasonsmatrix = t(matrix(data = raw$withouttimetrend, nrow = 12))
```
计算每一列的平均值，以得到该特定月份去趋势访客的平均值：
```
seasons = colMeans(seasonsmatrix, na.rm = T)
```
这给我们一个包含 12 个值的向量——每年一个月。由于我们有 20 年的数据，我们将这个向量重复 20 次：
```
raw$seasons<-c(rep(seasons,20))
```
最后，我们可以获得去趋势、去周期的数据，我们将它命名为error，因为在移除了数据中的时间趋势和年度周期后，随机误差是唯一剩下的东西：
```
raw$error<-raw$visitors-raw$timetrend-raw$seasons
```
我们可以绘制这个图表来查看其外观：
```
plot(raw$error,type='o')
```
绘图如下：

图 6.12：去趋势数据的绘图
将季节性建模的所有元素一起绘制。最后，我们可以展示我们从季节性数据中分离出的所有元素：时间趋势、年度周期和随机误差：
```
par(mfrow=c(3,1))
plot(raw$timetrend,type='o')
plot(raw$seasons,type='o')
plot(raw$error,type='o')
```
输出如下：

图 6.13：季节性建模元素的绘图

图 6.13 中一起显示的三个图显示了原始数据的分解。第一个图显示了时间趋势，换句话说，就是数据中观察到的整体模式，即每个月的访问量往往比上个月多。当然，这是对数据的过度简化，因为也存在季节性模式：某些月份的访问量往往比其他月份多或少，无论整体趋势如何。第二个图显示了这些季节性模式，这些模式每年都以相同的方式重复。最后，第三个图显示了误差，这是数据中未被整体时间趋势或每年内的季节性模式捕获的所有变化。如果你将这三个图中呈现的所有数据相加，你将恢复原始数据集。但在这些三个图中，我们可以看到这些元素被分解，这种分解使我们更好地理解了时间趋势、季节性和误差如何相互作用，从而构成观察到的数据。

这个练习使我们能够创建一个名为 error 的变量，并将其添加到原始数据框中。现在我们已经创建了误差向量，我们可以使用标准的异常检测来在数据框中找到异常。

练习 44：使用参数方法在季节性数据中寻找异常

在这个练习中，我们将通过寻找与预期值的最大偏差来进行异常检测：

计算前一个练习中计算出的误差数据的标准差：
```
stdev<-sd(raw$error)
```

找出哪些数据点距离平均值超过两个标准差：

high_outliers<-which(raw$error>(mean(raw$error)+2*sd(raw$error)))
low_outliers<-which(raw$error<(mean(raw$error)-2*sd(raw$error)))

检查我们已分类为异常的观察结果：

raw[high_outliers,]

输出如下：

    visitors period timetrend withouttimetrend   seasons    error
1       75.7      1  67.18931         8.510688 -55.94655 64.45724
130    392.7    130 305.93840        86.761602  24.35847 62.40313
142    408.0    142 323.44067        84.559332  24.35847 60.20086
147    397.4    147 330.70509        66.694909  11.55558 55.13933
188    559.9    188 389.78579       170.114205  91.11673 78.99748

低异常被分类如下：

raw[low_outliers,]

输出如下：

    visitors period timetrend withouttimetrend   seasons      error
80     266.8     80  231.4934         35.30663  91.11673  -55.81010
216    321.5    216  429.7569       -108.25691 -20.46137  -87.79553
217    260.9    217  431.1801       -170.28007 -55.94655 -114.33352
218    308.3    218  432.6029       -124.30295 -42.40371  -81.8992

我们可以如下绘制这些点：

plot(raw$period,raw$visitors,type='o')
points(raw$period[high_outliers],raw$visitors[high_outliers],pch=19,col='red')
points(raw$period[low_outliers],raw$visitors[low_outliers],pch=19,col='blue')

图形如下所示：

图 6.14：被分类为异常的数据的图形

该图显示了所有我们的数据和我们已经分类为异常的点，高异常用红色表示，低异常用蓝色表示。你应该注意，并非所有这些点都立即明显是异常值。我们只是在之前进行的季节性建模练习中，确定了异常偏离的预期值后，才将它们识别为异常。

接下来，我们将介绍两种更多类型的异常：上下文异常和集体异常。为了介绍这些概念，我们将生成一个包含上下文和集体异常的人工数据集。

上下文和集体异常

在图 6.15 中大约 x=1 的位置，你可以看到一个与邻居相当远的单独点。这是一个上下文异常的例子，我们将首先讨论这意味着什么以及如何检测这些类型的异常。在大约 x=3.6 的位置，你可以看到一个 y 值平坦的区域，每个值都等于 0。在这个数据中，零值并不异常，但这么多零值放在一起就是异常。因此，这些数据被称作集体异常。

我们将首先考虑上下文异常。上下文异常是指仅因为其邻居而被认为是异常的观测值。在我们刚刚生成的数据集中，有一个点在 x=1 处，y=0。然而，x=0.99 和 x=1.01 处的 y 值接近 0.84，在这个上下文中与 0 相差甚远。上下文异常可以通过找到与邻居距离异常的观测值来检测，正如我们将在以下练习中看到的那样。

练习 45：检测上下文异常

以下练习展示了如何检测我们刚刚介绍的数据集中的上下文异常。由于上下文异常是与邻居非常不同的观测值，我们需要对每个观测值与其邻居进行显式比较。为了做到这一点，我们计算第一差分。第一差分简单地是观测值的值减去其前一个观测值的值。

本练习的预期结果将是数据集中那些具有上下文异常的观测值：

我们将生成一个包含上下文和集体异常的人工数据集。该数据集可以通过在 R 控制台中运行以下代码生成：
```
x<-1:round(2*pi*100+100)/100
y<-rep(0,round(2*pi*100)+100)
y[1:314]<-sin(x[1:314])
y[415:728]<-sin(x[315:628])
y[100]<-0
```
你可以通过以下方式绘制x和y来查看这些数据的外观：
```
plot(x,y,type='o')
```
输出如下：

图 6.15：生成数据集的绘图

此图显示了一个正弦曲线：一种从 0 开始，向上和向下轻轻倾斜，最终再次回到 0 的简单曲线。这些数据是人工生成的；然而，我们可以想象它代表了温度的观测值：在某些月份较低，在某些月份上升较高，在某些月份较高，在其他月份下降。温度和天气数据通常遵循可以用正弦或余弦等三角曲线建模的规律。我们已经修改了我们的正弦曲线，使其包含一些异常数据。
同时为每个观测值找到第一差分如下：
```
difference_y<-y[2:length(y)]-y[1:(length(y)-1)]
```
创建第一差分数据的箱线图：
```
boxplot(difference_y)
```
结果箱线图如下所示：

图 6.16：第一差分数据的箱线图

此箱线图显示，几乎所有第一差分都非常接近零，而两个异常差分则远离其他数据。我们可以从第一差分箱线图中看到，单个高异常值大于 0.5。

你可能会注意到图 6.16 中有两个明显的异常值，但在图 6.15 中只有一个明显的异常值。这是因为图 6.16 显示的是一阶差分数据，数据中的单个异常值导致了两个大的第一阶差异：第 99 个和第 100 个值之间的差异，以及第 100 个和第 101 个值之间的差异。原始数据中的一个异常观察值导致了一阶差分数据中的两个异常观察值。
使用 R 的有用 which 函数确定哪个观察值对应于这个异常值：
```
which(difference_y>0.5)
```
which 返回值为 100，表示它是第 100 个观察值，是一个上下文异常。如果我们检查，y[100] 等于 0，而其邻居不是，所以我们已经成功找到了上下文异常。

接下来，我们将讨论集体异常。在我们的 x-y 图中，我们指出了在 x=3.64 处所有观察值都等于 0 的 100 个观察点。在这个数据集中，0 值本身并不是异常，但是有 100 个观察值都等于 0 则是异常的。这里 100 个零值一起被称为集体异常。集体异常比上下文异常更难检测，但我们将尝试在下面的练习中检测它们。

练习 46：检测集体异常

在接下来的练习中，我们将检测我们之前创建的数据集中的集体异常。这个练习的预期结果是构成集体异常的观察值列表：

为了检测这种异常，我们需要寻找包含没有变化或只有微小变化的观察值组或邻域。以下循环实现了这一点。它创建了一个向量，该向量由两个差异的最大值组成：一个观察值与其 50 个周期前的观察值之间的差异，以及一个观察值与其 50 个周期后的观察值之间的差异。如果这些差异中的最大值是零或非常小，那么我们就检测到了一个极端平坦的邻域，这是这种集体异常的迹象。以下是我们将使用的循环：
```
changes<-NULL
ks<-NULL
k<-51
while(k<(length(y)-50)){
changes<-c(changes,max(abs(y[k-50]),abs(y[k+50])))
ks<-c(ks,k)
k<-k+1
}
```
这个循环创建了一个名为 changes 的向量。这个向量的每个元素都测量了一个观察值与其邻居之间观察到的最大差异，距离为 50 个观察点。changes 向量中特别小的值将表明我们可能有一个由平坦邻域组成的集体异常。
现在我们有一个测量邻域变化的向量，我们可以创建一个简单的箱线图，就像我们在之前的练习中所做的那样：
```
boxplot(changes)
```
输出如下：

图 6.17：邻域变化的箱线图

我们可以看到有很多被分类为异常值的观察值，这表明存在非常低的邻域变化。
我们可以按照以下方式找到导致这种集体异常的观察值：
```
print(ks[which(changes==min(changes))])
```
输出为 364。
我们可以通过检查该索引处的 y 值来验证这是对应于集体异常的索引：
```
print(y[ks[which(changes==min(changes))]])
```
因此，y 的值为 0，这是集体异常处的 y 值，这为我们找到了正确的点提供了证据。

核密度

为了结束本章，我们将讨论如何使用核密度估计在一系列血液样本上执行异常值检测。核密度估计提供了一种自然的方式来测试特定的血液检测结果是否异常，即使没有特定血液测试或医学的一般专业知识。

假设你在一个诊所工作，你的老板要求你对患者进行一种新的血液测试。你的老板想知道是否有任何患者的测试结果异常。然而，你对这种新的血液测试不熟悉，也不知道正常和异常结果应该是什么样子。你所拥有的只是老板向你保证来自正常患者的先前血液测试记录。假设这些测试的结果如下：

normal_results<-c(100,95,106,92,109,190,210,201,198)

现在假设你的老板要求你找出以下新血液测试结果中的异常值（如果有的话）：

new_results<-c(98,35,270,140,200)

在核密度估计中，我们使用一组核来建模数据集。对于我们的目的，核将只是具有不同均值和方差的正态分布。你可以在第三章，概率分布中了解更多关于正态分布的信息。我们将假设我们的数据具有由正态分布之和捕获的密度。然后，任何看起来与我们指定的正态分布之和不一致的数据都可以被分类为异常。

当我们计算核密度时，我们必须指定一个叫做带宽的东西。在这里，带宽将是我们用来建模数据的正态分布的方差的度量。如果我们指定一个高的带宽，我们假设数据广泛分散，如果我们指定一个低的带宽，我们假设数据主要包含在一个相对较窄的范围内。随着你完成以下练习，这应该会变得更加清晰。

练习 47：使用核密度估计寻找异常值

在这个练习中，我们将介绍如何使用核密度估计寻找异常值。这个练习的预期输出是按照核密度估计方法构成异常的一组观测值列表：

指定我们的数据和参数。对于我们的数据，我们将使用之前指定的正常结果和新结果：
```
normal_results<-c(100,95,106,92,109,190,210,201,198)
new_results<-c(98,35,270,140,200)
```
核密度估计依赖于一个称为带宽的参数。因此，我们将从 25 开始。带宽的选择将取决于您的个人偏好以及原始数据。如果您不确定，可以选择一个与您数据的标准差大致相同的带宽。在这种情况下，我们将选择 25，这低于我们数据的标准差。您可以如下设置带宽为 25：
```
bandwidth<-25
```
注意

如果您喜欢，可以自由地尝试其他带宽值。

使用 R 的density函数来获得核密度估计：

our_estimate<-density(normal_results, bw=bandwidth)

绘制密度估计：
```
plot(our_estimate)
```
生成的图表如下所示：

图 6.18：密度估计图
图 6.18 中的图形形状代表我们原始数据的分布。您可以在第三章，概率分布中了解更多关于概率分布的信息。这个分布被称为双峰分布，这意味着数据似乎主要围绕两个点聚集：大多数观察值接近 100 或 200。如果观察值在图 6.18 所示的分布中的对应点显示它特别不可能，我们将将其解释为异常。
我们可以为我们的每个新结果获得密度估计。对于new_results中的每个观察，我们将根据图 6.18 中展示的核来计算密度估计。我们将按照以下方式将这些密度估计存储在新变量中：
```
new_density_1<-density(normal_results,bw=25,n=1,from=new_results[1],to=new_results[1])$y
new_density_2<-density(normal_results,bw=25,n=1,from=new_results[2],to=new_results[2])$y
new_density_3<-density(normal_results,bw=25,n=1,from=new_results[3],to=new_results[3])$y
new_density_4<-density(normal_results,bw=25,n=1,from=new_results[4],to=new_results[4])$y
new_density_5<-density(normal_results,bw=25,n=1,from=new_results[5],to=new_results[5])$y
```
输出将是图表上从步骤 3对应于new_results向量中指定的每个 x 值的的高度。我们可以通过打印来观察这些值。按照以下方式打印new_density_1：
```
print(new_density_1)
```
输出如下：
```
[1] 0.00854745
```
按照以下方式打印new_density_2：
```
print(new_density_2)
```
输出如下：
```
[1] 0.0003474778
```
按照以下方式打印new_density_3：
```
print(new_density_3)
```
输出如下：
```
[1] 0.0001787185
```
按照以下方式打印new_density_4：
```
print(new_density_4)
```
输出如下：
```
[1] 0.003143966
```
按照以下方式打印new_density_5：
```
print(new_density_5)
```
输出如下：
```
[1] 0.006817359
```

我们可以将这些点绘制在我们的密度图上如下：

plot(our_estimate)
points(new_results[1],new_density_1,col='red',pch=19)
points(new_results[2],new_density_2,col='red',pch=19)
points(new_results[3],new_density_3,col='red',pch=19)
points(new_results[4],new_density_4,col='red',pch=19)
points(new_results[5],new_density_5,col='red',pch=19)

执行这些命令的结果如下所示：![图 6.19：密度图上的点映射

图 6.19：密度图上的点映射

这显示了之前我们检查过的相同密度图，增加了五个点——每个新结果对应一个，我们正在检查这些新结果以寻找异常值。每个点都显示了每个特定观察的相对可能性：具有更高估计密度值的点更有可能被观察到，而具有较低估计密度值的点不太可能被观察到，因此更有可能是异常值。关于哪些观察是异常值，哪些不是，没有严格的规则，但一般来说，估计密度值最接近零的观察最可能是异常值。
解释结果以分类异常。

每个密度值看起来都相当接近零。然而，有些比其他更接近零。在这种情况下，密度估计值越接近零，我们越有信心认为观察结果是异常的。在这种情况下，我们将选择一个阈值值，并说那些核密度估计值低于阈值的血液检测结果是异常的，而那些核密度估计值高于阈值的检测结果不是异常的。看起来 0.001 是一个合理的阈值，因为它将高密度值与最低密度值分开——密度值低于 0.001 的观察结果在步骤 6所示的图表中看起来非常不可能。因此，我们将血液检测结果 35 和 270 归类为异常结果，而其他所有结果都视为合理，因为我们看到在步骤 4 中，35 和 270 对应的是低于 0.001 的密度估计值。

因此，我们练习的最终结论是血液检测结果 35 和 270 是异常的，而其他所有血液检测结果都是合理的。

继续学习异常检测

如果你继续学习异常检测，你会发现有大量的不同异常检测技术。然而，所有这些技术都遵循我们在季节性建模示例中看到的基本模式。具体来说，高级异常检测通常包括以下内容：

指定预期的模型
计算基于模型预期的值与观察到的值之间的差异——这被称为误差
使用单变量异常值检测对误差向量进行检测以确定异常

最大的困难在于第一步：指定一个有用的预期模型。在这种情况下，我们指定了一个季节性模型。在其他情况下，将需要指定考虑多维图像数据、音频记录、复杂的经济指标或其他复杂属性的模型。为这些情况设置模型的方法将需要研究数据来源的特定领域。然而，在每种情况下，异常检测都将遵循前面列出的三个步骤的要点模式。

活动十四：使用参数方法和非参数方法寻找单变量异常

活动的目的是使用参数方法来寻找单变量异常。为此活动，我们将使用 R 中内置的数据集，称为islands。如果你在 R 中执行?islands，你可以找到这个数据集的文档。在这个文档中，你可以注意到这个数据集包含了地球上面积超过 10,000 平方英里的陆地面积。

这可能是一个对研究地球陆地形成过程的地质学家感兴趣的数据集。根据科学家的说法，岛屿可以通过多种方式形成：有时是通过火山活动，有时是通过珊瑚生长，有时是通过其他方式。地质学家可能对寻找异常大或小的岛屿感兴趣——这些岛屿可能是进行进一步研究以尝试理解岛屿自然形成过程的最佳地点。在这个活动中，我们将寻找islands数据集中的异常值。

这些步骤将帮助您完成活动：

在 R 的内置数据集中加载名为islands的数据，并创建该数据的箱线图。您注意到数据分布和异常值有什么特点？
对islands数据进行对数变换，并创建变换后数据的箱线图。这如何改变了被分类为异常值的数据点？
使用非参数方法（该方法将异常值定义为位于第一四分位数和第三四分位数以上或以下 1.5 倍四分位距的点）手动计算islands数据集中的异常值。对islands数据的对数变换也进行同样的操作。
使用参数方法通过计算数据的均值和标准差来对islands数据集中的异常值进行分类，将异常值分类为距离均值超过两个标准差的观测值。对islands数据的对数变换也进行同样的操作。
比较这些异常值检测方法的结果。

注意

本活动的解决方案可在第 234 页找到。

活动十五：使用马氏距离寻找异常值

在接下来的活动中，我们将检查与汽车速度和停车距离相关的数据。这些数据可能对试图研究哪些汽车表现最佳的汽车工程师有所帮助。那些相对于其速度具有特别低停车距离的汽车可以用作高性能汽车的例子，而那些相对于其速度具有异常高停车距离的汽车可能是进一步研究以寻找改进领域的候选者。在这个活动中，我们将基于速度和停车距离来寻找异常值。因为我们正在处理多变量数据，所以使用马氏距离来寻找异常值是有意义的。

这些步骤将帮助您完成这项活动：

从 R 的内置数据集中加载cars数据集。这个数据集包含一些非常古老的汽车的速度以及在该速度下所需的停车距离。绘制数据图表。
计算这些数据的质心，并计算每个点与质心之间的马氏距离。找出异常值（与质心的马氏距离最高的点）。绘制一个显示所有观测值和异常值的图表。

图表将如下所示：

图 6.20：标记了异常值的图表

注意

本活动的解决方案在第 237 页。

摘要

在本章中，我们讨论了异常检测。我们首先介绍了单变量异常检测，包括非参数和参数方法。我们讨论了进行数据转换以获得更好的异常分类。然后，我们讨论了使用马氏距离的多变量异常检测。我们完成了更多高级练习，以分类与季节性变化数据相关的异常。我们讨论了集体和上下文异常，并在本章结束时讨论了如何在异常检测中使用核密度估计。

第八章：附录

关于

本节包含的内容旨在帮助学生完成书中的活动。它包括学生为实现活动目标需要执行的详细步骤。

第一章：聚类方法简介

活动一：具有三个聚类的 k-means 聚类

解答：

在 iris_data 变量中加载 Iris 数据集：
```
iris_data<-iris
```
创建一个 t_color 列并将其默认值设为 red。将两种物种的值改为 green 和 blue，这样第三个就保持为 red：
```
iris_data$t_color='red'
iris_data$t_color[which(iris_data$Species=='setosa')]<-'green'
iris_data$t_color[which(iris_data$Species=='virginica')]<-'blue'
```
注意

在这里，我们只更改那些物种为 setosa 或 virginica 的值的 color 列）
选择任意三个随机的聚类中心：
```
k1<-c(7,3)
k2<-c(5,3)
k3<-c(6,2.5)
```

通过在 plot() 函数中输入花萼长度和花萼宽度以及颜色来绘制 x，y 图：

plot(iris_data$Sepal.Length,iris_data$Sepal.Width,col=iris_data$t_color)
points(k1[1],k1[2],pch=4)
points(k2[1],k2[2],pch=5)
points(k3[1],k3[2],pch=6)

这里是输出：

图 1.36：给定聚类中心的散点图

选择迭代次数：
```
number_of_steps<-10
```
选择 n 的初始值：
```
n<-1
```
开始 while 循环以找到聚类中心：
```
while(n<number_of_steps){
```

计算每个点到当前聚类中心的距离。这里我们使用 sqrt 函数计算欧几里得距离：

iris_data$distance_to_clust1 <- sqrt((iris_data$Sepal.Length-k1[1])²+(iris_data$Sepal.Width-k1[2])²)
iris_data$distance_to_clust2 <- sqrt((iris_data$Sepal.Length-k2[1])²+(iris_data$Sepal.Width-k2[2])²)
iris_data$distance_to_clust3 <- sqrt((iris_data$Sepal.Length-k3[1])²+(iris_data$Sepal.Width-k3[2])²)

将每个点分配到最近的聚类中心：

  iris_data$clust_1 <- 1*(iris_data$distance_to_clust1<=iris_data$distance_to_clust2 & iris_data$distance_to_clust1<=iris_data$distance_to_clust3)
  iris_data$clust_2 <- 1*(iris_data$distance_to_clust1>iris_data$distance_to_clust2 & iris_data$distance_to_clust3>iris_data$distance_to_clust2)
  iris_data$clust_3 <- 1*(iris_data$distance_to_clust3<iris_data$distance_to_clust1 & iris_data$distance_to_clust3<iris_data$distance_to_clust2)

通过使用 R 中的 mean() 函数计算每个中心的平均 x 和 y 坐标来计算新的聚类中心：

  k1[1]<-mean(iris_data$Sepal.Length[which(iris_data$clust_1==1)])
  k1[2]<-mean(iris_data$Sepal.Width[which(iris_data$clust_1==1)])
  k2[1]<-mean(iris_data$Sepal.Length[which(iris_data$clust_2==1)])
  k2[2]<-mean(iris_data$Sepal.Width[which(iris_data$clust_2==1)])
  k3[1]<-mean(iris_data$Sepal.Length[which(iris_data$clust_3==1)])
  k3[2]<-mean(iris_data$Sepal.Width[which(iris_data$clust_3==1)])
  n=n+1
}

为每个中心选择颜色以绘制散点图：

iris_data$color='red'
iris_data$color[which(iris_data$clust_2==1)]<-'blue'
iris_data$color[which(iris_data$clust_3==1)]<-'green'

绘制最终的图表：

plot(iris_data$Sepal.Length,iris_data$Sepal.Width,col=iris_data$color)
points(k1[1],k1[2],pch=4)
points(k2[1],k2[2],pch=5)
points(k3[1],k3[2],pch=6)

输出如下：

图 1.37：用不同颜色表示不同物种的散点图

活动二：使用 k-means 进行客户细分

解答：

从 github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson01/Activity02/wholesale_customers_data.csv 下载数据。

将数据读入 ws 变量：

ws<-read.csv('wholesale_customers_data.csv')

只存储 ws 变量中的第 5 和第 6 列，丢弃其他列：
```
ws<-ws[5:6]
```
导入 factoextra 库：
```
library(factoextra)
```
计算两个中心的聚类中心：
```
clus<-kmeans(ws,2)
```
绘制两个聚类的图表：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.38：两个聚类的图表

注意到异常值也成为了两个聚类的一部分。
计算三个聚类的聚类中心：
```
clus<-kmeans(ws,3)
```
绘制三个聚类的图表：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.39：三个聚类的图表

注意现在一些异常值已经成为一个单独聚类的部分。
计算四个中心的聚类中心：
```
clus<-kmeans(ws,4)
```
绘制四个聚类的图表：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.40：四个聚类的图表

注意到异常值已经开始分离成两个不同的聚类。
计算五个聚类的聚类中心：
```
clus<-kmeans(ws,5)
```
绘制五个聚类的图表：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.41：五个聚类的图表

注意到异常值如何明显地形成了两个分开的红色和蓝色聚类，而其余数据被分类到三个不同的聚类中。
计算六个聚类的聚类中心：
```
clus<-kmeans(ws,6)
```
绘制六个聚类的图表：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.42：六个聚类的图表

活动 3：使用 k-medoids 聚类进行客户细分

解决方案：

将 CSV 文件读入ws变量：

ws<-read.csv('wholesale_customers_data.csv')

只在ws变量中存储第 5 和第 6 列：
```
ws<-ws[5:6]
```
导入factoextra库进行可视化：
```
library(factoextra)
```
导入cluster库进行 PAM 聚类：
```
library(cluster)
```
通过在pam函数中输入数据和聚类数量来计算聚类：
```
clus<-pam(ws,4)
```
绘制聚类可视化图：
```
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.43：聚类 K-medoid 图
再次，使用 k-means 聚类计算聚类并绘制输出，以与pam聚类的输出进行比较：
```
clus<-kmeans(ws,4)
fviz_cluster(clus,data=ws)
```
输出如下：

图 1.44：聚类 K-means 图

活动 4：寻找理想的市场细分数量

解决方案：

将下载的数据集读入ws变量：

ws<-read.csv('wholesale_customers_data.csv')

通过丢弃其他列，只在变量中存储第 5 和第 6 列：
```
ws<-ws[5:6]
```
使用轮廓分数计算最佳聚类数量：
```
fviz_nbclust(ws, kmeans, method = "silhouette",k.max=20)
```
这里是输出：

图 1.45：表示轮廓分数最佳聚类数量的图表

根据轮廓分数，最佳聚类数量是两个。
使用 WSS 分数计算最佳聚类数量：
```
fviz_nbclust(ws, kmeans, method = "wss", k.max=20)
```
这里是输出：

图 1.46：WSS 分数下的最佳聚类数量

根据 WSS 肘部方法，最佳聚类数量大约是六个。
使用 Gap 统计量计算最佳聚类数量：
```
fviz_nbclust(ws, kmeans, method = "gap_stat",k.max=20)
```
这里是输出：

图 1.47：Gap 统计量下的最佳聚类数量

根据 Gap 统计量，最佳聚类数量是一个。

第二章：高级聚类方法

活动 5：在蘑菇数据集上实现 k-modes 聚类

解决方案：

从 github.com/TrainingByP… 下载mushrooms.csv。
下载后，在 R 中加载mushrooms.csv文件：
```
ms<-read.csv('mushrooms.csv')
```
检查数据集的维度：
```
dim(ms)
```
输出如下：
```
[1] 8124   23
```
检查所有列的分布：
```
summary.data.frame(ms)
```
输出如下：

图 2.29：所有列分布摘要的屏幕截图

每一列包含所有唯一的标签及其计数。
将数据集的所有列（除了最终标签）存储在一个新变量ms_k中：
```
ms_k<-ms[,2:23]
```
导入具有kmodes函数的klaR库：
```
install.packages('klaR')
library(klaR)
```
计算kmodes聚类并将它们存储在kmodes_ms变量中。将不带true标签的数据集作为第一个参数输入，将聚类数量作为第二个参数输入：
```
kmodes_ms<-kmodes(ms_k,2)
```

通过创建一个包含true标签和cluster标签的表格来检查结果：

result = table(ms$class, kmodes_ms$cluster)
result

输出如下：

       1    2
  e   80 4128
  p 3052  864

如你所见，大多数可食用蘑菇都在第 2 个聚类中，而大多数有毒蘑菇都在第 1 个聚类中。因此，使用 k-modes 聚类已经合理地完成了识别每种蘑菇是否可食用或有毒的工作。

活动 6：实现 DBSCAN 并可视化结果

解决方案：

导入dbscan和factoextra库：
```
library(dbscan)
library(factoextra)
```
导入multishapes数据集：
```
data(multishapes)
```
将multishapes数据集的列放入ms变量中：
```
ms<-multishapes[,1:2]
```
按以下方式绘制数据集：
```
plot(ms)
```
输出如下：

图 2.30：多形状数据集的折线图
对数据集执行 k-means 聚类并绘制结果：
```
km.res<-kmeans(ms,4)
fviz_cluster(km.res, ms,ellipse = FALSE)
```
输出如下：

图 2.31：多形状数据集上的 k-means 折线图

对ms变量执行 DBSCAN 并绘制结果：

db.res<-dbscan(ms,eps = .15)
fviz_cluster(db.res, ms,ellipse = FALSE,geom = 'point')

输出如下：

图 2.32：多形状数据集上的 DBCAN 折线图

在这里，你可以看到所有黑色点都是异常值，并且不在任何聚类中，而 DBSCAN 形成的聚类在其他任何聚类方法中都是不可能的。这些聚类具有所有类型的形状和大小，而 k-means 中的所有聚类都是球形的。

活动 7：对种子数据集进行层次聚类分析

解决方案：

将下载的文件读入sd变量：
```
sd<-read.delim('seeds_dataset.txt')
```
注意

根据系统上文件的位置修改路径。
首先，将数据集的所有列（除了最终标签）放入sd_c变量中：
```
sd_c<-sd[,1:7]
```
导入cluster库：
```
library(cluster)
```
计算层次聚类并绘制树状图：
```
h.res<-hclust(dist(sd_c),"ave")
plot(h.res)
```
输出如下：

图 2.33：聚类树状图
在k=3处切割树并绘制一个表格，以查看聚类结果在分类三种种子类型方面的表现：
```
memb <- cutree(h.res, k = 3)
results<-table(sd$X1,memb)
results
```
输出如下：

图 2.34：分类三种种子类型的表格
对sd_c数据集执行划分聚类并绘制树状图：
```
d.res<-diana(sd_c,metric ="euclidean",)
plot(d.res)
```
输出如下：

图 2.35：分裂聚类的树状图
在k=3处切割树，并绘制一个表格以查看聚类结果在分类三种种子类型时的表现：
```
memb <- cutree(h.res, k = 3)
results<-table(sd$X1,memb)
results
```
输出如下：

图 2.36：分类三种种子类型的表格

您可以看到两种聚类方法都产生了相同的结果。这些结果还表明，分裂聚类是层次聚类的逆过程。

第三章：概率分布

活动 8：寻找与鸢尾花数据集中变量分布最接近的标准分布

解答：

将df变量加载到 Iris 数据集中：
```
df<-iris
```
仅选择对应 setosa 物种的行：
```
df=df[df$Species=='setosa',]
```
导入kdensity库：
```
library(kdensity)
```
使用kdensity函数计算并绘制花瓣长度的核密度估计图：
```
dist <- kdensity(df$Sepal.Length)
plot(dist)
```
输出如下：

图 3.36：花瓣长度的核密度估计图

这个分布与我们之前章节中研究的正态分布最接近。在这里，平均值和中位数都在 5 左右。
使用kdensity函数从df中计算并绘制花瓣宽度的核密度估计图：
```
dist <- kdensity(df$Sepal.Width)
plot(dist)
```
输出如下：

图 3.37：花瓣宽度的核密度估计图

这个分布也最接近正态分布。我们可以通过柯尔莫哥洛夫-斯米尔诺夫检验来形式化这种相似性。

活动 9：使用正态分布计算累积分布函数和执行柯尔莫哥洛夫-西蒙诺夫检验

解答：

将 Iris 数据集加载到df变量中：
```
df<-iris
```
仅保留 setosa 物种的行：
```
df=df[df$Species=='setosa',]
```

计算df中花瓣长度列的平均值和标准差：

sdev<-sd(df$Sepal.Length)
mn<-mean(df$Sepal.Length)

使用花瓣长度列的标准差和平均值生成一个新的分布：
```
xnorm<-rnorm(100,mean=mn,sd=sdev)
```
绘制xnorm和花瓣长度列的累积分布函数：
```
plot(ecdf(xnorm),col='blue')
plot(ecdf(df$Sepal.Length),add=TRUE,pch = 4,col='red')
```
输出如下：

图 3.38：xnorm 和花瓣长度的累积分布函数

在分布中，样本看起来非常接近。让我们看看，在下一个测试中，花瓣长度样本是否属于正态分布。
对两个样本进行柯尔莫哥洛夫-斯米尔诺夫检验，如下所示：
```
ks.test(xnorm,df$Sepal.Length)
```
输出如下：
```
    Two-sample Kolmogorov-Smirnov test
data: xnorm and df$Sepal.Length
D = 0.14, p-value = 0.5307
alternative hypothesis: two-sided
```
这里，p-value非常高，而D值很低，因此我们可以假设花瓣长度的分布接近正态分布。

对df中花瓣宽度列重复相同的步骤：

sdev<-sd(df$Sepal.Width)
mn<-mean(df$Sepal.Width)
xnorm<-rnorm(100,mean=mn,sd=sdev)
plot(ecdf(xnorm),col='blue')
plot(ecdf(df$Sepal.Width),add=TRUE,pch = 4,col='red')

输出如下：

图 3.39：xnorm 和花瓣宽度的累积分布函数

按如下进行柯尔莫哥洛夫-斯米尔诺夫检验：

ks.test(xnorm,df$Sepal.Length)

输出如下：

    Two-sample Kolmogorov-Smirnov test
data: xnorm and df$Sepal.Width
D = 0.12, p-value = 0.7232
alternative hypothesis: two-sided

这里，花瓣宽度的样本分布也接近正态分布。

第四章：降维

活动 10：在新数据集上执行 PCA 和市场篮子分析

解答：

在开始我们的主要分析之前，我们将移除一个对我们无关的变量：
```
Boston<-Boston[,-12]
```
我们将创建虚拟变量。最终我们将得到一个原始数据集和一个虚拟变量数据集。我们这样做如下：
```
Boston_original<-Boston
```
接下来，我们将为原始数据集中的每个测量值创建虚拟变量。您可以在 MASS 包的文档中找到每个变量的含义，该文档可在cran.r-project.org/web/packages/MASS/MASS.pdf找到。

为是否一个城镇的人均犯罪率高低创建虚拟变量：

Boston$highcrim<-1*(Boston$indus>median(Boston$crim))
Boston$lowcrim<-1*(Boston$indus<=median(Boston$crim))

为是否一个城镇划定的地块面积超过 25,000 平方英尺的比例高或低创建虚拟变量：

Boston$highzn<-1*(Boston$zn>median(Boston$zn))
Boston$lowzn<-1*(Boston$zn<=median(Boston$zn))

为是否一个城镇非零售商业用地比例高或低创建虚拟变量：

Boston$highindus<-1*(Boston$indus>median(Boston$indus))
Boston$lowindus<-1*(Boston$indus<=median(Boston$indus))

为是否一个城镇毗邻查尔斯河创建虚拟变量：

Boston$highchas<-(Boston$chas)
Boston$lowchas<-(1-Boston$chas)

为是否一个城镇氮氧化物浓度高或低创建虚拟变量：

Boston$highnox<-1*(Boston$nox>median(Boston$nox))
Boston$lownox<-1*(Boston$nox<=median(Boston$nox))

为是否一个城镇每户平均房间数量高或低创建虚拟变量：

Boston$highrm<-1*(Boston$rm>median(Boston$rm))
Boston$lowrm<-1*(Boston$rm<=median(Boston$rm))

为是否一个城镇拥有 1940 年之前建造的业主自住单元比例高或低创建虚拟变量：

Boston$highage<-1*(Boston$age>median(Boston$age))
Boston$lowage<-1*(Boston$age<=median(Boston$age))

为是否一个城镇距离波士顿五个就业中心平均距离高或低创建虚拟变量：

Boston$highdis<-1*(Boston$dis>median(Boston$dis))
Boston$lowdis<-1*(Boston$dis<=median(Boston$dis))

为是否一个城镇对放射状高速公路的可达性指数高或低创建虚拟变量：

Boston$highrad<-1*(Boston$rad>median(Boston$rad))
Boston$lowrad<-1*(Boston$rad<=median(Boston$rad))

为是否一个城镇的全值财产税税率高或低创建虚拟变量：

Boston$hightax<-1*(Boston$tax>median(Boston$tax))
Boston$lowtax<-1*(Boston$tax<=median(Boston$tax))

为是否一个城镇的学生-教师比例高或低创建虚拟变量：

Boston$highptratio<-1*(Boston$ptratio>median(Boston$ptratio))
Boston$lowptratio<-1*(Boston$ptratio<=median(Boston$ptratio))

为是否一个城镇低阶层人口比例高或低创建虚拟变量：

Boston$highlstat<-1*(Boston$lstat>median(Boston$lstat))
Boston$lowlstat<-1*(Boston$lstat<=median(Boston$lstat))

为是否一个城镇的中位家庭价值高或低创建虚拟变量：

Boston$highmedv<-1*(Boston$medv>median(Boston$medv))
Boston$lowmedv<-1*(Boston$medv<=median(Boston$medv))

创建一个完全由我们刚刚创建的虚拟变量组成的数据集：
```
Bostondummy<-Boston[,14:ncol(Boston)]
```
最后，我们将我们的Boston_2数据集恢复到添加所有虚拟变量之前的原始形式：
```
Boston<-Boston_original
```
按如下方式计算数据集协方差矩阵的特征值和特征向量：
```
Boston_cov<-cov(Boston)
Boston_eigen<-eigen(Boston_cov)
print(Boston_eigen$vectors)
```
输出如下：

图 4.17：协方差矩阵的特征向量
按如下方式打印特征值：
```
print(Boston_eigen$values)
```
输出如下：

图 4.18：协方差矩阵的特征值
对于第三部分，我们将基于特征值创建一个简单的斯克奇图：
```
plot(Boston_eigen$values,type='o')
```
输出如下：

图 4.19：特征值的绘图
接下来，我们选择我们将使用的特征向量的数量（我选择了 10），并将数据集转换为 10 维，如下所示：
```
neigen<-10
transformed<-t(t(as.matrix(Boston_eigen$vectors[,1:neigen])) %*% t(as.matrix(Boston)))
```

然后，我们将尽可能恢复数据集：

restored<- t(as.matrix(Boston_eigen$vectors[,1:neigen]) %*% t(as.matrix(transformed)))

最后，我们可以检查我们的恢复与原始数据集的接近程度，如下所示：
```
print(head(restored-Boston))
```
在这里，我们需要指定一个support阈值（例如，20%），并完成对数据的第一次遍历：
```
support_thresh<-0.2
firstpass<-unname(which(colMeans(Bostondummy,na.rm=TRUE)>support_thresh))
```

在这里，我们完成对数据的第二次遍历：

secondcand<-t(combn(firstpass,2))
secondpass<-NULL
k<-1
while(k<=nrow(secondcand)){
support<-mean(Bostondummy[,secondcand[k,1]]*Bostondummy[,secondcand[k,2]],na.rm=TRUE)
if(support>support_thresh){
secondpass<-rbind(secondpass,secondcand[k,])
}
k<-k+1
}

在这里，我们完成第三次遍历，然后根据confidence和lift阈值进行过滤：

thirdpass<-NULL
k<-1
while(k<=nrow(secondpass)){
j<-1
while(j<=length(firstpass)){
n<-1
product<-1
while(n<=ncol(secondpass)){
product<-product*Bostondummy[,secondpass[k,n]]
n<-n+1
}
if(!(firstpass[j] %in% secondpass[k,])){
product<-product*Bostondummy[,firstpass[j]]
support<-mean(product,na.rm=TRUE)
if(support>support_thresh){
thirdpass<-rbind(thirdpass,c(secondpass[k,],firstpass[j]))
}
}
j<-j+1
}
k<-k+1
}
thirdpass_conf<-NULL
k<-1
while(k<=nrow(thirdpass)){
support<-mean(Bostondummy[,thirdpass[k,1]]*Bostondummy[,thirdpass[k,2]]*Bostondummy[,thirdpass[k,3]],na.rm=TRUE)
confidence<-mean(Bostondummy[,thirdpass[k,1]]*Bostondummy[,thirdpass[k,2]]*Bostondummy[,thirdpass[k,3]],na.rm=TRUE)/mean(Bostondummy[,thirdpass[k,1]]*Bostondummy[,thirdpass[k,2]],na.rm=TRUE)
lift<-confidence/mean(Bostondummy[,thirdpass[k,3]],na.rm=TRUE)
thirdpass_conf<-rbind(thirdpass_conf,unname(c(thirdpass[k,],support,confidence,lift)))
k<-k+1
}

我们最终的输出是那些通过了support、confidence和lift阈值的三个项目篮子的列表：
```
print(head(thirdpass_conf))
```
输出结果如下：

图 4.20：三项篮子的输出

第五章：数据比较方法

活动 11：为人物照片创建图像签名

解决方案：

将 Borges 照片下载到你的电脑上，并保存为borges.jpg。确保它保存在 R 的工作目录中。如果它不在 R 的工作目录中，那么使用setwd()函数更改 R 的工作目录。然后，你可以将这张图片加载到一个名为im（代表图像）的变量中，如下所示：
```
install.packages('imager')
library('imager')
filepath<-'borges.jpg'
im <- imager::load.image(file =filepath) 
```
我们将要探索的其余代码将使用这张图片，称为im。在这里，我们已经将阿拉莫的图片加载到im中。然而，你可以通过将图片保存到你的工作目录并在filepath变量中指定其路径来运行其余的代码。
我们正在开发的签名旨在用于灰度图像。因此，我们将使用imager包中的函数将此图像转换为灰度：
```
im<-imager::rm.alpha(im)
im<-imager::grayscale(im)
im<-imager::imsplit(im,axis = "x", nb = 10)   
```
代码的第二行是转换为灰度。最后一行执行了将图像分割成 10 等份的操作。
以下代码创建了一个空矩阵，我们将用关于我们 10x10 网格每个部分的详细信息来填充它：
```
matrix <- matrix(nrow = 10, ncol = 10)
```
接下来，我们将运行以下循环。这个循环的第一行使用了imsplit命令。这个命令之前也被用来将 x 轴分成 10 等份。这次，对于 x 轴的每个 10 等份，我们将沿着 y 轴进行分割，也将它分成 10 等份：
```
for (i in 1:10) {
  is <- imager::imsplit(im = im[[i]], axis = "y", nb = 10)
  for (j in 1:10) {
    matrix[j,i] <- mean(is[[j]])
  }
}
```
到目前为止的输出是matrix变量。我们将在步骤 4中使用它。
通过运行以下代码获取 Borges 照片的签名：
```
borges_signature<-get_signature(matrix)
borges_signature
```
输出结果如下：

图 5.12：borges_signature 矩阵
接下来，我们将开始使用 9x9 矩阵而不是 10x10 矩阵来计算一个签名。我们开始使用之前使用过的相同过程。以下代码行加载我们的 Borges 图像，就像我们之前做的那样。代码的最后一行将图像分割成相等的部分，但这次我们设置nb=9，以便将图像分割成 9 等份：
```
filepath<-'borges.jpg'
im <- imager::load.image(file =filepath) 
im<-imager::rm.alpha(im)
im<-imager::grayscale(im)
im<-imager::imsplit(im,axis = "x", nb = 9)
```
以下代码创建了一个空矩阵，我们将用关于我们 9x9 网格每个部分的详细信息来填充它：
```
matrix <- matrix(nrow = 9, ncol = 9)
```
注意我们使用nrow=9和ncol=9，这样我们就有了一个 9x9 的矩阵来填充我们的亮度测量值。
接下来，我们将运行以下循环。此循环的第一行使用 imsplit 命令。此命令之前也用于将 x 轴分割成 9 个相等的部分。这次，对于每个 9 个 x 轴分割，我们将沿着 y 轴进行分割，也将它分割成 9 个相等的部分：
```
for (i in 1:9) {
  is <- imager::imsplit(im = im[[i]], axis = "y", nb = 9)
  for (j in 1:9) {
    matrix[j,i] <- mean(is[[j]])
  }
}
```
到目前为止的输出是 matrix 变量。我们将重复 步骤 4。

通过运行以下代码获取博尔赫斯照片的 9x9 签名：

borges_signature_ninebynine<-get_signature(matrix)
borges_signature_ninebynine

输出如下：

图 5.13：borges_signature_ninebynine 矩阵

活动 12：为水印图像创建图像签名

解决方案：

将带水印的图片下载到您的计算机上，并将其保存为 alamo_marked.jpg。请确保它保存在 R 的工作目录中。如果它不在 R 的工作目录中，请使用 setwd() 函数更改 R 的工作目录。然后，您可以将此图像加载到名为 im（代表图像）的变量中，如下所示：
```
install.packages('imager')
library('imager')
filepath<-'alamo_marked.jpg'
im <- imager::load.image(file =filepath) 
```
我们将要探索的其余代码将使用名为 im 的此图像。在这里，我们已经将阿拉莫的水印图片加载到 im 中。然而，您只需将图片保存到您的工作目录中，并在 filepath 变量中指定其路径，就可以在任意图片上运行其余的代码。
我们正在开发的签名旨在用于灰度图像。因此，我们将使用 imager 包中的函数将此图像转换为灰度：
```
im<-imager::rm.alpha(im)
im<-imager::grayscale(im)
im<-imager::imsplit(im,axis = "x", nb = 10)   
```
此代码的第二行是转换为灰度。最后一行执行将图像分割成 10 个相等部分的操作。
以下代码创建了一个空矩阵，我们将用有关我们 10x10 网格每个部分的信息来填充它：
```
matrix <- matrix(nrow = 10, ncol = 10)
```
接下来，我们将运行以下循环。此循环的第一行使用 imsplit 命令。此命令之前也用于将 x 轴分割成 10 个相等的部分。这次，对于每个 10 个 x 轴分割，我们将沿着 y 轴进行分割，也将它分割成 10 个相等的部分：
```
for (i in 1:10) {
  is <- imager::imsplit(im = im[[i]], axis = "y", nb = 10)
  for (j in 1:10) {
    matrix[j,i] <- mean(is[[j]])
  }
}
```
到目前为止的输出是 matrix 变量。我们将在 步骤 4 中使用此变量。
我们可以通过运行以下代码获取水印照片的签名：
```
watermarked_signature<-get_signature(matrix)
watermarked_signature
```
输出如下：

图 5.14：水印图像的签名

此活动的最终输出是 watermarked_signature 变量，它是水印阿拉莫照片的分析签名。如果您已经完成了迄今为止的所有练习和活动，那么您应该有三个分析签名：一个名为 building_signature，一个名为 borges_signature，还有一个名为 watermarked_signature。
完成此活动后，我们将此签名存储在名为 watermarked_signature 的变量中。现在，我们可以将其与我们的原始阿拉莫签名进行比较，如下所示：
```
comparison<-mean(abs(watermarked_signature-building_signature))
comparison
```
在这种情况下，我们得到的结果是 0.015，这表明原始图像签名与这个新图像的签名非常接近。

我们所看到的是，我们的分析签名方法对相似图像返回相似的签名，对不同的图像返回不同的签名。这正是我们希望签名所做的，因此我们可以判断这种方法是成功的。

活动 13：执行因子分析

解决方案：

数据文件可以从github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-R/tree/master/Lesson05/Data/factor.csv下载。将其保存到您的计算机上，并确保它在 R 的工作目录中。如果您将其保存为factor.csv，则可以通过执行以下命令在 R 中加载它：
```
factor<-read.csv('factor.csv')
```
按如下方式加载psych包：
```
library(psych)
```
我们将对用户评分进行因子分析，这些评分记录在数据文件的 2 至 11 列中。我们可以如下选择这些列：
```
ratings<-factor[,2:11]
```
按如下方式创建评分数据的相关矩阵：
```
ratings_cor<-cor(ratings)
```
通过创建灰度图来确定我们应该使用多少个因素。灰度图是以下命令的输出之一：
```
parallel <- fa.parallel(ratings_cor, fm = 'minres', fa = 'fa')
```
灰度图如下所示：

图 5.15：平行分析灰度图

灰度图显示了其中一个特征值远高于其他特征的特征。虽然我们在分析中可以自由选择任意数量的特征，但其中一个特征值远大于其他特征，这为我们使用一个特征进行分析提供了充分的理由。
我们可以如下进行因子分析，指定nfactors参数中的因子数量：
```
factor_analysis<-fa(ratings_cor, nfactors=1)
```
这将我们的因子分析结果存储在一个名为factor_analysis的变量中：
我们可以如下检查我们的因子分析结果：
```
print(factor_analysis)
```
输出如下所示：

图 5.16：因子分析结果

在MR1下的数字显示了我们对每个类别的单特征因子载荷。由于我们只有一个解释性因素，所有在这个因素上有正载荷的类别之间都是正相关。我们可以将这个因素解释为普遍的积极性，因为它会表明，如果人们对一个类别评价很高，他们也会对其他类别评价很高；如果他们对一个类别评价很低，他们很可能也会对其他类别评价很低。

唯一的例外是“类别 10”，它记录了用户对宗教机构的平均评分。在这种情况下，因子载荷大且为负。这表明，那些对大多数其他类别评分很高的人往往对宗教机构评分较低，反之亦然。因此，我们可以将我们发现的积极因子解释为对休闲活动的积极态度，而不是对宗教机构的积极态度，因为宗教机构可能不是休闲场所，而是崇拜场所。似乎在这个数据集中，对休闲活动持积极态度的人对崇拜持消极态度，反之亦然。对于接近 0 的因子载荷，我们也可以得出关于休闲活动积极性的规则不太适用的结论。你可以看到，因子分析使我们能够找到我们数据中观测值之间的关系，这是我们之前未曾怀疑的。

第六章：异常检测

活动 14：使用参数方法和非参数方法寻找单变量异常值

解答：

按以下方式加载数据：
```
data(islands)
```
按以下方式绘制箱线图：
```
boxplot(islands)
```
图 6.21：岛屿数据集的箱线图

你应该注意到数据具有极端的厚尾分布，这意味着中位数和四分位距在图中所占比例相对较小，与 R 软件分类为异常值的许多观测值相比。
按以下方式创建一个新的对数转换后的数据集：
```
log_islands<-log(islands)
```
按以下方式创建对数转换数据的箱线图：
```
boxplot(log_islands)
```
图 6.22：对数转换后的数据集箱线图

你应该注意到对数转换后只有五个异常值。

计算四分位距：

interquartile_range<-quantile(islands,.75)-quantile(islands,.25)

将四分位距的 1.5 倍加到第三四分位数上，以得到非异常数据的上限：
```
upper_limit<-quantile(islands,.75)+1.5*interquartile_range
```
将异常值定义为任何高于此上限的观测值：
```
outliers<-islands[which(islands>upper_limit)]
```

计算对数转换数据的四分位距：

interquartile_range_log<-quantile(log_islands,.75)-quantile(log_islands,.25)

将四分位距的 1.5 倍加到第三四分位数上，以得到非异常数据的上限：
```
upper_limit_log<-quantile(log_islands,.75)+1.5*interquartile_range_log
```

将异常值定义为任何高于此上限的观测值：

outliers_log<-islands[which(log_islands>upper_limit_log)]

按以下方式打印未转换的异常值：
```
print(outliers)
```
对于未转换的异常值，我们得到以下结果：

图 6.23：未转换的异常值

按以下方式打印对数转换后的异常值：
```
print(outliers_log)
```
对于对数转换后的异常值，我们得到以下结果：

图 6.24：对数转换后的异常值

计算数据的均值和标准差：

island_mean<-mean(islands)
island_sd<-sd(islands)

选择距离均值超过两个标准差的观测值：
```
outliers<-islands[which(islands>(island_mean+2*island_sd))]
outliers
```
我们得到以下异常值：

图 6.25：异常值的截图

首先，我们计算对数转换数据的均值和标准差：

island_mean_log<-mean(log_islands)
island_sd_log<-sd(log_islands)

选择距离平均值超过两个标准差的观测值：

outliers_log<-log_islands[which(log_islands>(island_mean_log+2*island_sd_log))]

我们如下打印对数变换后的异常值：
```
print(outliers_log)
```
输出如下：

图 6.26：对数变换后的异常值

活动 15：使用马氏距离查找异常值

解答：

您可以按照以下步骤加载数据并绘制图表：
```
data(cars)
plot(cars)
```
输出图表如下：

图 6.27：汽车数据集的绘图

计算质心：

centroid<-c(mean(cars$speed),mean(cars$dist))

计算协方差矩阵：
```
cov_mat<-cov(cars)
```
计算协方差矩阵的逆：
```
inv_cov_mat<-solve(cov_mat)
```
创建一个NULL变量，它将保存我们计算出的每个距离：
```
all_distances<-NULL
```

我们可以遍历每个观测值，并计算它们与数据集质心的马氏距离：

k<-1
while(k<=nrow(cars)){
the_distance<-cars[k,]-centroid
mahalanobis_dist<-t(matrix(as.numeric(the_distance)))%*% matrix(inv_cov_mat,nrow=2) %*% matrix(as.numeric(the_distance))
all_distances<-c(all_distances,mahalanobis_dist)
k<-k+1
}

绘制具有特别高马氏距离的所有观测值，以查看我们的异常值：

plot(cars)
points(cars$speed[which(all_distances>quantile(all_distances,.9))], cars$dist[which(all_distances>quantile(all_distances,.9))],col='red',pch=19)

我们可以如下查看输出图表，异常点用红色表示：

图 6.28：标记异常值的绘图

R-应用无监督学习-二-