R-机器学习秘籍-四-

布客飞龙
83 阅读39分钟

R 机器学习秘籍（四）

原文：annas-archive.org/md5/167951e4c9873bc214756c25a63518c8

译者：飞龙

协议：CC BY-NC-SA 4.0

第十章. 关联分析和序列挖掘

在本章中，我们将涵盖以下主题：

  • 将数据转换为交易

  • 显示交易和关联

  • 使用 Apriori 规则挖掘关联

  • 剪枝冗余规则

  • 可视化关联规则

  • 使用 Eclat 算法挖掘频繁项集

  • 创建具有时间信息的交易

  • 使用 cSPADE 挖掘频繁顺序模式

简介

企业从日常运营中积累大量交易数据（例如，来自零售商的销售订单、发票和运输文件），发现数据中的隐藏关系可能很有用，例如，“哪些产品经常一起购买？”或“购买手机后的后续购买是什么？”为了回答这两个问题，我们需要在事务数据集上执行关联分析和频繁顺序模式挖掘。

关联分析是一种在事务数据集中寻找有趣关系的方法。一个著名的关于产品的关联是购买尿布的客户也会购买啤酒。虽然这种关联听起来可能有些不寻常，但如果零售商能够利用这类信息或规则向客户交叉销售产品，那么他们提高销售额的可能性很高。

关联分析用于寻找项集之间的相关性，但如果你想知道频繁购买的项目顺序，该怎么办？为了实现这一点，你可以采用频繁顺序模式挖掘，从具有时间信息的事务数据集中找到频繁子序列。然后，你可以使用挖掘出的频繁子序列来预测客户购物顺序、网页点击流、生物序列以及其他应用中的使用情况。

在本章中，我们将涵盖创建和检查事务数据集的配方，使用 Apriori 算法执行关联分析，以各种图形格式可视化关联，以及使用 Eclat 算法查找频繁项集。最后，我们将创建具有时间信息的交易，并使用 cSPADE 算法发现频繁顺序模式。

将数据转换为交易

在创建挖掘关联规则之前，您需要将数据转换为交易。在下面的配方中，我们将介绍如何将列表、矩阵或数据框转换为交易。

准备工作

在这个配方中，我们将在一个列表、矩阵或数据框中生成三个不同的数据集。然后，我们可以将生成的数据集转换为交易。

如何操作...

执行以下步骤将不同格式的数据转换为交易：

  1. 首先，您必须安装并加载arule包：

    > install.packages("arules")
> library(arules)

  2. 您可以创建一个包含三个向量的购买记录列表：

    > tr_list = list(c("Apple", "Bread", "Cake"),
+                c("Apple", "Bread", "Milk"),
+                c("Bread", "Cake", "Milk"))
> names(tr_list) = paste("Tr",c(1:3), sep = "")

  3. 接下来，您可以使用as函数将数据框转换为交易：

    > trans = as(tr_list, "transactions")
> trans
transactions in sparse format with
 3 transactions (rows) and
 4 items (columns)

  4. 您还可以将矩阵格式的数据转换为交易：

    > tr_matrix = matrix(
+   c(1,1,1,0,
+     1,1,0,1,
+     0,1,1,1), ncol = 4)
> dimnames(tr_matrix) =  list(
+   paste("Tr",c(1:3), sep = ""),
+   c("Apple","Bread","Cake", "Milk")
+   )
> trans2 =  as(tr_matrix, "transactions")
> trans2
transactions in sparse format with
 3 transactions (rows) and
 4 items (columns)

  5. 最后，你可以将数据框格式的数据集转换为事务：

    > Tr_df = data.frame(
+   TrID= as.factor(c(1,2,1,1,2,3,2,3,2,3)),
+   Item = as.factor(c("Apple","Milk","Cake","Bread",
+                      "Cake","Milk","Apple","Cake",
+                      "Bread","Bread")) 
+ )
> trans3 = as(split(Tr_df[,"Item"], Tr_df[,"TrID"]), "transactions")
> trans3
transactions in sparse format with
 3 transactions (rows) and
 4 items (columns)

它是如何工作的...

在挖掘频繁项集或使用关联规则之前，通过事务类别准备数据集非常重要。在本菜谱中，我们演示如何将数据集从列表、矩阵和数据框格式转换为事务。在第一步中，我们生成一个包含三个购买记录向量的列表格式的数据集。然后，在为每个交易分配交易 ID 之后，我们使用 as 函数将数据转换为事务。

接下来，我们将演示如何将矩阵格式中的数据转换为事务。为了表示项目的购买情况，应使用二元发生矩阵来记录每个交易相对于不同购买项目的购买行为。同样，我们可以使用 as 函数将数据集从矩阵格式转换为事务。

最后，我们将展示如何将数据集从数据框格式转换为事务。数据框包含两个因子类型的向量：一个是名为 TrID 的交易 ID，另一个显示不同交易中的购买项目（命名为 Item）。此外，可以使用 as 函数将数据框格式的数据转换为事务。

参考内容

  • transactions 类用于表示规则或频繁模式挖掘的事务数据。它是 itemMatrix 类的扩展。如果您对如何使用这两个不同的类来表示事务数据感兴趣，请使用 help 函数参考以下文档：

    > help(transactions)
> help(itemMatrix)

显示事务和关联

arule 包使用其自己的 transactions 类来存储事务数据。因此，我们必须使用 arule 提供的通用函数来显示事务和关联规则。在本菜谱中，我们将展示如何通过 arule 包中的各种函数显示事务和关联规则。

准备工作

确保您已通过生成事务并将这些存储在变量 trans 中完成前面的菜谱。

如何做到这一点...

执行以下步骤以显示事务和关联：

  1. 首先，您可以获取事务数据的 LIST 表示：

    > LIST(trans)
$Tr1
[1] "Apple" "Bread" "Cake" 

$Tr2
[1] "Apple" "Bread" "Milk" 

$Tr3
[1] "Bread" "Cake"  "Milk"

  2. 接下来，您可以使用 summary 函数显示事务的统计和详细摘要：

    > summary(trans)
transactions as itemMatrix in sparse format with
 3 rows (elements/itemsets/transactions) and
 4 columns (items) and a density of 0.75 

most frequent items:
 Bread   Apple    Cake    Milk (Other) 
 3       2       2       2       0 

element (itemset/transaction) length distribution:
sizes
3 
3 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 3       3       3       3       3       3 

includes extended item information - examples:
 labels
1  Apple
2  Bread
3   Cake

includes extended transaction information - examples:
 transactionID
1           Tr1
2           Tr2
3           Tr3

  3. 然后，您可以使用 inspect 函数显示事务：

    > inspect(trans)
 items   transactionID
1 {Apple, 
 Bread, 
 Cake}            Tr1
2 {Apple, 
 Bread, 
 Milk}            Tr2
3 {Bread, 
 Cake, 
 Milk}            Tr3

  4. 此外，您还可以根据大小筛选事务：

    > filter_trains = trans[size(trans) >=3]
> inspect(filter_trains)
 items   transactionID
1 {Apple, 
 Bread, 
 Cake}            Tr1
2 {Apple, 
 Bread, 
 Milk}            Tr2
3 {Bread, 
 Cake, 
 Milk}            Tr3

  5. 此外，您还可以使用图像函数来直观检查事务：

    > image(trans)

    如何做到这一点...

    事务的可视检查

  6. 为了直观地显示频率/支持条形图，可以使用 itemFrequenctPlot：

    > itemFrequencyPlot (trans)

    如何做到这一点...

    事务项目频率条形图

它是如何工作的...

由于事务数据是挖掘关联和频繁模式的基石，我们必须学习如何显示关联以获得洞察力并确定关联是如何构建的。arules包提供了各种方法来检查事务。首先，我们使用LIST函数获取事务数据的列表表示。然后，我们可以使用summary函数获取信息，例如基本描述、最频繁的项以及事务长度分布。

接下来，我们使用inspect函数显示事务。除了显示所有事务外，您还可以首先通过大小过滤事务，然后使用inspect函数显示关联。此外，我们可以使用image函数对事务进行视觉检查。最后，我们说明如何使用频率/支持条形图显示每个项目的相对频率。

相关内容

  • 除了使用itemFrequencyPlot显示频率/条形图外，您还可以使用itemFrequency函数显示支持分布。有关更多详细信息，请使用help函数查看以下文档：

    > help(itemFrequency)

使用 Apriori 规则挖掘关联

关联挖掘是一种技术，可以挖掘出隐藏在事务数据集中的有趣关系。这种方法首先找到所有频繁项集，然后从频繁项集中生成强关联规则。Apriori 是最著名的关联挖掘算法，它首先识别频繁的个体项，然后执行广度优先搜索策略，将个体项扩展到更大的项集，直到无法找到更大的频繁项集。在本食谱中，我们将介绍如何使用 Apriori 规则进行关联分析。

准备工作

在本食谱中，我们将使用内置的事务数据集Groceries，演示如何在arules包中使用 Apriori 算法进行关联分析。请确保首先安装并加载arules包。

如何操作...

执行以下步骤以分析关联规则：

  1. 首先，您需要加载数据集Groceries：

    > data(Groceries)

  2. 您可以检查Groceries数据集的摘要：

    > summary(Groceries)

  3. 接下来，您可以使用itemFrequencyPlot检查项集的相对频率：

    > itemFrequencyPlot(Groceries, support = 0.1, cex.names=0.8, topN=5)

    如何操作...

    食物事务的前五项频率条形图

  4. 使用apriori发现支持率超过 0.001 且置信度超过 0.5 的规则：

    > rules = apriori(Groceries, parameter = list(supp = 0.001, conf = 0.5, target= "rules"))
> summary(rules)
set of 5668 rules

rule length distribution (lhs + rhs):sizes
 2    3    4    5    6 
 11 1461 3211  939   46 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 2.00    3.00    4.00    3.92    4.00    6.00 

summary of quality measures:
 support           confidence          lift 
 Min.   :0.001017   Min.   :0.5000   Min.   : 1.957 
 1st Qu.:0.001118   1st Qu.:0.5455   1st Qu.: 2.464 
 Median :0.001322   Median :0.6000   Median : 2.899 
 Mean   :0.001668   Mean   :0.6250   Mean   : 3.262 
 3rd Qu.:0.001729   3rd Qu.:0.6842   3rd Qu.: 3.691 
 Max.   :0.022267   Max.   :1.0000   Max.   :18.996 

mining info:
 data ntransactions support confidence
 Groceries          9835   0.001        0.5

  5. 然后，我们可以检查前几条规则：

    > inspect(head(rules))
 lhs                    rhs              support confidence     lift
1 {honey}             => {whole milk} 0.001118454  0.7333333 2.870009
2 {tidbits}           => {rolls/buns} 0.001220132  0.5217391 2.836542
3 {cocoa drinks}      => {whole milk} 0.001321810  0.5909091 2.312611
4 {pudding powder}    => {whole milk} 0.001321810  0.5652174 2.212062
5 {cooking chocolate} => {whole milk} 0.001321810  0.5200000 2.035097
6 {cereals}           => {whole milk} 0.003660397  0.6428571 2.515917

  6. 您可以根据置信度对规则进行排序并检查前几条规则：

    > rules=sort(rules, by="confidence", decreasing=TRUE)
> inspect(head(rules))
 lhs                     rhs                    support confidence     lift
1 {rice, 
 sugar}              => {whole milk}       0.001220132          1 3.913649
2 {canned fish, 
 hygiene articles}   => {whole milk}       0.001118454          1 3.913649
3 {root vegetables, 
 butter, 
 rice}               => {whole milk}       0.001016777          1 3.913649
4 {root vegetables, 
 whipped/sour cream, 
 flour}              => {whole milk}       0.001728521          1 3.913649
5 {butter, 
 soft cheese, 
 domestic eggs}      => {whole milk}       0.001016777          1 3.913649
6 {citrus fruit, 
 root vegetables, 
 soft cheese}        => {other vegetables} 0.001016777          1 5.168156

工作原理...

关联挖掘的目的是从事务数据库中挖掘项目之间的关联。通常，关联挖掘的过程是通过查找支持度大于最小支持度的项集来进行的。接下来，过程使用频繁项集生成具有大于最小置信度的强规则（例如，milk => bread；购买牛奶的客户很可能也会购买面包）。根据定义，关联规则可以表示为 X=>Y 的形式，其中 X 和 Y 是不相交的项集。我们可以衡量两个术语之间的关联强度：支持度和置信度。支持度表示规则在数据集中适用的百分比，而置信度表示 X 和 Y 同时出现在同一交易中的概率：

  • 支持度 = 如何工作...

  • 置信度 = 如何工作...

在这里，如何工作... 指的是特定项集的频率；N 表示人口数量。

由于支持度和置信度是衡量规则强度的指标，你仍然可能会获得许多具有高支持度和置信度的冗余规则。因此，我们可以使用第三个度量，即提升度，来评估规则的质量（排名）。根据定义，提升度表明规则相对于 X 和 Y 的随机共现的强度，因此我们可以以下形式制定提升度：

升力 = 如何工作...

Apriori 是挖掘关联规则最著名的算法，它执行一种按层级的广度优先算法来计数候选项集。Apriori 的过程从按层级查找频繁项集（具有最小支持度的项集集合）开始。例如，过程从查找频繁的 1-项集开始。然后，过程继续使用频繁的 1-项集来查找频繁的 2-项集。过程迭代地从频繁的 k-项集中发现新的频繁 k+1-项集，直到没有找到频繁项集。

最后，过程利用频繁项集生成关联规则：

如何工作...

Apriori 算法示意图（其中支持度 = 2）

在这个菜谱中，我们使用 Apriori 算法在事务中查找关联规则。我们使用内置的 Groceries 数据集，它包含来自典型杂货店的现实世界一个月的销售额交易数据。然后，我们使用 summary 函数来获取 Groceries 数据集的摘要统计信息。摘要统计显示，数据集包含 9,835 笔交易，分为 169 个类别。此外，摘要还显示了数据集中的信息，例如最频繁的项、项集分布以及数据集中的示例扩展项信息。然后，我们可以使用 itemFrequencyPlot 来可视化支持率超过 0.1 的五个最频繁项。

接下来，我们应用 Apriori 算法来搜索支持度超过 0.001 且置信度超过 0.5 的规则。然后，我们使用 summary 函数来检查生成规则的详细信息。从输出摘要中，我们发现 Apriori 算法生成了 5,668 条支持度超过 0.001 且置信度超过 0.5 的规则。进一步，我们可以找到规则长度分布、质量度量摘要和挖掘信息。在质量度量摘要中，我们发现三个度量（支持度、置信度和提升度）的描述性统计。支持度是包含特定项集的交易比例。置信度是规则的正确百分比。提升度是响应目标关联规则除以平均响应。

要探索一些生成的规则，我们可以使用 inspect 函数来查看 5,668 条生成规则中的前六条规则。最后，我们可以根据置信度对规则进行排序，并列出置信度最高的规则。因此，我们发现与 whole milk 相关的 rich sugar 是置信度最高的规则，其支持度为 0.001220132，置信度为 1，提升度为 3.913649。

参考以下内容

对于对使用 Groceries 数据集的研究结果感兴趣的人，以及支持度、置信度和提升度测量是如何定义的，您可以参考以下论文：

  • Michael Hahsler，Kurt Hornik 和 Thomas Reutterer（2006）概率数据建模对挖掘关联规则的影响。在 M. Spiliopoulou，R. Kruse，C. Borgelt，A 中

  • Nuernberger 和 W. Gaul 编著，从数据和信息分析到知识工程，分类研究，数据分析与知识组织，第 598-605 页。Springer-Verlag

此外，除了使用 summary 和 inspect 函数来检查关联规则外，您还可以使用 interestMeasure 获取额外的兴趣度量：

> head(interestMeasure(rules, c("support", "chiSquare", "confidence", "conviction","cosine", "coverage", "leverage", "lift","oddsRatio"), Groceries))

剪枝冗余规则

在生成的规则中，我们有时会发现重复或冗余的规则（例如，一条规则是另一条规则的超级规则或子集）。在这个菜谱中，我们将向您展示如何剪枝（或删除）重复或冗余的规则。

准备工作

在这个菜谱中，您必须通过生成规则并存储在变量 rules 中来完成前面的菜谱。

如何做...

执行以下步骤以剪枝冗余规则：

  1. 首先，按照以下步骤查找冗余规则：

    > rules.sorted = sort(rules, by="lift")
> subset.matrix = is.subset(rules.sorted, rules.sorted)
> subset.matrix[lower.tri(subset.matrix, diag=T)] = NA
> redundant = colSums(subset.matrix, na.rm=T) >= 1

  2. 然后，您可以移除冗余规则：

    > rules.pruned = rules.sorted[!redundant]
> inspect(head(rules.pruned))
 lhs                        rhs                  support confidence     lift
1 {Instant food products, 
 soda}                  => {hamburger meat} 0.001220132  0.6315789 18.99565
2 {soda, 
 popcorn}               => {salty snack}    0.001220132  0.6315789 16.69779
3 {flour, 
 baking powder}         => {sugar}          0.001016777  0.5555556 16.40807
4 {ham, 
 processed cheese}      => {white bread}    0.001931876  0.6333333 15.04549
5 {whole milk, 
 Instant food products} => {hamburger meat} 0.001525165  0.5000000 15.03823
6 {other vegetables, 
 curd, 
 yogurt, 
 whipped/sour cream}    => {cream cheese }  0.001016777  0.5882353 14.83409

它是如何工作的...

关联挖掘的两个主要约束是在支持度和置信度之间进行选择。例如，如果你使用高的支持度阈值，你可能会移除稀有的项目规则，而不考虑这些规则是否具有高的置信度值。另一方面，如果你选择使用低的支持度阈值，关联挖掘可能会产生大量的冗余关联规则，这使得这些规则难以利用和分析。因此，我们需要剪枝冗余规则，以便从这些生成的规则中找到有意义的信息。

在这个菜谱中，我们演示了如何剪枝冗余规则。首先，我们搜索冗余规则。我们根据提升度对规则进行排序，然后使用is.subset函数找到排序规则的子集，这将生成一个itemMatrix对象。然后，我们可以将矩阵的下三角设置为 NA。最后，我们计算生成的矩阵的colSums，其中colSums >=1表示该特定规则是冗余的。

在我们找到冗余规则后，我们可以从排序规则中剪除这些规则。最后，我们可以使用inspect函数检查剪枝规则。

参见

  • 为了找到规则的子集或超集，你可以在关联规则上使用is.superset和is.subset函数。这两种方法可能会生成一个itemMatrix对象来显示哪个规则是其他规则的子集或超集。你可以参考help函数获取更多信息：

    > help(is.superset)
> help(is.subset)

关联规则的可视化

除了以文本形式列出规则外，你还可以可视化关联规则，这有助于找到项目集之间的关系。在下面的菜谱中，我们将介绍如何使用aruleViz包来可视化关联规则。

准备工作

在这个菜谱中，我们将继续使用Groceries数据集。你需要完成上一个菜谱，生成剪枝规则rules.pruned。

如何操作...

执行以下步骤以可视化关联规则：

  1. 首先，你需要安装并加载arulesViz包：

    > install.packages("arulesViz")
> library(arulesViz)

  2. 然后，你可以从剪枝规则中制作散点图：

    > plot(rules.pruned)

    如何操作...

    剪枝关联规则的散点图

  3. 此外，为了防止重叠，你可以在散点图中添加抖动：

    > plot(rules.pruned, shading="order", control=list(jitter=6))

    如何操作...

    关联规则剪枝散点图带有抖动

  4. 然后，我们使用 Apriori 算法在左侧使用soda生成新规则：

    > soda_rule=apriori(data=Groceries, parameter=list(supp=0.001,conf = 0.1, minlen=2), appearance = list(default="rhs",lhs="soda"))

  5. 接下来，你可以在图形图中绘制soda_rule：

    > plot(sort(soda_rule, by="lift"), method="graph", control=list(type="items"))

    如何操作...

    关联规则的图形图

  6. 此外，关联规则还可以在气球图中可视化：

    > plot(soda_rule, method="grouped")

    如何操作...

    关联规则的气球图

工作原理...

除了以文本形式呈现关联规则之外，还可以使用arulesViz来可视化关联规则。arulesViz是arules的扩展包，它提供了许多可视化技术来探索关联规则。要开始使用arulesViz，首先安装并加载arulesViz包。然后我们使用前一个菜谱中生成的剪枝规则制作散点图。根据步骤 2 中的图，我们发现规则在散点图中显示为点，x 轴表示支持度，y 轴表示置信度。颜色的深浅显示了规则的提升度；颜色越深，提升度越高。接下来，为了防止点重叠，我们可以在控制列表中包含 jitter 作为参数。添加了 jitter 的图在步骤 3 中的图中提供。

除了在散点图中绘制规则之外，arulesViz还允许您在图中和分组矩阵中绘制规则。我们不是在单个图中打印所有规则，而是选择使用soda在左侧生成新规则。然后我们使用提升排序规则，并在步骤 4 中的图中可视化规则。从图中可以看出，左侧的soda到右侧的全脂牛奶的规则具有最大的支持度，因为节点的大小最大。此外，该规则显示左侧的soda到右侧的瓶装水具有最大的提升度，因为圆圈的颜色最深。然后我们可以使用左侧带有soda的相同数据生成一个分组矩阵，这是一个在步骤 5 中的图中显示的气球图，左侧规则作为列标签，右侧作为行标签。类似于步骤 4 中的图中的图，步骤 5 中的图中的气球大小显示了规则的支持度，气球的颜色显示了规则的提升度。

另请参阅

  • 在这个菜谱中，我们介绍了三种可视化方法来绘制关联规则。然而，arulesViz也提供了绘制平行坐标图、双层图、马赛克图和其他相关图表的功能。对于那些对如何实现这些图表感兴趣的人，您可以参考：Hahsler, M. 和 Chelluboina, S. (2011). 可视化关联规则：R 扩展包 arulesViz 的介绍。R 项目模块。

  • 除了生成静态图之外，您可以通过以下步骤通过将交互设置为 TRUE 来生成交互式图：

    > plot(rules.pruned,interactive=TRUE)

    另请参阅

    交互式散点图

使用 Eclat 挖掘频繁项集

除了 Apriori 算法外，您还可以使用 Eclat 算法来生成频繁项集。由于 Apriori 算法执行广度优先搜索以扫描整个数据库，支持度计数相当耗时。如果数据库可以放入内存中，您可以使用 Eclat 算法，它执行深度优先搜索来计数支持度。因此，Eclat 算法比 Apriori 算法运行得更快。在这个食谱中，我们将介绍如何使用 Eclat 算法生成频繁项集。

准备工作

在这个食谱中，我们将继续使用数据集Groceries作为我们的输入数据源。

如何操作...

执行以下步骤以使用 Eclat 算法生成频繁项集：

  1. 与 Apriori 方法类似，我们可以使用eclat函数来生成频繁项集：

    > frequentsets=eclat(Groceries,parameter=list(support=0.05,maxlen=10))

  2. 然后，我们可以从生成的频繁项集中获取汇总信息：

    > summary(frequentsets)
set of 31 itemsets

most frequent items:
 whole milk other vegetables           yogurt 
 4                2                2 
 rolls/buns      frankfurter          (Other) 
 2                1               23 

element (itemset/transaction) length distribution:sizes
 1  2 
28  3 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 1.000   1.000   1.000   1.097   1.000   2.000 

summary of quality measures:
 support 
 Min.   :0.05236 
 1st Qu.:0.05831 
 Median :0.07565 
 Mean   :0.09212 
 3rd Qu.:0.10173 
 Max.   :0.25552 

includes transaction ID lists: FALSE 

mining info:
 data ntransactions support
 Groceries          9835    0.05

  3. 最后，我们可以检查前十个支持度最高的频繁项集：

    > inspect(sort(frequentsets,by="support")[1:10])
 items                 support
1  {whole milk}       0.25551601
2  {other vegetables} 0.19349263
3  {rolls/buns}       0.18393493
4  {soda}             0.17437722
5  {yogurt}           0.13950178
6  {bottled water}    0.11052364
7  {root vegetables}  0.10899847
8  {tropical fruit}   0.10493137
9  {shopping bags}    0.09852567
10 {sausage}          0.09395018

工作原理...

在这个食谱中，我们介绍另一个算法，Eclat，用于执行频繁项集生成。尽管 Apriori 是一种简单易懂的关联挖掘方法，但该算法的缺点是生成巨大的候选集，在支持度计数方面效率低下，因为它需要多次扫描数据库。相比之下，Eclat 使用等价类、深度优先搜索和集合交集，这大大提高了支持度计数的速度。

在 Apriori 算法中，算法使用水平数据布局来存储事务。另一方面，Eclat 使用垂直数据布局来存储每个项目的交易 ID 列表（tid）。然后，Eclat 通过交集两个 k-item 集的 tid 列表来确定任何 k+1-item 集的支持度。最后，Eclat 利用频繁项集来生成关联规则：

工作原理...

Eclat 算法的说明

与使用 Apriori 算法的食谱类似，我们可以使用eclat函数来生成具有给定支持度（本例中假设支持度为 2）和最大长度的频繁项集。

工作原理...

生成频繁项集

然后，我们可以使用summary函数来获取汇总统计信息，包括：最频繁的项目、项集长度分布、质量度量摘要和挖掘信息。最后，我们可以按支持度对频繁项集进行排序，并检查前十个支持度最高的频繁项集。

相关内容

  • 除了 Apriori 和 Eclat，另一个流行的关联挖掘算法是FP-Growth。与 Eclat 类似，它也采用深度优先搜索来计数支持度。然而，目前没有可以从 CRAN 下载的 R 包包含此算法。但是，如果您对在事务数据集中应用 FP-growth 算法感兴趣，您可以参考 Christian Borgelt 的页面www.borgelt.net/fpgrowth.html以获取更多信息。

创建带有时间信息的交易

除了在事务数据库中挖掘有趣的关联之外，我们还可以使用带有时间信息的交易挖掘有趣的序列模式。在下面的菜谱中，我们展示了如何创建带有时间信息的交易。

准备工作

在这个菜谱中，我们将生成带有时间信息的交易。我们可以使用生成的交易作为频繁序列模式挖掘的输入源。

如何操作...

执行以下步骤以创建带有时间信息的交易：

  1. 首先，你需要安装并加载包arulesSequences：

    > install.packages("arulesSequences")
> library(arulesSequences)

  2. 你可以先创建一个包含购买记录的列表：

    > tmp_data=list(c("A"),
+                c("A","B","C"),
+                c("A","C"),
+                c("D"),
+                c("C","F"),
+                c("A","D"),
+                c("C"),
+                c("B","C"),
+                c("A","E"),
+                c("E","F"),
+                c("A","B"),
+                c("D","F"),
+                c("C"),
+                c("B"),
+                c("E"),
+                c("G"),
+                c("A","F"),
+                c("C"),
+                c("B"),
+                c("C"))

  3. 然后，你可以将列表转换为交易并添加时间信息：

    >names(tmp_data) = paste("Tr",c(1:20), sep = "")
>trans =  as(tmp_data,"transactions")
>transactionInfo(trans)$sequenceID=c(1,1,1,1,1,2,2,2,2,3,3,3,3,3,4,4,4,4,4,4)
>transactionInfo(trans)$eventID=c(10,20,30,40,50,10,20,30,40,10,20,30,40,50,10,20,30,40,50,60)
> trans
transactions in sparse format with
 20 transactions (rows) and
 7 items (columns)

  4. 接下来，你可以使用inspect函数来检查交易：

    > inspect(head(trans))
 items transactionID sequenceID eventID
1 {A}             Tr1          1      10
2 {A, 
 B, 
 C}             Tr2          1      20
3 {A, 
 C}             Tr3          1      30
4 {D}             Tr4          1      40
5 {C, 
 F}             Tr5          1      50
6 {A, 
 D}             Tr6          2      10

  5. 然后，你可以获取带有时间信息的交易的摘要信息：

    > summary(trans)
transactions as itemMatrix in sparse format with
 20 rows (elements/itemsets/transactions) and
 7 columns (items) and a density of 0.2214286 

most frequent items:
 C       A       B       F       D (Other) 
 8       7       5       4       3       4 

element (itemset/transaction) length distribution:
sizes
 1  2  3 
10  9  1 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 1.00    1.00    1.50    1.55    2.00    3.00 

includes extended item information - examples:
 labels
1      A
2      B
3      C

includes extended transaction information - examples:
 transactionID sequenceID eventID
1           Tr1          1      10
2           Tr2          1      20
3           Tr3          1      30

  6. 你也可以以篮子格式读取事务数据：

    > zaki=read_baskets(con = system.file("misc", "zaki.txt", package = "arulesSequences"), info = c("sequenceID","eventID","SIZE"))
> as(zaki, "data.frame")
 transactionID.sequenceID transactionID.eventID transactionID.SIZE     items
1                         1                    10                  2     {C,D}
2                         1                    15                  3   {A,B,C}
3                         1                    20                  3   {A,B,F}
4                         1                    25                  4 {A,C,D,F}
5                         2                    15                  3   {A,B,F}
6                         2                    20                  1       {E}
7                         3                    10                  3   {A,B,F}
8                         4                    10                  3   {D,G,H}
9                         4                    20                  2     {B,F}
10                        4                    25                  3   {A,G,H}

工作原理...

在挖掘频繁序列模式之前，你必须创建带有时间信息的事务。在这个菜谱中，我们介绍了两种获取带有时间信息的事务的方法。在第一种方法中，我们创建一个事务列表，并为每个事务分配一个事务 ID。我们使用as函数将列表数据转换为事务数据集。然后我们添加eventID和sequenceID作为时间信息；sequenceID是事件所属的序列，eventID表示事件发生的时间。在生成带有时间信息的交易后，可以使用这个数据集进行频繁序列模式挖掘。

除了创建自己的带有时间信息的交易之外，如果你已经有一个存储在文本文件中的数据，你可以使用arulesSequences中的read_basket函数将事务数据读取到篮子格式。我们还可以读取事务数据集以进行进一步的频繁序列模式挖掘。

相关内容

  • arulesSequences函数提供了两个额外的数据结构，sequences和timedsequences，用于表示纯序列数据和带时间信息的序列数据。对于对这些两个集合感兴趣的人，请使用帮助函数查看以下文档：

    > help("sequences-class")
> help("timedsequences-class")

使用 cSPADE 挖掘频繁序列模式

与仅发现项集之间关系的关联挖掘不同，我们可能对探索在一系列项集按顺序发生的事务中共享的模式感兴趣。

最著名的频繁序列模式挖掘算法之一是使用等价类的序列模式发现（SPADE）算法，该算法利用垂直数据库的特性在 ID 列表上执行高效格搜索，并允许我们对挖掘的序列施加约束。在本食谱中，我们将演示如何使用 cSPADE 来挖掘频繁序列模式。

准备工作

在本食谱中，您必须通过生成带有时间信息的交易并存储在变量trans中来完成前面的食谱。

如何操作...

执行以下步骤以挖掘频繁序列模式：

  1. 首先，您可以使用cspade函数生成频繁序列模式：

    > s_result=cspade(trans,parameter = list(support = 0.75),control = list(verbose = TRUE))

  2. 然后，您可以检查频繁序列模式的摘要：

    > summary(s_result)
set of 14 sequences with

most frequent items:
 C       A       B       D       E (Other) 
 8       5       5       2       1       1 

most frequent elements:
 {C}     {A}     {B}     {D}     {E} (Other) 
 8       5       5       2       1       1 

element (sequence) size distribution:
sizes
1 2 3 
6 6 2 

sequence length distribution:
lengths
1 2 3 
6 6 2 

summary of quality measures:
 support 
 Min.   :0.7500 
 1st Qu.:0.7500 
 Median :0.7500 
 Mean   :0.8393 
 3rd Qu.:1.0000 
 Max.   :1.0000 

includes transaction ID lists: FALSE 

mining info:
 data ntransactions nsequences support
 trans            20          4    0.75

  3. 将生成的序列格式数据转换回数据框：

    > as(s_result, "data.frame")
 sequence support
1          <{A}>    1.00
2          <{B}>    1.00
3          <{C}>    1.00
4          <{D}>    0.75
5          <{E}>    0.75
6          <{F}>    0.75
7      <{A},{C}>    1.00
8      <{B},{C}>    0.75
9      <{C},{C}>    0.75
10     <{D},{C}>    0.75
11 <{A},{C},{C}>    0.75
12     <{A},{B}>    1.00
13     <{C},{B}>    0.75
14 <{A},{C},{B}>    0.75

它是如何工作的...

序列模式挖掘的目标是在事务中发现序列关系或模式。您可以使用模式挖掘结果来预测未来事件，或向用户推荐项目。

序列模式挖掘的一种流行方法是 SPADE。SPADE 使用垂直数据布局来存储 ID 列表。在这些列表中，数据库中的每个输入序列被称为 SID，给定输入序列中的每个事件被称为 EID。SPADE 的过程是通过按层次生成 Apriori 候选来进行的。具体来说，SPADE 通过从 ID 列表的交集生成(n-1)-序列来生成后续的 n-序列。如果序列的数量大于最小支持度（minsup），我们可以认为序列足够频繁。算法停止，直到无法找到更多频繁序列为止：

它是如何工作的...

SPADE 算法的说明

在本食谱中，我们说明了如何使用频繁序列模式挖掘算法 cSPADE 来挖掘频繁序列模式。首先，由于我们已经将带有时间信息的交易加载到变量trans中，我们可以使用cspade函数，支持度超过 0.75，以sequences格式生成频繁序列模式。然后我们可以获取总结信息，例如最频繁的项目、序列大小分布、质量度量摘要和挖掘信息。最后，我们可以将生成的sequence信息转换回数据框格式，这样我们就可以检查支持度超过 0.75 的频繁序列模式的序列和支持度。

参考内容

  • 如果您对 SPADE 算法的概念和设计感兴趣，您可以参考原始发表的论文：M. J. Zaki. (2001). SPADE: An Efficient Algorithm for Mining Frequent Sequences. Machine Learning Journal, 42, 31–60。

第十一章. 降维

在本章中，我们将涵盖以下主题：

  • 使用 FSelector 进行特征选择

  • 使用 PCA 进行降维

  • 使用 scree 测试确定主成分的数量

  • 使用 Kaiser 方法确定主成分的数量

  • 使用双图可视化多元数据

  • 使用 MDS 进行降维

  • 使用奇异值分解（SVD）降维

  • 使用 SVD 压缩图像

  • 使用 ISOMAP 进行非线性降维

  • 使用局部线性嵌入进行非线性降维

简介

大多数数据集都包含高度冗余的特征（如属性或变量）。为了去除无关和冗余的数据，以降低计算成本并避免过拟合，你可以将特征减少到较小的子集，而不会造成信息损失。减少特征的数学过程称为降维。

特征的减少可以提高数据处理效率。因此，降维在模式识别、文本检索和机器学习等领域得到广泛应用。降维可以分为两个部分：特征提取和特征选择。特征提取是一种使用低维空间来表示高维空间中数据的技巧。特征选择用于找到原始变量的子集。

特征选择的目的是选择一组相关特征来构建模型。特征选择的技术可以分为特征排名和特征选择。特征排名根据一定标准对特征进行排名，然后选择超过定义阈值的特征。另一方面，特征选择从特征子集空间中搜索最优子集。

在特征提取中，问题可以分为线性或非线性。线性方法搜索一个最佳解释数据分布变化的仿射空间。相比之下，非线性方法对于分布在高度非线性曲面上数据的处理是一个更好的选择。在此，我们列出一些常见的线性和非线性方法。

这里有一些常见的线性方法：

  • PCA：主成分分析将数据映射到低维空间，使得低维表示中的数据方差最大化。

  • MDS：多维尺度分析是一种方法，它允许你可视化对象之间的接近程度（模式邻近性），并且可以使用低维空间表示你的数据。如果 MDS 中使用的距离测量等于数据的协方差，则 PCA 可以被视为 MDS 的最简单形式。

  • SVD：奇异值分解从线性代数的角度移除了线性相关的冗余特征。PCA 也可以被视为 SVD 的特例。

这里有一些常见的非线性方法：

  • ISOMAP：ISOMAP 可以看作是多维尺度分析（MDS）的扩展，它使用测地距离的距离度量。在此方法中，测地距离是通过绘制最短路径距离来计算的。

  • LLE：局部线性嵌入（Locally linear embedding）执行局部 PCA 和全局特征分解。LLE 是一种局部方法，它涉及为每个类特征类别选择特征。相比之下，ISOMAP 是一种全局方法，它涉及为所有特征选择特征。

在本章中，我们将首先讨论如何进行特征排名和选择。接下来，我们将聚焦于特征提取的主题，并涵盖使用线性和非线性方法进行降维的食谱。对于线性方法，我们将介绍如何执行主成分分析（PCA），确定主成分的数量及其可视化。然后，我们将继续介绍多维尺度分析（MDS）和奇异值分解（SVD）。此外，我们将介绍 SVD 在图像压缩中的应用。对于非线性方法，我们将介绍如何使用 ISOMAP 和 LLE 进行降维。

使用 FSelector 进行特征选择

FSelector包提供了两种从原始特征集中选择最具影响力的特征的方法。首先，根据某些标准对特征进行排名，并选择那些高于定义阈值的特征。其次，从一个特征子集的空间中搜索最佳特征子集。在本食谱中，我们将介绍如何使用FSelector包进行特征选择。

准备工作

在这个食谱中，我们将继续使用电信churn数据集作为输入数据源来训练支持向量机。对于那些尚未准备数据集的人，请参阅第五章, 分类（I） – 树、懒惰和概率，以获取详细信息。

如何操作...

执行以下步骤以在churn数据集上执行特征选择：

  1. 首先，安装并加载包，FSelector：

    > install.packages("FSelector")
> library(FSelector)

  2. 然后，我们可以使用random.forest.importance计算每个属性的权重，我们将重要性类型设置为 1：

    > weights = random.forest.importance(churn~., trainset, importance.type = 1)
> print(weights)
 attr_importance
international_plan                 96.3255882
voice_mail_plan                    24.8921239
number_vmail_messages              31.5420332
total_day_minutes                  51.9365357
total_day_calls                    -0.1766420
total_day_charge                   53.7930096
total_eve_minutes                  33.2006078
total_eve_calls                    -2.2270323
total_eve_charge                   32.4317375
total_night_minutes                22.0888120
total_night_calls                   0.3407087
total_night_charge                 21.6368855
total_intl_minutes                 32.4984413
total_intl_calls                   51.1154046
total_intl_charge                  32.4855194
number_customer_service_calls     114.2566676

  3. 接下来，我们可以使用cutoff函数获取前五个权重的属性：

    > subset = cutoff.k(weights, 5)
> f = as.simple.formula(subset, "Class")
> print(f)
Class ~ number_customer_service_calls + international_plan + 
 total_day_charge + total_day_minutes + total_intl_calls
<environment: 0x00000000269a28e8>

  4. 接下来，我们可以创建一个评估器来选择特征子集：

    > evaluator = function(subset) {
+   k = 5 
+   set.seed(2)
+   ind = sample(5, nrow(trainset), replace = TRUE)
+   results = sapply(1:k, function(i) {
+     train = trainset[ind ==i,]
+     test  = trainset[ind !=i,]
+     tree  = rpart(as.simple.formula(subset, "churn"), trainset)
+     error.rate = sum(test$churn != predict(tree, test, type="class")) / nrow(test)
+     return(1 - error.rate)
+   })
+   return(mean(results))
+ }

  5. 最后，我们可以使用爬山搜索找到最佳特征子集：

    > attr.subset = hill.climbing.search(names(trainset)[!names(trainset) %in% "churn"], evaluator)
> f = as.simple.formula(attr.subset, "churn")
> print(f)
churn ~ international_plan + voice_mail_plan + number_vmail_messages + 
 total_day_minutes + total_day_calls + total_eve_minutes + 
 total_eve_charge + total_intl_minutes + total_intl_calls + 
 total_intl_charge + number_customer_service_calls
<environment: 0x000000002224d3d0>

工作原理...

在这个示例中，我们展示了如何使用FSelector包来选择最有影响力的特征。我们首先演示了如何使用特征排名方法。在特征排名方法中，算法首先使用一个权重函数为每个特征生成权重。在这里，我们使用随机森林算法，以平均准确度下降（其中importance.type = 1）作为重要性度量来获取每个属性的权重。除了随机森林算法之外，您还可以从FSelector包中选择其他特征排名算法（例如，chi.squared，information.gain）。然后，该过程按权重对属性进行排序。最后，我们可以使用cutoff函数从排序的特征列表中获得前五个特征。在这种情况下，number_customer_service_calls，international_plan，total_day_charge，total_day_minutes和total_intl_calls是五个最重要的特征。

接下来，我们将说明如何搜索最佳特征子集。首先，我们需要创建一个五折交叉验证函数来评估特征子集的重要性。然后，我们使用爬山搜索算法从原始特征集中找到最佳特征子集。除了爬山法之外，您还可以从FSelector包中选择其他特征选择算法（例如，forward.search）。最后，我们可以发现international_plan + voice_mail_plan + number_vmail_messages + total_day_minutes + total_day_calls + total_eve_minutes + total_eve_charge + total_intl_minutes + total_intl_calls + total_intl_charge + number_customer_service_calls是最佳特征子集。

参考内容

  • 您也可以使用caret包来进行特征选择。正如我们在模型评估章节中讨论了相关食谱，您可以参考第七章，模型评估，以获取更详细的信息。

  • 对于特征排名和最佳特征选择，您可以探索FSelector包以获取更多相关函数：

    > help(package="FSelector")

使用 PCA 进行降维

主成分分析（PCA）是处理降维问题中最广泛使用的线性方法。当数据包含许多特征，并且这些特征之间存在冗余（相关性）时，它非常有用。为了去除冗余特征，PCA 通过将特征减少到更少的、解释原始特征大部分方差的主成分，将高维数据映射到低维。在这个示例中，我们将介绍如何使用 PCA 方法进行降维。

准备工作

在这个示例中，我们将使用swiss数据集作为我们的目标来执行 PCA。swiss数据集包括来自瑞士 47 个法语省份大约 1888 年的标准化生育指标和社会经济指标。

如何操作...

执行以下步骤以对swiss数据集进行主成分分析：

  1. 首先，加载swiss数据集：

    > data(swiss)

  2. 排除swiss数据的第一列：

    > swiss = swiss[,-1]

  3. 然后，你可以对swiss数据进行主成分分析：

    > swiss.pca = prcomp(swiss,
+ center = TRUE,
+ scale  = TRUE)
> swiss.pca
Standard deviations:
[1] 1.6228065 1.0354873 0.9033447 0.5592765 0.4067472

Rotation:
 PC1         PC2          PC3        PC4         PC5
Agriculture      0.52396452 -0.25834215  0.003003672 -0.8090741  0.06411415
Examination  -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
Education       -0.49150243  0.19028476  0.539337412 -0.3321615  0.56656945
Catholic            0.38530580  0.36956307  0.725888143 0.1007965 -0.42176895
Infant.Mortality 0.09167606 0.87197641 -0.424976789 -0.2154928 0.06488642

  4. 从 PCA 结果中获取摘要：

    > summary(swiss.pca)
Importance of components:
 PC1    PC2    PC3     PC4     PC5
Standard deviation     1.6228 1.0355 0.9033 0.55928 0.40675
Proportion of Variance 0.5267 0.2145 0.1632 0.06256 0.03309
Cumulative Proportion  0.5267 0.7411 0.9043 0.96691 1.00000

  5. 最后，你可以使用predict函数输出数据第一行的主成分值：

    > predict(swiss.pca, newdata=head(swiss, 1))
 PC1       PC2        PC3      PC4       PC5
Courtelary -0.9390479 0.8047122 -0.8118681 1.000307 0.4618643

它是如何工作的...

由于特征选择方法可能会移除一些相关但具有信息量的特征，你必须考虑使用特征提取方法将这些相关特征组合成一个单一的特征。PCA 是特征提取方法之一，它通过正交变换将可能相关的变量转换为主成分。此外，你可以使用这些主成分来识别方差的方向。

PCA 的过程通过以下步骤进行：首先，找到均值向量，它是如何工作的...，其中它是如何工作的...表示数据点，n表示点的数量。其次，通过方程计算协方差矩阵，它是如何工作的...。第三，计算特征向量，它是如何工作的...，以及相应的特征值。在第四步，我们按顺序选择前k个特征向量。在第五步，我们构建一个d x k维的特征向量矩阵，U。在这里，d是原始维度的数量，k是特征向量的数量。最后，我们可以在方程中变换数据样本到新的子空间，它是如何工作的...。

在以下图中，说明了我们可以使用两个主成分，它是如何工作的...和它是如何工作的...，将数据点从二维空间转换到新的二维子空间：

它是如何工作的...

PCA 的一个示例说明

在这个菜谱中，我们使用stats包中的prcomp函数对swiss数据集进行 PCA。首先，我们移除标准化的生育率指标，并将剩余的预测变量作为输入传递给函数prcomp。此外，我们将swiss作为输入数据集；变量应通过指定center=TRUE移至零中心；通过选项scale=TRUE将变量缩放到单位方差，并将输出存储在变量swiss.pca中。

然后，当我们打印出存储在 swiss.pca 中的值时，我们可以找到主成分的标准差和旋转。标准差表示协方差/相关矩阵的特征值的平方根。另一方面，主成分的旋转显示了输入特征的线性组合的系数。例如，PC1 等于 Agriculture * 0.524 + Examination * -0.572 + Education * -0.492 + Catholic 0.385 + Infant.Mortality * 0.092*。在这里，我们可以发现属性 Agriculture 对 PC1 的贡献最大，因为它具有最高的系数。

此外，我们可以使用 summary 函数获取成分的重要性。第一行显示每个主成分的标准差，第二行显示每个成分解释的方差比例，第三行显示解释的方差累积比例。最后，你可以使用 predict 函数从输入特征中获取主成分。在这里，我们输入数据集的第一行，并检索五个主成分。

还有更多...

另一个主成分分析函数是 princomp。在这个函数中，计算是通过在相关或协方差矩阵上使用特征值而不是 prcomp 函数中使用的单一值分解来进行的。在一般实践中，使用 prcomp 是首选的；然而，我们在这里介绍了如何使用 princomp：

  1. 首先，使用 princomp 执行主成分分析（PCA）：

    > swiss.princomp = princomp(swiss,
+ center = TRUE,
+ scale  = TRUE)
> swiss.princomp
Call:
princomp(x = swiss, center = TRUE, scale = TRUE)

Standard deviations:
 Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
42.896335 21.201887  7.587978  3.687888  2.721105 

 5 variables and 47 observations.

  2. 然后，你可以获取摘要信息：

    > summary(swiss.princomp)
Importance of components:
 Comp.1     Comp.2     Comp.3      Comp.4      Comp.5
Standard deviation     42.8963346 21.2018868 7.58797830 3.687888330 2.721104713
Proportion of Variance  0.7770024  0.1898152 0.02431275 0.005742983 0.003126601
Cumulative Proportion   0.7770024  0.9668177 0.99113042 0.996873399 1.000000000

  3. 你可以使用 predict 函数从输入特征中获取主成分：

    > predict(swiss.princomp, swiss[1,])
 Comp.1    Comp.2   Comp.3   Comp.4   Comp.5
Courtelary -38.95923 -20.40504 12.45808 4.713234 -1.46634

除了 stats 包中的 prcomp 和 princomp 函数外，你还可以使用 psych 包中的 principal 函数：

  1. 首先，安装并加载 psych 包：

    > install.packages("psych")
> install.packages("GPArotation")
> library(psych)

  2. 然后，你可以使用 principal 函数检索主成分：

    > swiss.principal = principal(swiss, nfactors=5, rotate="none")
> swiss.principal
Principal Components Analysis
Call: principal(r = swiss, nfactors = 5, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
 PC1   PC2   PC3   PC4   PC5 h2       u2
Agriculture      -0.85 -0.27  0.00  0.45 -0.03  1 -6.7e-16
Examination       0.93 -0.01 -0.04  0.24  0.29  1  4.4e-16
Education         0.80  0.20  0.49  0.19 -0.23  1  2.2e-16
Catholic         -0.63  0.38  0.66 -0.06  0.17  1 -2.2e-16
Infant.Mortality -0.15  0.90 -0.38  0.12 -0.03  1 -8.9e-16

 PC1  PC2  PC3  PC4  PC5
SS loadings           2.63 1.07 0.82 0.31 0.17
Proportion Var        0.53 0.21 0.16 0.06 0.03
Cumulative Var        0.53 0.74 0.90 0.97 1.00
Proportion Explained  0.53 0.21 0.16 0.06 0.03
Cumulative Proportion 0.53 0.74 0.90 0.97 1.00

Test of the hypothesis that 5 components are sufficient.

The degrees of freedom for the null model are 10 and the objective function was 2.13
The degrees of freedom for the model are -5  and the objective function was  0 
The total number of observations was  47  with MLE Chi Square =  0  with prob <  NA 

Fit based upon off diagonal values = 1

使用斯皮尔曼测试确定主成分的数量

由于我们只需要保留解释原始特征大部分方差的主成分，我们可以使用凯撒方法、斯皮尔曼测试或解释的方差百分比作为选择标准。斯皮尔曼测试的主要目的是将成分分析结果以斯皮尔曼图的形式绘制出来，并找到斜率（肘部）发生明显变化的位置。在这个菜谱中，我们将演示如何使用斯皮尔曼图确定主成分的数量。

准备工作

确保你已经完成了前面的菜谱，通过生成主成分对象并将其保存在变量 swiss.pca 中。

如何操作...

执行以下步骤以使用斯皮尔曼图确定主成分的数量：

  1. 首先，你可以使用 screeplot 生成条形图：

    > screeplot(swiss.pca, type="barplot")

    如何操作...

    以条形图形式展示的斯皮尔曼图

  2. 你也可以使用 screeplot 生成线形图：

    > screeplot(swiss.pca, type="line")

    如何操作...

    以线形图形式展示的斯皮尔曼图

工作原理...

在本菜谱中，我们展示了如何使用碎石图来确定主成分的数量。在碎石图中，有两种类型的图形，即柱状图和线形图。正如生成的碎石图所揭示的，斜率（所谓的拐点或膝部）的明显变化发生在第 2 个成分处。因此，我们应该保留第 1 个成分，在第 2 个成分之前，该成分处于陡峭的曲线，这是平坦线趋势开始的地方。然而，由于这种方法可能存在歧义，您可以使用其他方法（如凯撒方法）来确定成分的数量。

更多内容...

默认情况下，如果您在生成的主成分对象上使用plot函数，您也可以检索碎石图。有关screeplot的更多详细信息，请参阅以下文档：

> help(screeplot)

您还可以使用nfactors进行并行分析和 Cattell 碎石测试的非图形解决方案：

> install.packages("nFactors")
> library(nFactors)
> ev = eigen(cor(swiss))
> ap = parallel(subject=nrow(swiss),var=ncol(swiss),rep=100,cent=.05)
> nS = nScree(x=ev$values, aparallel=ap$eigen$qevpea)
> plotnScree(nS)

更多内容...

碎石测试的非图形解决方案

使用凯撒方法确定主成分的数量

除了碎石测试外，您还可以使用凯撒方法来确定主成分的数量。在此方法中，选择标准是保留大于1的特征值。在本菜谱中，我们将展示如何使用凯撒方法确定主成分的数量。

准备工作

确保您已经完成了前面的菜谱，通过生成主成分对象并将其保存在变量swiss.pca中。

如何操作...

按照以下步骤使用凯撒方法确定主成分的数量：

  1. 首先，您可以从swiss.pca中获取标准差：

    > swiss.pca$sdev 
[1] 1.6228065 1.0354873 0.9033447 0.5592765 0.4067472

  2. 接着，您可以从swiss.pca中获取方差：

    > swiss.pca$sdev ^ 2
[1] 2.6335008 1.0722340 0.8160316 0.3127902 0.1654433

  3. 选择方差大于 1 的成分：

    > which(swiss.pca$sdev ^ 2> 1)
[1] 1 2

  4. 您还可以使用碎石图来选择方差大于 1 的成分：

    > screeplot(swiss.pca, type="line")
> abline(h=1, col="red", lty= 3)

    如何操作...

    选择方差大于 1 的成分

工作原理...

您还可以使用凯撒方法来确定成分的数量。由于计算出的主成分对象包含每个成分的标准差，我们可以将方差计算为标准差，即方差的平方根。从计算出的方差中，我们发现第 1 个和第 2 个成分的方差都大于 1。因此，我们可以确定主成分的数量为 2（第 1 个和第 2 个成分）。此外，我们可以在碎石图上画一条红线（如图所示）来表示在这种情况下我们需要保留第 1 个和第 2 个成分。

参考资料也

为了确定要保留哪些主成分，请参考以下内容：

  • Ledesma, R. D. 和 Valero-Mora, P. (2007). 在 EFA 中确定要保留的因素数量：一个用于执行并行分析的易于使用的计算机程序. 实用评估、研究和评估, 12(2), 1-11.

使用双图可视化多元数据

为了找出数据变量如何映射到主成分上，您可以使用biplot，它将数据以及原始特征在第一两个成分上的投影绘制出来。在这个菜谱中，我们将演示如何使用biplot在同一张图上绘制变量和数据。

准备工作

确保您已经完成了之前的菜谱，通过生成主成分对象并将其保存在变量swiss.pca中。

如何操作...

执行以下步骤以创建 biplot：

  1. 您可以使用成分 1 和 2 创建散点图：

    >  plot(swiss.pca$x[,1], swiss.pca$x[,2], xlim=c(-4,4))
> text(swiss.pca$x[,1], swiss.pca$x[,2], rownames(swiss.pca$x), cex=0.7, pos=4, col="red")

    如何操作...

    PCA 结果的前两个成分的散点图

  2. 如果您想在图上添加特征，您可以使用生成的主成分对象创建 biplot：

    > biplot(swiss.pca)

    如何操作...

    使用 PCA 结果的 biplot

工作原理...

在这个菜谱中，我们演示了如何使用biplot绘制数据和原始特征在第一两个成分上的投影。在第一步中，我们演示了实际上我们可以使用前两个成分来创建散点图。此外，如果您想在同一张图上添加变量，您可以使用biplot。在biplot中，您可以看到在农业变量中指标较高的省份，在教育变量中指标较低的省份，以及在 PC1 中得分较高的考试变量。另一方面，婴儿死亡率指标较高且农业指标较低的省份在 PC2 中得分较高。

还有更多...

除了stats包中的biplot，您还可以使用ggbiplot。但是，您可能无法从 CRAN 找到这个包；您必须首先安装devtools，然后从 GitHub 安装ggbiplot：

> install.packages("devtools")

> library(ggbiplot)
> g = ggbiplot(swiss.pca, obs.scale = 1, var.scale = 1, 
+ ellipse = TRUE, 
+ circle = TRUE)
> print(g)

还有更多...

使用 PCA 结果的 ggbiplot

使用 MDS 进行降维

多维尺度分析（MDS）是一种创建多个对象相似性或差异性（距离）的视觉表示的技术。多前缀表示可以在一维、二维或更多维度中创建一个表示图。然而，我们最常使用 MDS 来展示一维或二维数据点之间的距离。

在 MDS 中，您可以使用度量或非度量解决方案。两种解决方案之间的主要区别在于，度量解决方案试图重现原始度量，而非度量解决方案假设距离的排名是已知的。在这个菜谱中，我们将说明如何在swiss数据集上执行 MDS。

准备工作

在这个菜谱中，我们将继续使用swiss数据集作为我们的输入数据源。

如何操作...

执行以下步骤以使用度量方法执行多维尺度分析：

  1. 首先，您可以使用最多两个维度执行度量 MDS：

    > swiss.dist =dist(swiss)
> swiss.mds = cmdscale(swiss.dist, k=2)

  2. 然后，您可以将swiss数据绘制在二维散点图中：

    > plot(swiss.mds[,1], swiss.mds[,2], type = "n", main = "cmdscale (stats)")
> text(swiss.mds[,1], swiss.mds[,2], rownames(swiss), cex = 0.9, xpd = TRUE)

    如何操作...

    从 cmdscale 对象得到的二维散点图

  3. 此外，您可以使用isoMDS执行非度量 MDS：

    > library(MASS)
> swiss.nmmds = isoMDS(swiss.dist, k=2)
initial  value 2.979731 
iter   5 value 2.431486
iter  10 value 2.343353
final  value 2.338839 
converged

  4. 你也可以在二维散点图上绘制数据点：

    > plot(swiss.nmmds$points, type = "n", main = "isoMDS (MASS)")
> text(swiss.nmmds$points, rownames(swiss), cex = 0.9, xpd = TRUE)

    如何操作...

    从 isoMDS 对象生成的二维散点图

  5. 你还可以在二维散点图上绘制数据点：

    > swiss.sh = Shepard(swiss.dist, swiss.mds)
> plot(swiss.sh, pch = ".")
> lines(swiss.sh$x, swiss.sh$yf, type = "S")

    如何操作...

    从 isoMDS 对象生成的谢泼德图

如何工作...

MDS 通过提供一组对象之间相似性的视觉表示来揭示数据的结构。更详细地说，MDS 将一个对象放置在 n 维空间中，其中点对之间的距离对应于对象对之间的相似性。通常，维空间是二维欧几里得空间，但它可能是非欧几里得空间，并且具有超过两个维度。根据输入矩阵的意义，MDS 可以主要分为两种类型：度量 MDS，其中输入矩阵基于度量，非度量 MDS，其中输入矩阵基于非度量。

度量 MDS 也称为主坐标分析，它首先将距离转换为相似性。在最简单的情况下，该过程通过在相似性上执行主成分分析，将原始数据点线性投影到子空间。另一方面，该过程也可以通过最小化应力值如何工作...对相似性进行非线性投影，其中如何工作...是两点之间的距离测量，如何工作...和如何工作...，而如何工作...是两个投影点的相似性度量，如何工作...和如何工作...。因此，我们可以在欧几里得空间中表示对象之间的关系。

与使用基于度量的输入矩阵的度量 MDS 相比，当数据在序数级别测量时，使用基于非度量的 MDS。由于只有向量之间距离的秩序是有意义的，非度量 MDS 在原始距离上应用单调递增函数 f，并将距离投影到新的值，以保留秩序。归一化方程可以表示为如何工作...。

在这个菜谱中，我们展示了如何在swiss数据集上执行度量和非度量 MDS。要执行度量 MDS，我们首先需要从swiss数据中获得距离度量。在这个步骤中，你可以将距离度量替换为任何度量，只要它产生数据点的相似性/不相似性度量。你可以使用cmdscale来执行度量多维度缩放。在这里，我们指定k = 2，因此生成的最大维度等于2。你还可以在二维散点图上直观地表示数据点的距离。

接下来，你可以使用isoMDS执行非度量 MDS。在非度量 MDS 中，我们不匹配距离，而只是按顺序排列它们。我们还将swiss设置为输入数据集，最大维度为两个。与度量 MDS 示例类似，我们可以在二维散点图上绘制数据点之间的距离。然后，我们使用谢泼德图，它显示了投影距离与距离矩阵中的距离匹配得有多好。根据步骤 4 中的图，投影距离在距离矩阵中匹配得很好。

还有更多...

另一种可视化方法是将以图形的形式表示 MDS 对象。以下是一个示例代码：

> library(igraph)
> swiss.sample = swiss[1:10,]

> g = graph.full(nrow(swiss.sample))
> V(g)$label = rownames(swiss.sample)
> layout = layout.mds(g, dist = as.matrix(dist(swiss.sample)))
> plot(g, layout = layout, vertex.size = 3)

还有更多...

MDS 对象的图形表示

你还可以比较 MDS 和 PCA 生成的结果之间的差异。你可以通过在同一散点图上绘制投影维度来比较它们的差异。如果你在 MDS 中使用欧几里得距离，则投影维度与从 PCA 投影的维度完全相同：

> swiss.dist = dist(swiss)
> swiss.mds = cmdscale(swiss.dist, k=2)
> plot(swiss.mds[,1], swiss.mds[,2], type="n")
> text(swiss.mds[,1], swiss.mds[,2], rownames(swiss), cex = 0.9, xpd = TRUE)
> swiss.pca = prcomp(swiss)
> text(-swiss.pca$x[,1],-swiss.pca$x[,2], rownames(swiss), 
+      ,col="blue", adj = c(0.2,-0.5),cex = 0.9, xpd = TRUE)

还有更多...

MDS 与 PCA 的比较

使用 SVD 降维

奇异值分解（SVD）是一种矩阵分解（分解）类型，可以将矩阵分解为两个正交矩阵和对角矩阵。你可以使用这三个矩阵将原始矩阵乘回。从线性代数的角度来看，SVD 可以减少线性相关的冗余数据。因此，它可以应用于特征选择、图像处理、聚类以及许多其他领域。在这个菜谱中，我们将说明如何使用 SVD 进行降维。

准备工作

在这个菜谱中，我们将继续使用数据集swiss作为我们的输入数据源。

如何操作...

执行以下步骤以使用 SVD 进行降维：

  1. 首先，你可以在swiss数据集上执行svd：

    > swiss.svd = svd(swiss)

  2. 然后，你可以根据 SVD 列绘制解释的方差百分比和累积解释的方差百分比：

    > plot(swiss.svd$d²/sum(swiss.svd$d²), type="l", xlab=" Singular vector", ylab = "Variance explained")

    如何操作...

    解释的方差百分比

    > plot(cumsum(swiss.svd$d²/sum(swiss.svd$d²)), type="l", xlab="Singular vector", ylab = "Cumulative percent of variance explained")

    如何操作...

    累积方差百分比

  3. 接下来，你可以使用仅一个奇异向量来重建数据：

    > swiss.recon = swiss.svd$u[,1] %*% diag(swiss.svd$d[1], length(1), length(1)) %*% t(swiss.svd$v[,1])

  4. 最后，你可以在图像中比较原始数据集与构建的数据集：

    > par(mfrow=c(1,2))
> image(as.matrix(swiss), main="swiss data Image")
> image(swiss.recon,  main="Reconstructed Image")

    如何操作...

    原始数据集与重建数据集之间的比较

如何工作...

SVD 是实数或复数矩阵的分解。具体来说，m x n 矩阵 A 的 SVD 是将 A 分解为三个矩阵的乘积如何工作...。在这里，U 是一个 m x m 的正交矩阵，D 包含奇异值，是一个 m x n 的对角矩阵，V^T 是一个 n x n 的正交矩阵。

在这个食谱中，我们展示了如何使用 SVD 进行降维。首先，你可以在swiss数据集上应用svd函数以获得分解矩阵。然后你可以生成两个图表：一个显示了根据奇异向量解释的方差，另一个显示了根据奇异向量解释的累积方差。

前面的图表明第一个奇异向量可以解释 80%的方差。我们现在想比较原始数据集和仅使用一个奇异向量重建的数据集之间的差异。因此，我们使用单个奇异向量重建数据，并使用image函数将原始和重建的数据集并排展示，以查看它们之间的差异。下一个图显示这两个图像非常相似。

参考阅读

  • 如我们之前提到的，PCA 可以被视为 SVD 的一个特例。在这里，我们从 SVD 生成的swiss数据的正交向量，并从prcomp获得旋转。我们可以看到这两个生成的矩阵是相同的：

    > svd.m = svd(scale(swiss))
> svd.m$v
 [,1]        [,2]         [,3]       [,4]        [,5]
[1,]  0.52396452 -0.25834215  0.003003672 -0.8090741  0.06411415
[2,] -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
[3,] -0.49150243  0.19028476  0.539337412 -0.3321615  0.56656945
[4,]  0.38530580  0.36956307  0.725888143  0.1007965 -0.42176895
[5,]  0.09167606  0.87197641 -0.424976789 -0.2154928  0.06488642
> pca.m = prcomp(swiss,scale=TRUE)
> pca.m$rotation
 PC1         PC2          PC3        PC4         PC5
Agriculture      0.52396452 -0.25834215  0.003003672 -0.8090741  0.06411415
Examination  -0.57185792 -0.01145981 -0.039840522 -0.4224580 -0.70198942
Education       -0.49150243  0.19028476  0.539337412 -0.3321615  0.56656945
Catholic          0.38530580  0.36956307  0.725888143  0.1007965 -0.42176895
Infant.Mortality 0.09167606 0.87197641 -0.424976789 -0.2154928 0.06488642

使用 SVD 压缩图像

在上一个食谱中，我们展示了如何使用 SVD 分解矩阵，然后通过乘以分解矩阵来重建数据集。此外，矩阵分解的应用可以应用于图像压缩。在这个食谱中，我们将展示如何对经典图像处理材料 Lenna 执行 SVD。

准备工作

在这个食谱中，你应该事先下载 Lenna 的图像（参考www.ece.rice.edu/~wakin/images/lena512.bmp），或者你可以准备你自己的图像来查看图像压缩是如何工作的。

如何操作...

执行以下步骤以使用 SVD 压缩图像：

  1. 首先，安装并加载bmp：

    > install.packages("bmp")
> library(bmp)

  2. 你可以使用read.bmp函数将 Lenna 的图像读取为一个数值矩阵。当读取器下载图像时，默认名称是lena512.bmp：

    > lenna = read.bmp("lena512.bmp")

  3. 旋转并绘制图像：

    > lenna = t(lenna)[,nrow(lenna):1]
> image(lenna)

    如何操作...

    Lenna 的图片

  4. 接下来，你可以对读取的数值矩阵执行 SVD 并绘制解释方差的百分比：

    > lenna.svd = svd(scale(lenna))
> plot(lenna.svd$d²/sum(lenna.svd$d²), type="l", xlab=" Singular vector", ylab = "Variance explained")

    如何操作...

    解释方差的百分比

  5. 接下来，你可以获得重建图像所需的维度数：

    > length(lenna.svd$d)
[1] 512

  6. 获得奇异向量可以解释超过 90%方差的点：

    > min(which(cumsum(lenna.svd$d²/sum(lenna.svd$d²))> 0.9))
[1] 18

  7. 你也可以将代码封装成一个函数，lenna_compression，然后你可以使用这个函数来绘制压缩后的 Lenna：

    > lenna_compression = function(dim){
+     u=as.matrix(lenna.svd$u[, 1:dim])
+     v=as.matrix(lenna.svd$v[, 1:dim])
+     d=as.matrix(diag(lenna.svd$d)[1:dim, 1:dim])
+     image(u%*%d%*%t(v))
+ }

  8. 此外，你可以使用 18 个向量来重建图像：

    > lenna_compression(18)

    如何操作...

    由 18 个组件重建的图像

  9. 你可以获得奇异向量可以解释超过 99%方差的点；

    > min(which(cumsum(lenna.svd$d²/sum(lenna.svd$d²))> 0.99))
[1] 92
> lenna_compression(92)

    如何操作...

    由 92 个组件重建的图像

它是如何工作的...

在本菜谱中，我们展示了如何使用 SVD 压缩图像。在第一步中，我们使用bmp包将图像 Lenna 加载到 R 会话中。然后，由于读取的图像已旋转，我们可以将图像旋转回来并使用plot函数在 R 中绘制 Lenna（如图 3 步中的图所示）。接下来，我们对图像矩阵进行 SVD 分解，然后绘制与奇异向量数量相关的方差解释百分比。

此外，我们发现可以使用 18 个成分解释 90%的方差，因此我们使用这 18 个成分来重建 Lenna。因此，我们创建了一个名为lenna_compression的函数，其目的是通过矩阵乘法重建图像。结果，我们将 18 作为函数的输入，返回一个相当模糊的 Lenna 图像（如图 8 步中的图所示）。然而，我们至少可以看到图像的轮廓。为了获得更清晰的图像，我们发现可以使用 92 个成分解释 99%的方差。因此，我们将函数lenna_compression的输入设置为 92。第 9 步中的图显示，这比仅使用 18 个成分构建的图像更清晰。

参见

  • Lenna 图片是压缩算法最广泛使用的标准测试图像之一。有关 Lenna 图片的更多详细信息，请参阅www.cs.cmu.edu/~chuck/lennapg/。

使用 ISOMAP 进行非线性降维

ISOMAP 是流形学习的一种方法，它将线性框架推广到非线性数据结构。与 MDS 类似，ISOMAP 创建了一系列对象相似性或差异性（距离）的视觉表示。然而，由于数据以非线性格式结构化，MDS 的欧几里得距离度量在 ISOMAP 中被数据流形的地形距离所取代。在本菜谱中，我们将说明如何使用 ISOMAP 进行非线性降维。

准备工作

在本菜谱中，我们将使用来自RnavGraphImageData的digits数据作为我们的输入源。

如何操作...

执行以下步骤以使用 ISOMAP 进行非线性降维：

  1. 首先，安装并加载RnavGraphImageData和vegan包：

    > install.packages("RnavGraphImageData")
> install.packages("vegan")
> library(RnavGraphImageData)
> library(vegan)

  2. 然后，你可以加载数据集digits：

    > data(digits)

  3. 旋转并绘制图像：

    > sample.digit = matrix(digits[,3000],ncol = 16, byrow=FALSE)
> image(t(sample.digit)[,nrow(sample.digit):1])

    如何操作...

    来自数字数据集的一个示例图像

  4. 接下来，你可以从总体中随机抽取 300 个数字：

    > set.seed(2)
> digit.idx = sample(1:ncol(digits),size = 600)
> digit.select = digits[,digit.idx]

  5. 将选定的数字数据转置，然后使用vegdist计算对象之间的差异：

    > digits.Transpose = t(digit.select)
> digit.dist = vegdist(digits.Transpose, method="euclidean")

  6. 接下来，你可以使用isomap进行降维：

    > digit.isomap = isomap(digit.dist,k = 8, ndim=6, fragmentedOK = TRUE)
> plot(digit.isomap)

    如何操作...

    ISOMAP 对象的一个二维散点图

  7. 最后，你可以将散点图与红色标记的最小生成树叠加；

    > digit.st = spantree(digit.dist)
> digit.plot = plot(digit.isomap, main="isomap k=8")
> lines(digit.st, digit.plot, col="red")

    如何操作...

    二维散点图与最小生成树的叠加

它是如何工作的...

ISOMAP 是一种非线性降维方法，也是等距映射方法的代表。ISOMAP 可以被视为度量 MDS 的扩展，其中数据点之间的成对欧几里得距离被由邻域图诱导的测地距离所取代。

ISOMAP 算法的描述分为四个步骤。首先，确定每个点的邻居。其次，构建一个邻域图。第三，计算两个节点之间的最短距离路径。最后，通过执行 MDS 找到数据的一个低维嵌入。

在这个菜谱中，我们展示了如何使用 ISOMAP 进行非线性降维。首先，我们从 RnavGraphImageData 加载数字数据。然后，在我们选择一个数字并绘制其旋转图像后，我们可以看到一个手写数字的图像（第 3 步中的图中的数字 3）。

接下来，我们随机抽取 300 个数字作为我们的输入数据给 ISOMAP。然后我们将数据集转置以计算每个图像对象之间的距离。一旦数据准备就绪，我们计算每个对象之间的距离并进行降维。在这里，我们使用 vegdist 通过欧几里得度量计算每个对象之间的不相似性。然后我们使用 ISOMAP 对 digits 数据进行非线性降维，维度设置为 6，保留一个点的最短不相似性数量为 8，并确保通过指定 fragmentedOK 为 TRUE 来分析最大的连通组。

最后，我们可以使用生成的 ISOMAP 对象制作一个二维散点图（第 6 步中的图），并在散点图上用红色线条叠加最小生成树（第 7 步中的图）。

还有更多...

你还可以使用 RnavGraph 包通过图形作为导航基础设施来可视化高维数据（在这个例子中是数字）。有关更多信息，请参阅 www.icesi.edu.co/CRAN/web/packages/RnavGraph/vignettes/RnavGraph.pdf。

这里是关于如何使用 RnavGraph 在图中可视化高维数据的描述：

  1. 首先，安装并加载 RnavGraph 和 graph 包：

    > install.packages("RnavGraph")
> source("http://bioconductor.org/biocLite.R")
> biocLite("graph")
> library(RnavGraph)

  2. 然后，你可以从 digit 数据创建一个 NG_data 对象：

    > digit.group = rep(c(1:9,0), each = 1100)
> digit.ng_data = ng_data(name = "ISO_digits",
+ data = data.frame(digit.isomap$points),
+ shortnames = paste('i',1:6, sep = ''),
+ group = digit.group[digit.idx],
+ labels = as.character(digits.group[digit.idx]))

  3. 从 NG_data 创建一个 NG_graph 对象：

    >  V = shortnames(digit.ng_data)
>  G = completegraph(V)
>  LG =linegraph(G)
>  LGnot = complement(LG)
>  ng.LG = ng_graph(name = "3D Transition", graph = LG)
> ng.LGnot = ng_graph(name = "4D Transition", graph = LGnot)

  4. 最后，你可以在 tk2d 绘图中可视化图形：

    > ng.i.digits = ng_image_array_gray('USPS Handwritten Digits',
+ digit.select,16,16,invert = TRUE,
+ img_in_row = FALSE)
> vizDigits1 = ng_2d(data = digit.ng_data, graph = ng.LG, images = ng.i.digits)
> vizDigits2 = ng_2d(data = digit.ng_data, graph = ng.LGnot, images = ng.i.digits)
> nav = navGraph(data = digit.ng_data, graph = list(ng.LG, ng.LGnot), viz = list(vizDigits1, vizDigits2))

    还有更多...

    一个三维转换图

  5. 还可以查看一个四维转换图：还有更多...

    一个四维转换图

使用局部线性嵌入进行非线性降维

局部线性嵌入（LLE）是 PCA 的扩展，它将嵌入在高维空间中的流形上的数据降低到低维空间。与 ISOMAP 不同，ISOMAP 是一种用于非线性降维的全局方法，LLE 是一种局部方法，它使用 k 近邻的线性组合来保留数据的局部属性。在本菜谱中，我们将简要介绍如何在 s 曲线数据上使用 LLE。

准备工作

在本菜谱中，我们将使用 lle 包中的 lle_scurve_data 中的数字数据作为我们的输入源。

如何操作...

执行以下步骤以使用 LLE 进行非线性降维：

  1. 首先，您需要安装并加载包，lle：

    > install.packages("lle")
> library(lle)

  2. 您可以从 lle 中加载 ll_scurve_data：

    > data( lle_scurve_data )

  3. 然后，对 lle_scurve_data 执行 lle：

    > X = lle_scurve_data
> results = lle( X=X , m=2, k=12,  id=TRUE)
finding neighbours
calculating weights
intrinsic dim: mean=2.47875, mode=2
computing coordinates

  4. 使用 str 和 plot 函数检查结果：

    > str( results )
List of 4
 $ Y     : num [1:800, 1:2] -1.586 -0.415 0.896 0.513 1.477 ...
 $ X     : num [1:800, 1:3] 0.955 -0.66 -0.983 0.954 0.958 ...
 $ choise: NULL
 $ id    : num [1:800] 3 3 2 3 2 2 2 3 3 3 ...
>plot( results$Y, main="embedded data", xlab=expression(y[1]), ylab=expression(y[2]) )

    如何操作...

    嵌入数据的 2-D 散点图

  5. 最后，您可以使用 plot_lle 来绘制 LLE 结果：

    > plot_lle( results$Y, X, FALSE, col="red", inter=TRUE )

    如何操作...

    LLE 结果的 LLE 图

如何工作...

LLE 是一种非线性降维方法，它计算高维数据的低维、邻域保留嵌入。LLE 的算法可以按以下步骤说明：首先，LLE 计算每个数据点的 k 个邻居，如何工作...。其次，它为每个点计算一组权重，这些权重最小化了残差误差之和，可以从其邻居中最佳地重建每个数据点。残差误差之和可以描述为 如何工作...，其中 如何工作... 如果 如何工作... 不是 如何工作... 的 k 个最近邻之一，并且对于每个 i，如何工作...。最后，找到由权重 W 最佳重建的向量 Y。成本函数可以表示为 如何工作...，约束条件为 如何工作...，和 如何工作...。

在本菜谱中，我们展示了如何使用 LLE 进行非线性降维。首先，我们从 lle 中加载 lle_scurve_data。然后，我们使用两个维度和 12 个邻居执行 lle，并通过指定 id = TRUE 列出每个数据点的维度。LLE 有三个步骤，包括：为数据中的每个点构建邻域，找到在该邻域中线性逼近数据的权重，以及找到低维坐标。

接下来，我们可以使用 str 和 plot 函数来检查数据。str 函数返回 X、Y、choice 和 ID。在这里，X 代表输入数据，Y 表示嵌入数据，choice 表示保留数据的索引向量，而子集选择和 ID 显示每个数据输入的维度。plot 函数返回嵌入数据的散点图。最后，我们使用 plot_lle 来绘制结果。在这里，我们通过将交互模式设置为 TRUE 来启用交互模式。

参见

  • 另一个用于非线性降维的有用软件包是 RDRToolbox，这是一个包含 ISOMAP 和 LLE 的非线性降维软件包。您可以使用以下命令安装 RDRToolbox：

    > source("http://bioconductor.org/biocLite.R")
> biocLite("RDRToolbox")
> library(RDRToolbox)

第十二章：大数据分析（R 和 Hadoop）

在本章中，我们将涵盖以下主题：

  • 准备 RHadoop 环境

  • 安装 rmr2

  • 安装 rhdfs

  • 使用 rhdfs 操作 HDFS

  • 使用 RHadoop 实现词频统计问题

  • 比较 R MapReduce 程序和标准 R 程序的性能

  • 测试和调试 rmr2 程序

  • 安装 plyrmr

  • 使用 plyrmr 操作数据

  • 使用 RHadoop 进行机器学习

  • 在 Amazon EMR 上配置 RHadoop 集群

简介

RHadoop 是一组 R 包，使用户能够使用 Hadoop 处理和分析大数据。在了解如何设置 RHadoop 并将其付诸实践之前，我们必须知道为什么我们需要使用机器学习来处理大数据规模。

在前面的章节中，我们提到了 R 在进行数据分析和机器学习时的有用性。在传统的统计分析中，重点是分析历史样本（小数据），这可能会忽略很少发生但很有价值的事件和结果，导致不确定的结论。

云技术的出现使得客户与业务之间的实时互动变得更加频繁；因此，机器学习的重点现在已转向为各种客户开发准确的预测。例如，企业可以通过使用实时预测模型，根据个人行为提供实时个性化推荐或在线广告。

然而，如果数据（例如，所有在线用户的行为）太大，无法适应单台机器的内存，你就不得不使用超级计算机或其他可扩展的解决方案。最流行的可扩展大数据解决方案是 Hadoop，它是一个开源框架，能够在集群之间存储和执行并行计算。因此，你可以使用 RHadoop，它允许 R 利用 Hadoop 的可扩展性，帮助处理和分析大数据。在 RHadoop 中，有五个主要包，它们是：

  • rmr：这是 R 和 Hadoop MapReduce 之间的接口，通过调用 Hadoop 流式 MapReduce API 在 Hadoop 集群中执行 MapReduce 作业。要开发 R MapReduce 程序，你只需要关注 map 和 reduce 函数的设计，其余的可扩展性问题将由 Hadoop 本身处理。

  • rhdfs：这是 R 和 HDFS 之间的接口，通过调用 HDFS API 来访问存储在 HDFS 中的数据。使用rhdfs的方式与使用 Hadoop shell 非常相似，它允许用户从 R 控制台轻松地操作 HDFS。

  • rhbase：这是 R 和 HBase 之间的接口，通过 Thrift 服务器在集群中访问 HBase。你可以使用rhbase来读写数据并操作存储在 HBase 中的表。

  • plyrmr：这是 MapReduce 的高级抽象，允许用户以 plyr-like 语法执行常见的数据操作。这个包大大降低了大数据操作的学习曲线。

  • ravro：这允许用户在 R 中读取avro文件或写入avro文件。它允许 R 与 HDFS 交换数据。

在本章中，我们将首先准备 Hadoop 环境，以便你可以安装 RHadoop。然后，我们将介绍三个主要包的安装：rmr、rhdfs和plyrmr。接下来，我们将介绍如何使用rmr从 R 执行 MapReduce，通过rhdfs操作 HDFS 文件，并使用plyrmr执行常见的数据操作。进一步，我们将探讨如何使用 RHadoop 进行机器学习。最后，我们将介绍如何在 Amazon EC2 上部署多个 RHadoop 集群。

准备 RHadoop 环境

由于 RHadoop 需要一个 R 和 Hadoop 集成环境，我们必须首先准备一个安装了 R 和 Hadoop 的环境。我们不必构建一个新的 Hadoop 系统，可以使用Cloudera QuickStart VM（该 VM 免费），它包含一个节点 Apache Hadoop 集群和 R。在这个菜谱中，我们将演示如何下载 Cloudera QuickStart VM。

准备工作

要使用 Cloudera QuickStart VM，建议你准备一个 64 位虚拟机操作系统，安装 VMWare 或 VirtualBox，或者安装 KVM。

如果你选择使用 VMWare，你应该准备一个与 WorkStation 8.x 或更高版本兼容的播放器：4.x 或更高版本，ESXi 5.x 或更高版本，或者 Fusion 4.x 或更高版本。

注意，启动 VM 需要 4 GB 的 RAM，至少有 3 GB 的可用磁盘空间。

如何操作...

使用 Cloudera QuickStart VM 设置 Hadoop 环境的以下步骤：

  1. 访问 Cloudera QuickStart VM 下载站点（你可能需要更新链接，因为 Cloudera 升级了其 VMs，当前 CDH 版本为 5.3），请参阅www.cloudera.com/content/cloudera/en/downloads/quickstart_vms/cdh-5-3-x.html。如何操作...

    Cloudera QuickStart VM 下载站点的截图

  2. 根据你操作系统上安装的虚拟机平台，选择适当的链接（你可能需要更新链接，因为 Cloudera 升级了其 VMs）以下载 VM 文件：

  3. 接下来，您可以使用安装在您的操作系统上的虚拟机平台启动 QuickStart 虚拟机。您应该在几分钟内看到 Centos 6.2 的桌面。如何操作...

    Cloudera QuickStart 虚拟机的截图。

  4. 您可以打开一个终端并输入 hadoop，这将显示可以操作 Hadoop 集群的一组功能。如何操作...

    输入 hadoop 后的终端截图

  5. 打开一个终端并输入 R。访问 R 会话并检查版本 3.1.1 是否已在 Cloudera QuickStart 虚拟机中安装。如果您在虚拟机中找不到已安装的 R，请使用以下命令安装 R：

    $ yum install R R-core R-core-devel R-devel

它是如何工作的...

您不必自己构建 Hadoop 系统，可以使用 Cloudera 提供的 Hadoop VM 应用程序（虚拟机是免费的）。QuickStart 虚拟机在 CentOS 6.2 上运行，包含单个节点 Apache Hadoop 集群、Hadoop 生态系统模块和已安装的 R。这可以帮助您节省时间，而不是需要您学习如何安装和使用 Hadoop。

QuickStart 虚拟机要求您拥有一个具有 64 位客户操作系统的计算机，至少 4 GB 的 RAM、3 GB 的磁盘空间，并且已安装 VMWare、VirtualBox 或 KVM。因此，您可能无法在某些计算机上使用此版本的虚拟机。作为替代方案，您可以考虑使用 Amazon 的 Elastic MapReduce。我们将在本章的最后一个小节中说明如何在 EMR 中准备 RHadoop 环境。

设置 Cloudera QuickStart 虚拟机很简单。从下载站点下载虚拟机，然后使用 VMWare、VirtualBox 或 KVM 打开构建的镜像。一旦您可以看到 CentOS 的桌面，您就可以访问终端并输入 hadoop 来查看 Hadoop 是否正在运行；然后，输入 R 来查看 R 是否在 QuickStart 虚拟机中运行。

参见

安装 rmr2

rmr2 包允许您通过 Hadoop 集群上的 MapReduce 执行大数据处理和分析。要在 Hadoop 集群上执行 MapReduce，您必须在每个任务节点上安装 R 和 rmr2。在本菜谱中，我们将说明如何在 Hadoop 集群的单个节点上安装 rmr2。

准备工作

确保您已完成了前面的菜谱，通过启动 Cloudera QuickStart 虚拟机并将其连接到互联网，以便您可以继续下载和安装 rmr2 包。

如何操作...

执行以下步骤在 QuickStart 虚拟机上安装 rmr2：

  1. 首先，在 Cloudera QuickStart 虚拟机内部打开终端。

  2. 使用 root 权限进入 R 会话：

    $ sudo R

  3. 您可以在安装 rmr2 之前安装依赖包：

    > install.packages(c("codetools", "Rcpp", "RJSONIO", "bitops", "digest", "functional", "stringr", "plyr", "reshape2", "rJava", "caTools"))

  4. 退出 R 会话：

    > q()

  5. 接下来，您可以将 rmr-3.3.0 下载到 QuickStart VM 上。如果 Revolution Analytics 升级了 rmr2 的版本，您可能需要更新链接：

    $ wget --no-check-certificate https://raw.githubusercontent.com/RevolutionAnalytics/rmr2/3.3.0/build/rmr2_3.3.0.tar.gz

  6. 然后，您可以将 rmr-3.3.0 安装到 QuickStart VM 上：

    $ sudo R CMD INSTALL rmr2_3.3.0.tar.gz

  7. 最后，您可以进入 R 会话并使用 library 函数来测试库是否已成功安装：

    $ R
> library(rmr2)

它是如何工作的...

为了在 Hadoop 集群上执行 MapReduce，您必须在每个任务节点上安装 R 和 RHadoop。在此，我们将说明如何在 Hadoop 集群的单个节点上安装 rmr2。首先，打开 Cloudera QuickStart VM 的终端。在安装 rmr2 之前，我们首先以 root 权限访问 R 会话并安装依赖的 R 包。

接下来，在所有依赖包安装完成后，退出 R 会话，并在 Linux shell 中使用 wget 命令从 GitHub 下载 rmr-3.3.0 到本地文件系统。然后，您可以开始安装 rmr2。最后，您可以通过 R 会话使用库函数来验证包是否已安装。

参考信息

安装 rhdfs

rhdfs 包是 R 和 HDFS 之间的接口，它允许用户从 R 控制台访问 HDFS。类似于 rmr2，应该在每个任务节点上安装 rhdfs，以便可以通过 R 控制台访问 HDFS 资源。在本菜谱中，我们将介绍如何在 Cloudera QuickStart VM 上安装 rhdfs。

准备工作

确保您已通过启动 Cloudera QuickStart VM 并将其连接到互联网来完成前面的菜谱，这样您就可以继续下载和安装 rhdfs 包。

如何操作...

执行以下步骤来安装 rhdfs：

  1. 首先，您可以从 GitHub 下载 rhdfs 1.0.8。如果 Revolution Analytics 升级了 rhdfs 的版本，您可能需要更新链接：

    $wget --no-check-certificate https://raw.github.com/RevolutionAnalytics/rhdfs/master/build/rhdfs_1.0.8.tar.gz

  2. 接下来，您可以在命令行模式下安装 rhdfs：

    $ sudo HADOOP_CMD=/usr/bin/hadoop  R CMD INSTALL rhdfs_1.0.8.tar.gz

  3. 然后，您可以设置 JAVA_HOME。JAVA_HOME 的配置取决于 VM 内安装的 Java 版本：

    $ sudo JAVA_HOME=/usr/java/jdk1.7.0_67-cloudera R CMD javareconf

  4. 最后，您可以设置系统环境和初始化 rhdfs。如果您使用的是不同版本的 QuickStart VM，可能需要更新环境设置：

    $ R
> Sys.setenv(HADOOP_CMD="/usr/bin/hadoop")
> Sys.setenv(HADOOP_STREAMING="/usr/lib/hadoop-mapreduce/hadoop-
streaming-2.5.0-cdh5.2.0.jar")
> library(rhdfs)
> hdfs.init()

它是如何工作的...

包 rhdfs 提供了函数，使用户可以使用 R 管理 HDFS。类似于 rmr2，您应该在每个任务节点上安装 rhdfs，以便可以通过 R 控制台访问 HDFS。

要安装 rhdfs，您首先应从 GitHub 下载 rhdfs。然后，您可以通过指定 HADOOP_CMD 的位置来在 R 中安装 rhdfs。您必须通过命令 javareconf 配置 R 以支持 Java。

接下来，您可以访问 R 并配置 HADOOP_CMD 和 HADOOP_STREAMING 的位置。最后，您可以通过 rhdfs.init 函数初始化 rhdfs，这允许您通过 rhdfs 开始操作 HDFS。

参见

  • 要查找 HADOOP_CMD 的位置，您可以在 Linux shell 中使用 which hadoop 命令。在大多数 Hadoop 系统中，HADOOP_CMD 位于 /usr/bin/hadoop。

  • 对于 HADOOP_STREAMING 的位置，流 JAR 文件通常位于 /usr/lib/hadoop-mapreduce/。但是，如果您在 Linux 系统中找不到 /usr/lib/Hadoop-mapreduce 目录，您可以使用 locate 命令搜索流 JAR。例如：

    $ sudo updatedb
$ locate streaming | grep jar | more

使用 rhdfs 操作 HDFS

rhdfs 包是 Hadoop 和 R 之间的接口，可以在后端调用 HDFS API 来操作 HDFS。因此，您可以通过使用 rhdfs 包轻松地从 R 控制台操作 HDFS。在下面的菜谱中，我们将演示如何使用 rhdfs 函数操作 HDFS。

准备工作

要继续本菜谱，您需要完成之前的菜谱，将 rhdfs 安装到 R 中，并验证您可以通过 hdfs.init 函数初始化 HDFS。

如何做...

执行以下步骤以操作存储在 HDFS 上的文件：

  1. 初始化 rhdfs 包：

    > Sys.setenv(HADOOP_CMD="/usr/bin/hadoop")
> Sys.setenv(HADOOP_STREAMING="/usr/lib/hadoop-mapreduce/hadoop-streaming-2.5.0-cdh5.2.0.jar")
> library(rhdfs)
> hdfs.init ()

  2. 然后，您可以按照以下方式操作存储在 HDFS 上的文件：

    • hdfs.put: 从本地文件系统复制文件到 HDFS：

      > hdfs.put('word.txt', './')

    • hdfs.ls: 读取 HDFS 中的目录列表：

      > hdfs.ls('./')

    • hdfs.copy: 将文件从一个 HDFS 目录复制到另一个目录：

      > hdfs.copy('word.txt', 'wordcnt.txt')

    • hdfs.move：将文件从一个 HDFS 目录移动到另一个目录：

      > hdfs.move('wordcnt.txt', './data/wordcnt.txt')

    • hdfs.delete: 从 R 删除 HDFS 目录：

      > hdfs.delete('./data/')

    • hdfs.rm: 从 R 删除 HDFS 目录：

      > hdfs.rm('./data/')

    • hdfs.get: 从 HDFS 下载文件到本地文件系统：

      > hdfs.get(word.txt', '/home/cloudera/word.txt')

    • hdfs.rename: 重命名 HDFS 上存储的文件：

      hdfs.rename('./test/q1.txt','./test/test.txt')

    • hdfs.chmod: 更改文件或目录的权限：

      > hdfs.chmod('test', permissions= '777')

    • hdfs.file.info: 读取 HDFS 文件的元信息：

      > hdfs.file.info('./')

  3. 此外，您还可以将流写入 HDFS 文件：

    > f = hdfs.file("iris.txt","w")
> data(iris)
> hdfs.write(iris,f)
> hdfs.close(f)

  4. 最后，您还可以从 HDFS 文件中读取流：

    > f = hdfs.file("iris.txt", "r")
> dfserialized = hdfs.read(f)
> df = unserialize(dfserialized)
> df
> hdfs.close(f)

工作原理...

在本菜谱中，我们演示如何使用 rhdfs 包操作 HDFS。通常，您可以使用 Hadoop shell 操作 HDFS，但如果您想从 R 访问 HDFS，则可以使用 rhdfs 包。

在开始使用 rhdfs 之前，您必须使用 hdfs.init() 初始化 rhdfs。初始化后，您可以通过 rhdfs 包中提供的函数操作 HDFS。

除了操作 HDFS 文件外，您还可以通过 hdfs.read 和 hdfs.write 将流交换到 HDFS。因此，我们将演示如何使用 hdfs.write 将 R 中的数据框写入 HDFS 文件，iris.txt。最后，您可以使用 hdfs.read 函数和 unserialize 函数将写入的文件恢复到数据框。

参见

  • 要初始化rhdfs，您必须在系统环境中设置HADOOP_CMD和HADOOP_STREAMING。您不必每次使用rhdfs时都设置配置，可以将配置放在.rprofile文件中。因此，每次您启动一个 R 会话时，配置将自动加载。

使用 RHadoop 实现单词计数问题

为了演示 MapReduce 是如何工作的，我们以单词计数为例，它计算给定输入集中每个单词的出现次数。在这个菜谱中，我们将演示如何使用rmr2实现单词计数问题。

准备工作

在这个菜谱中，我们需要一个输入文件作为我们的单词计数程序输入。您可以从github.com/ywchiu/ml_R_cookbook/tree/master/CH12下载示例输入文件。

如何做...

执行以下步骤以实现单词计数程序：

  1. 首先，您需要配置系统环境，然后将rmr2和rhdfs加载到一个 R 会话中。如果您使用的是不同版本的 QuickStart VM，可能需要更新 JAR 文件的使用：

    > Sys.setenv(HADOOP_CMD="/usr/bin/hadoop")
> Sys.setenv(HADOOP_STREAMING="/usr/lib/hadoop-mapreduce/hadoop-streaming-2.5.0-cdh5.2.0.jar ")
> library(rmr2)
> library(rhdfs)
> hdfs.init()

  2. 您可以在 HDFS 上创建一个目录并将输入文件放入新创建的目录中：

    > hdfs.mkdir("/user/cloudera/wordcount/data")
> hdfs.put("wc_input.txt", "/user/cloudera/wordcount/data")

  3. 接下来，您可以创建一个map函数：

    > map = function(.,lines) { keyval(
+   unlist(
+     strsplit(
+       x = lines, 
+       split = " +")),
+   1)}

  4. 创建一个reduce函数：

    > reduce = function(word, counts) { 
+   keyval(word, sum(counts)) 
+ }

  5. 调用 MapReduce 程序来计数文档中的单词：

    > hdfs.root = 'wordcount' > hdfs.data = file.path(hdfs.root, 'data')
> hdfs.out = file.path(hdfs.root, 'out')
> wordcount = function (input, output=NULL) { 
+  mapreduce(input=input, output=output, input.format="text", map=map, 
+  reduce=reduce) 
+ } 
> out = wordcount(hdfs.data, hdfs.out)

  6. 最后，您可以检索文档中出现的最频繁的 10 个单词：

    > results = from.dfs(out) 
> results$key[order(results$val, decreasing = TRUE)][1:10]

它是如何工作的...

在这个菜谱中，我们演示了如何使用rmr2包实现单词计数。首先，我们需要配置系统环境并将rhdfs和rmr2加载到 R 中。然后，我们通过hdfs.put函数将我们的单词计数程序的输入从本地文件系统指定到 HDFS 目录/user/cloudera/wordcount/data。

接下来，我们开始实现 MapReduce 程序。通常，我们可以将 MapReduce 程序分为map和reduce函数。在map函数中，我们首先使用strsplit函数将每一行拆分成单词。然后，由于strsplit函数返回一个单词列表，我们可以使用unlist函数将字符向量转换为列表。最后，我们可以返回键值对，其中每个单词作为键，值为一。由于reduce函数接收来自map函数生成的键值对，reduce函数将计数求和并返回每个单词（或键）的出现次数。

在我们实现了map和reduce函数之后，我们可以通过mapreduce函数提交我们的作业。通常，mapreduce函数需要四个输入，分别是 HDFS 输入路径、HDFS 输出路径、map 函数和 reduce 函数。在这种情况下，我们指定输入为wordcount/data，输出为wordcount/out，mapfunction 为map，reduce function 为reduce，并将mapreduce调用封装在wordcount函数中。最后，我们调用函数wordcount并将输出路径存储在变量out中。

我们可以使用from.dfs函数将 HDFS 数据加载到results变量中，该变量包含单词和出现次数的映射。然后我们可以从results变量中生成出现次数最多的前 10 个单词。

参见

比较 R MapReduce 程序和标准 R 程序的性能

对于不熟悉 Hadoop 工作原理的人来说，他们可能会经常将 Hadoop 视为大数据处理的一种补救措施。有些人可能认为 Hadoop 可以在几毫秒内为任何大小的数据返回处理结果。在这个菜谱中，我们将比较 R MapReduce 程序和标准 R 程序的性能，以证明 Hadoop 的运行速度并不像一些人认为的那样快。

准备工作

在这个菜谱中，你应该通过将rmr2安装到 R 环境中来完成之前的菜谱。

如何做...

执行以下步骤来比较标准 R 程序和 R MapReduce 程序的性能：

  1. 首先，你可以实现一个标准 R 程序来将所有数字平方：

    > a.time = proc.time() 
> small.ints2=1:100000 
> result.normal = sapply(small.ints2, function(x) x²) 
> proc.time() - a.time

  2. 为了比较性能，你可以实现一个 R MapReduce 程序来将所有数字平方：

    > b.time = proc.time() 
> small.ints= to.dfs(1:100000) 
> result = mapreduce(input = small.ints, map = function(k,v)       cbind(v,v²)) 
> proc.time() - b.time

它是如何工作的...

在这个菜谱中，我们实现了两个程序来平方所有数字。在第一个程序中，我们使用标准的 R 函数``sapply来平方从 1 到 100,000 的序列。为了记录程序执行时间，我们在执行前首先记录处理时间到a.time，然后从执行后的当前处理时间中减去a.time。通常，执行时间不超过 10 秒。在第二个程序中，我们使用rmr2`包在 R MapReduce 版本中实现程序。在这个程序中，我们也记录了执行时间。通常，这个程序需要几分钟来完成一个任务。

性能比较显示，当处理少量数据时，标准 R 程序的性能优于 MapReduce 程序。这是因为 Hadoop 系统通常需要时间来启动守护进程，守护进程之间的作业协调，以及从数据节点获取数据。因此，MapReduce 程序通常需要几分钟到几个小时来完成执行。因此，如果你可以将你的数据放入内存中，你应该编写一个标准 R 程序来解决问题。否则，如果数据太大而无法放入内存，你可以实现一个 MapReduce 解决方案。

参见

  • 为了检查一个作业在 Hadoop 中是否能够顺利高效地运行，你可以运行一个 MapReduce 基准测试，MRBench，来评估作业的性能：

    $ hadoop jar /usr/lib/hadoop-0.20-mapreduce/hadoop-test.jar mrbench -numRuns 50

测试和调试 rmr2 程序

由于运行 MapReduce 程序将需要相当多的时间，从几分钟到几小时不等，因此测试和调试变得非常重要。在这个菜谱中，我们将展示您可以用来调试 R MapReduce 程序的一些技术。

准备工作

在这个菜谱中，您应该通过将 rmr2 安装到 R 环境中来完成前面的菜谱。

如何操作...

执行以下步骤以测试和调试 R MapReduce 程序：

  1. 首先，您可以在 rmr.options 中将后端配置为本地：

    > rmr.options(backend = 'local')

  2. 再次，您可以执行前面菜谱中提到的数字平方 MapReduce 程序：

    > b.time = proc.time() 
> small.ints= to.dfs(1:100000) 
> result = mapreduce(input = small.ints, map = function(k,v)       cbind(v,v²)) 
> proc.time() - b.time

  3. 此外，如果您想在 MapReduce 程序中打印任何变量的结构信息，您可以使用 rmr.str 函数：

    > out = mapreduce(to.dfs(1), map = function(k, v) rmr.str(v))
Dotted pair list of 14
 $ : language mapreduce(to.dfs(1), map = function(k, v) rmr.str(v))
 $ : language mr(map = map, reduce = reduce, combine = combine, vectorized.reduce, in.folder = if (is.list(input)) {     lapply(input, to.dfs.path) ...
 $ : language c.keyval(do.call(c, lapply(in.folder, function(fname) {     kv = get.data(fname) ...
 $ : language do.call(c, lapply(in.folder, function(fname) {     kv = get.data(fname) ...
 $ : language lapply(in.folder, function(fname) {     kv = get.data(fname) ...
 $ : language FUN("/tmp/Rtmp813BFJ/file25af6e85cfde"[[1L]], ...)
 $ : language unname(tapply(1:lkv, ceiling((1:lkv)/(lkv/(object.size(kv)/10⁶))), function(r) {     kvr = slice.keyval(kv, r) ...
 $ : language tapply(1:lkv, ceiling((1:lkv)/(lkv/(object.size(kv)/10⁶))), function(r) {     kvr = slice.keyval(kv, r) ...
 $ : language lapply(X = split(X, group), FUN = FUN, ...)
 $ : language FUN(X[[1L]], ...)
 $ : language as.keyval(map(keys(kvr), values(kvr)))
 $ : language is.keyval(x)
 $ : language map(keys(kvr), values(kvr))
 $ :length 2 rmr.str(v)
 ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 34 1 58 34 58 1 1
 .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x3f984f0> 
v
 num 1

它是如何工作的...

在这个菜谱中，我们介绍了一些在实现 MapReduce 程序时可以使用的调试和测试技术。首先，我们介绍了在本地模式下测试 MapReduce 程序的技术。如果您想以伪分布式或完全分布式模式运行 MapReduce 程序，这将花费您几分钟到几小时的时间来完成，这将在调试您的 MapReduce 程序时造成大量时间的浪费。因此，您可以在 rmr.options 中将后端设置为本地模式，这样程序将以本地模式执行，执行时间更短。

另一种调试技术是在 map 或 reduce 函数中列出变量的内容。在 R 程序中，您可以使用 str 函数显示单个变量的紧凑结构。在 rmr2 中，该包还提供了一个名为 rmr.str 的函数，允许您将单个变量的内容打印到控制台。在这个例子中，我们使用 rmr.str 来打印 MapReduce 程序中变量的内容。

参考以下内容

  • 对于那些对 rmr2 包的 option 设置感兴趣的人，您可以参考 rmr.options 的帮助文档：

    > help(rmr.options)

安装 plyrmr

plyrmr 包为用户提供了一些常见操作（如在 plyr 或 reshape2 中找到的操作），以便用户能够通过 MapReduce 框架轻松进行数据处理。在这个菜谱中，我们将介绍如何在 Hadoop 系统上安装 plyrmr。

准备工作

通过启动 Cloudera QuickStart 虚拟机并将其连接到互联网，确保您已经完成了前面的菜谱。此外，您需要事先安装 rmr2 包。

如何操作...

执行以下步骤在 Hadoop 系统上安装 plyrmr：

  1. 首先，您应该在 Linux 命令行中安装 libxml2-devel 和 curl-devel：

    $ yum install libxml2-devel
$ sudo yum install curl-devel

  2. 然后，您可以访问 R 并安装依赖包：

    $ sudo R
> Install.packages(c(" Rcurl", "httr"),  dependencies = TRUE
> Install.packages("devtools", dependencies = TRUE)
> library(devtools)
> install_github("pryr", "hadley")
> install.packages(c(" R.methodsS3", "hydroPSO"),  dependencies = TRUE)
> q()

  3. 接下来，您可以在 Hadoop 虚拟机上下载 plyrmr 0.5.0 并安装它。如果 Revolution Analytics 升级了 plyrmr 的版本，您可能需要更新链接：

    $ wget -no-check-certificate https://raw.github.com/RevolutionAnalytics/plyrmr/master/build/plyrmr_0.5.0.tar.gz
$ sudo R CMD INSTALL plyrmr_0.5.0.tar.gz

  4. 最后，验证安装：

    $ R
> library(plyrmr)

它是如何工作的...

除了使用rmr2包编写 R MapReduce 程序外，您还可以使用plyrmr来操作数据。plyrmr包类似于 Hadoop 生态系统中的 hive 和 pig，它是 MapReduce 程序的高级抽象。因此，我们可以以plyr风格实现 R MapReduce 程序，而不是实现map和reduce函数。

要安装plyrmr，首先使用yum install命令安装libxml2-devel和curl-devel包。然后，进入 R 并安装依赖包。最后，从 GitHub 下载文件并在 R 中安装plyrmr。

参考以下内容

  • 要获取有关plyrmr的更多信息，您可以使用help函数参考以下文档：

    > help(package=plyrmr)

使用 plyrmr 操作数据

虽然使用rmr2编写 MapReduce 程序比编写原生 Java 版本要容易得多，但对于非开发者来说，编写 MapReduce 程序仍然很困难。因此，您可以使用plyrmr，它是 MapReduce 程序的高级抽象，这样您就可以使用类似于 plyr 的操作来操作大数据。在本菜谱中，我们将介绍一些您可以用来操作数据的操作。

准备工作

在本菜谱中，您应该已经通过在 R 中安装plyrmr和rmr2来完成之前的菜谱。

如何操作...

按照以下步骤使用plyrmr操作数据：

  1. 首先，您需要将plyrmr和rmr2都加载到 R 中：

    > library(rmr2)
> library(plyrmr)

  2. 然后，您可以将执行模式设置为本地模式：

    > plyrmr.options(backend="local")

  3. 接着，将泰坦尼克号数据集加载到 R 中：

    > data(Titanic)
> titanic = data.frame(Titanic)

  4. 通过过滤数据开始操作：

    > where(
+    Titanic, 
+ Freq >=100)

  5. 您还可以使用管道操作符过滤数据：

    > titanic %|% where(Freq >=100)

  6. 将泰坦尼克号数据放入 HDFS，并将数据路径加载到变量tidata中：

    > tidata = to.dfs(data.frame(Titanic), output = '/tmp/titanic')
> tidata

  7. 接下来，您可以从泰坦尼克号数据中生成频率的汇总：

    > input(tidata) %|% transmute(sum(Freq))

  8. 您还可以按性别分组频率：

    > input(tidata) %|% group(Sex) %|% transmute(sum(Freq))

  9. 您可以从总体中抽取 10 条记录作为样本：

    > sample(input(tidata), n=10)

  10. 此外，您还可以使用 plyrmr 将两个数据集连接起来：

    > convert_tb = data.frame(Label=c("No","Yes"), Symbol=c(0,1))
ctb = to.dfs(convert_tb, output = 'convert')
> as.data.frame(plyrmr::merge(input(tidata), input(ctb), by.x="Survived", by.y="Label"))
> file.remove('convert')

工作原理...

在本菜谱中，我们介绍如何使用plyrmr操作数据。首先，我们需要将plyrmr包加载到 R 中。然后，类似于rmr2，您必须将plyrmr的后端选项设置为本地模式。否则，如果plyrmr在 Hadoop 模式下运行（默认设置），您可能需要等待几分钟到几个小时。

接下来，我们可以开始使用数据过滤进行数据操作。您可以选择在步骤 4 中调用嵌套在其他函数调用中的函数。另一方面，您可以使用管道操作符%|%来链式执行多个操作。因此，我们可以在步骤 5 中使用管道操作符执行类似于步骤 4 的数据过滤。

接下来，你可以使用to.dfs根据当前运行模式将数据集输入到 HDFS 或本地文件系统。该函数将生成数据集的路径并将其保存在变量tidata中。通过知道路径，你可以使用input函数访问数据。接下来，我们将展示如何使用transmute和sum函数从泰坦尼克号数据集中生成频率的总和。此外，plyrmr允许用户按性别汇总频率。

此外，为了从总体中采样数据，你也可以使用sample函数从泰坦尼克号数据集中选择 10 条记录。最后，我们演示了如何使用plyrmr的merge函数合并两个数据集。

参见

这里我们列出了一些可以使用plyrmr来操作数据的函数。你可以参考help函数以获取它们的使用和功能的更多详细信息：

  • 数据操作：

    • bind.cols: 这用于添加新列

    • select: 这用于选择列

    • where: 这用于选择行

    • transmute: 这使用上述所有内容及其汇总

  • 从reshape2：

    • melt和dcast：它将长宽数据框进行转换

  • 摘要：

    • count

    • quantile

    • sample

  • 提取：

    • top.k

    • bottom.k

使用 RHadoop 进行机器学习

在前面的章节中，我们已经展示了当使用 R 解决机器学习问题时，R 是多么强大。我们还展示了使用 Hadoop 可以使 R 并行处理大数据。此时，有些人可能认为使用 RHadoop 可以通过现有的许多机器学习包轻松解决大数据的机器学习问题。然而，你不能使用其中大多数来解决机器学习问题，因为它们不能在 MapReduce 模式下执行。在下面的菜谱中，我们将展示如何实现线性回归的 MapReduce 版本，并将其与使用lm函数的版本进行比较。

准备工作

在这个菜谱中，你应该通过将rmr2安装到 R 环境中来完成之前的菜谱。

如何操作...

执行以下步骤以在 MapReduce 中实现线性回归：

  1. 首先，从MASS包中加载cats数据集：

    > library(MASS)
> data(cats)
> X = matrix(cats$Bwt)
> y = matrix(cats$Hwt)

  2. 你可以通过调用lm函数来生成线性回归模型：

    > model = lm(y~X)
> summary(model)

Call:
lm(formula = y ~ X)

Residuals:
 Min      1Q  Median      3Q     Max 
-3.5694 -0.9634 -0.0921  1.0426  5.1238 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)  -0.3567     0.6923  -0.515    0.607 
X             4.0341     0.2503  16.119   <2e-16 ***
---
Signif. codes: 
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.452 on 142 degrees of freedom
Multiple R-squared:  0.6466,  Adjusted R-squared:  0.6441 
F-statistic: 259.8 on 1 and 142 DF,  p-value: < 2.2e-16

  3. 你现在可以使用给定的数据点和模型制作回归图：

    > plot(y~X)
> abline(model, col="red")

    如何操作...

    cats 数据集的线性回归图

  4. 将rmr2加载到 R 中：

    > Sys.setenv(HADOOP_CMD="/usr/bin/hadoop")
> Sys.setenv(HADOOP_STREAMING="/usr/lib/hadoop-mapreduce/hadoop-> streaming-2.5.0-cdh5.2.0.jar")
> library(rmr2)
> rmr.options(backend="local")

  5. 你可以设置X和y的值：

    > X = matrix(cats$Bwt)
> X.index = to.dfs(cbind(1:nrow(X), X))
> y = as.matrix(cats$Hwt)

  6. 创建一个Sum函数来汇总值：

    > Sum = 
+   function(., YY) 
+     keyval(1, list(Reduce('+', YY)))

  7. 在 MapReduce 的Job1中计算Xtx：

    > XtX = 
+    values(
+      from.dfs(
+        mapreduce(
+          input = X.index,
+          map = 
+            function(., Xi) {
+              Xi = Xi[,-1]
+              keyval(1, list(t(Xi) %*% Xi))},
+          reduce = Sum,
+          combine = TRUE)))[[1]]

  8. 你可以在 MapReduce 的Job2中计算Xty：

    Xty = 
+    values(
+      from.dfs(
+        mapreduce(
+          input = X.index,
+          map = function(., Xi) {
+            yi = y[Xi[,1],]
+            Xi = Xi[,-1]
+            keyval(1, list(t(Xi) %*% yi))},
+          reduce = Sum,
+          combine = TRUE)))[[1]]

  9. 最后，你可以从XtX和Xty中推导出系数：

    > solve(XtX, Xty)
 [,1]
[1,] 3.907113

它是如何工作的...

在这个菜谱中，我们展示了如何在 R 中用 MapReduce 风格实现线性逻辑回归。在我们开始实现之前，我们回顾了传统线性模型是如何工作的。我们首先从MASS包中检索cats数据集。然后我们将X作为体重（Bwt）和y作为心脏重量（Hwt）加载。

接下来，我们开始使用lm函数将数据拟合到线性回归模型中。然后我们可以计算拟合模型并获得模型的摘要。摘要显示系数为 4.0341，截距为-0.3567。此外，我们根据给定的数据点绘制散点图，然后在图上绘制回归线。

由于我们不能在 MapReduce 形式中使用lm函数进行线性回归，我们必须以 MapReduce 风格重写回归模型。在这里，我们希望分三步实现线性回归的 MapReduce 版本，这些步骤是：使用 MapReduce 计算Xtx值，任务 1，使用 MapReduce 计算Xty值，任务 2，然后推导出系数值：

  • 在第一步中，我们将矩阵X作为输入传递给map函数。然后map函数计算转置矩阵X和X的叉积。然后reduce函数执行之前定义的求和操作。

  • 在第二步中，计算Xty的过程与计算XtX的过程类似。这个过程计算了转置矩阵X和y的叉积。然后reduce函数执行之前定义的求和操作。

  • 最后，我们使用solve函数推导出系数，其值为 3.907113。

如结果所示，lm和 MapReduce 计算出的系数略有差异。一般来说，lm模型计算出的系数比 MapReduce 计算出的更准确。然而，如果你的数据太大而无法适应内存，你不得不选择在 MapReduce 版本中实现线性回归。

参考信息

在 Amazon EMR 上配置 RHadoop 集群

到目前为止，我们只演示了如何在单个 Hadoop 节点上运行 RHadoop 程序。为了在多节点集群上测试我们的 RHadoop 程序，你需要做的唯一一件事是在 Hadoop 集群的所有任务节点（无论是 mapreduce 版本 1 的任务跟踪器还是 map reduce 版本 2 的节点管理器）上安装 RHadoop。然而，部署和安装是耗时的。另一方面，你可以选择在 Amazon EMR 上部署你的 RHadoop 程序，这样你就可以在几分钟内部署多节点集群和 RHadoop 到每个任务节点。在下面的菜谱中，我们将演示如何在 Amazon EMR 服务上配置 RHadoop 集群。

准备工作

在这个菜谱中，你必须注册并创建 AWS 账户，并且在使用 Amazon EMR 之前，你必须知道如何生成 EC2 密钥对。

对于那些寻求如何开始使用 AWS 的更多信息的人，请参阅亚马逊提供的教程docs.aws.amazon.com/AWSEC2/latest/UserGuide/EC2_GetStarted.html。

如何操作...

在 Amazon EMR 上配置 RHadoop 的以下步骤

  1. 首先，你可以访问亚马逊网络服务的控制台（参考us-west-2.console.aws.amazon.com/console/)，在分析部分找到 EMR。然后，点击EMR。如何操作...

    从 AWS 控制台访问 EMR 服务。

  2. 你应该发现自己处于 EMR 仪表板的集群列表中（参考us-west-2.console.aws.amazon.com/elasticmapreduce/home?region=us-west-2#cluster-list::)；点击创建集群。如何操作...

    EMR 集群列表

  3. 然后，你应该发现自己处于创建集群页面（参考us-west-2.console.aws.amazon.com/elasticmapreduce/home?region=us-west-2#create-cluster:）。

  4. 接下来，你应在集群配置中指定集群名称和日志文件夹 S3 位置。如何操作...

    创建集群页面中的集群配置

  5. 然后，你可以在软件配置中配置 Hadoop 发行版。如何操作...

    配置软件和应用程序

  6. 接下来，你可以在 Hadoop 集群内配置节点数量。如何操作...

    配置 Hadoop 集群内的硬件

  7. 然后，你可以指定用于主节点登录的 EC2 密钥对。如何操作...

    EMR 集群主节点的安全和访问

  8. 要设置 RHadoop，必须在每个任务节点上执行引导操作来安装 RHadoop。请创建一个名为 bootstrapRHadoop.sh 的文件，并在文件中插入以下行：

    echo 'install.packages(c("codetools", "Rcpp", "RJSONIO", "bitops", "digest", "functional", "stringr", "plyr", "reshape2", "rJava", "caTools"), repos="http://cran.us.r-project.org")' > /home/hadoop/installPackage.R
sudo Rscript /home/hadoop/installPackage.R
wget --no-check-certificate https://raw.githubusercontent.com/RevolutionAnalytics/rmr2/master/build/rmr2_3.3.0.tar.gz
sudo R CMD INSTALL rmr2_3.3.0.tar.gz
wget --no-check-certificate https://raw.github.com/RevolutionAnalytics/rhdfs/master/build/rhdfs_1.0.8.tar.gz
sudo HADOOP_CMD=/home/hadoop/bin/hadoop R CMD INSTALL rhdfs_1.0.8.tar.gz

  9. 你应该将 bootstrapRHadoop.sh 上传到 S3。

  10. 现在，你需要使用自定义操作添加引导操作，并在 S3 位置添加 s3://<location>/bootstrapRHadoop.sh。如何操作...

    设置引导操作

  11. 接下来，你可以点击创建集群来启动 Hadoop 集群。如何操作...

    创建集群

  12. 最后，当集群准备就绪时，你应该看到主节点的公共 DNS。现在，你可以使用你的 EC2 密钥对访问主节点的终端：如何操作...

    创建集群的截图

它是如何工作的...

在这个菜谱中，我们展示了如何在 Amazon EMR 上设置 RHadoop。这样做的好处是，您只需点击几下，在几分钟内就可以快速创建一个可扩展的、按需的 Hadoop。这有助于节省您构建和部署 Hadoop 应用程序的时间。然而，您必须为每个实例的运行时间付费。在使用 Amazon EMR 之前，您应该创建一个 AWS 账户，并了解如何设置 EC2 密钥对和 S3。然后，您可以在 Amazon EMR 上开始安装 RHadoop。

在第一步中，访问 EMR 集群列表并点击 创建集群。您可以在 创建集群 页面上看到配置列表。然后，您应该在集群配置中的 S3 位置设置集群名称和日志文件夹。

接下来，您可以设置软件配置并选择您想要安装的 Hadoop 发行版。亚马逊提供自己的发行版和 MapR 发行版。通常情况下，除非您对默认的 Hadoop 发行版有顾虑，否则您会跳过这一部分。

然后，通过指定主节点、核心节点和任务节点来配置硬件。默认情况下，只有一个主节点和两个核心节点。如果您愿意，可以添加更多核心节点和任务节点。然后，您应该设置密钥对以登录主节点。

接下来，创建一个包含所有启动脚本名为 bootstrapRHadoop.sh 的文件。文件创建后，您应该将其保存在 S3 存储中。然后，您可以在 引导操作 中指定 custom action，并将 bootstrapRHadoop.sh 作为引导脚本。最后，您可以点击 创建集群 并等待集群就绪。一旦集群就绪，您可以看到主节点的公共 DNS，并可以使用 EC2 密钥对访问主节点的终端。

警告！如果您不想继续使用 EMR 服务，请终止正在运行的实例。否则，您将按实例每小时收费。

参见

  • Google 也提供了自己的云解决方案，即 Google 计算引擎。如果您想了解更多信息，请参阅 cloud.google.com/compute/。

附录 A. R 和机器学习资源

以下表格列出了 R 和机器学习的所有资源：

R 简介
标题链接作者
《R 动作》www.amazon.com/R-Action-Robert-Kabacoff/dp/1935182390Robert Kabacoff
《R 编程艺术：统计软件设计之旅》www.amazon.com/The-Art-Programming-Statistical-Software/dp/1593273843Norman Matloff
《R 简介》cran.r-project.org/doc/manuals/R-intro.pdfW. N. Venables, D. M. Smith, 和 R 核心团队
Quick-Rwww.statmethods.net/罗伯特·I·卡巴科夫，博士
在线课程
标题链接讲师
使用 R 进行数据分析（含 R）www.coursera.org/course/compdata罗杰·D·彭，约翰霍普金斯大学
数据分析www.coursera.org/course/dataanalysis杰夫·利克，约翰霍普金斯大学
数据分析与统计推断www.coursera.org/course/statistics矿·切廷卡亚-朗德尔，杜克大学
机器学习
标题链接作者
为黑客的机器学习www.amazon.com/dp/1449303714?tag=inspiredalgor-20德鲁·康威和约翰·迈尔斯·怀特
使用 R 进行机器学习www.packtpub.com/machine-learning-with-r/book布雷特·兰茨
在线博客
标题链接
R 博客www.r-bloggers.com/
R 杂志journal.r-project.org/
CRAN 任务视图
标题链接
CRAN 任务视图：机器学习和统计学习cran.r-project.org/web/views/MachineLearning.html

附录 B. 数据集 – 泰坦尼克号乘客的生存情况

在探索过程之前，我们想介绍这里采用的例子。它是关于 RMS 泰坦尼克号上乘客的人口统计信息，由 Kaggle 提供（www.kaggle.com/，一个数据预测竞赛的平台）。我们正在考察的结果是乘客是否能在船难中幸存。

应用这个数据集有两个原因：

  • RMS 泰坦尼克号被认为是历史上最臭名昭著的船难，死亡人数高达 1,502 人，占 2,224 名乘客和船员的近三分之二。然而，在船沉没后，乘客的生存机会并不仅仅取决于运气；实际上，舱位等级、性别、年龄和其他因素也可能影响了他们的生存机会。

  • 该数据集相对简单；你不需要花大部分时间在数据处理上（除非处理一些缺失值），但你可以专注于探索性分析的应用。

下面的图表是目标数据集的变量描述：

数据集 – 泰坦尼克号乘客的生存情况

从变量的描述来看，人们可能会有一些疑问，例如，“这个数据集中是否有缺失值？”、“泰坦尼克号乘客的平均年龄是多少？”、“有多少比例的乘客在灾难中幸存？”、“船上大多数乘客属于哪个社会阶层？”在这里提出的所有这些问题都将在 第二章 《使用 RMS 泰坦尼克号进行数据探索》 中得到解答。

除了与描述性统计相关的问题之外，第二章 《使用 RMS 泰坦尼克号进行数据探索》 的最终目标是生成一个模型，用以预测给定输入参数下的生存概率。除此之外，我们还将评估所生成模型的性能，以确定该模型是否适合该问题。

avatar
文章
阅读
粉丝