大数据分析 - 快速指南

大数据分析 - 概述

在过去十年中，需要处理的数据量激增到难以想象的程度，与此同时，数据存储的成本也系统性地降低了。私营公司和研究机构捕获了关于用户互动、业务、社交媒体以及来自移动电话和汽车等设备传感器的海量数据。这个时代的挑战在于理解这片数据海洋。这就是大数据分析发挥作用的地方。

大数据分析主要涉及从不同来源收集数据，对其进行整理，使其能够被分析人员使用，最终交付对组织业务有用的数据产品。

将从不同来源获取的大量非结构化原始数据转换为对组织有用的数据产品，构成了大数据分析的核心。

大数据分析 - 数据生命周期

传统数据挖掘生命周期

为了提供一个框架来组织组织所需的工作并从大数据中提取清晰的见解，将其视为一个具有不同阶段的周期非常有用。它绝不是线性的，这意味着所有阶段彼此相关。这个周期与CRISP方法论中描述的更传统的数据挖掘周期具有表面上的相似性。

CRISP-DM方法论

CRISP-DM方法论代表跨行业标准数据挖掘流程，是一个描述数据挖掘专家用来解决传统BI数据挖掘问题中常用方法的循环。它仍在传统BI数据挖掘团队中使用。

请查看下图。它显示了CRISP-DM方法论描述的循环的主要阶段以及它们之间是如何相互关联的。

CRISP-DM于1996年构思，次年作为欧盟在ESPRIT资助计划下的一个项目启动。该项目由五家公司牵头：SPSS、Teradata、戴姆勒股份公司、NCR公司和OHRA（一家保险公司）。该项目最终被整合到SPSS中。该方法论极其注重如何详细指定数据挖掘项目。

现在让我们更深入地了解CRISP-DM生命周期中涉及的每个阶段：

业务理解 - 此初始阶段侧重于从业务角度理解项目目标和需求，然后将此知识转化为数据挖掘问题定义。设计初步计划以实现目标。可以使用决策模型，特别是使用决策模型和符号标准构建的模型。
数据理解 - 数据理解阶段从初始数据收集开始，并继续进行活动，以便熟悉数据，识别数据质量问题，发现对数据的初步见解，或检测有趣的子集以形成隐藏信息的假设。
数据准备 - 数据准备阶段涵盖所有活动，以从初始原始数据构建最终数据集（将馈送到建模工具的数据）。数据准备任务很可能多次执行，并且不按任何规定的顺序执行。任务包括表、记录和属性选择，以及为建模工具转换和清理数据。
建模 - 在此阶段，选择和应用各种建模技术，并将其参数校准到最佳值。通常，对于相同的数据挖掘问题类型，有多种技术。一些技术对数据形式有特殊要求。因此，通常需要回到数据准备阶段。
评估 - 在项目的这个阶段，您已经构建了一个（或多个）从数据分析角度来看似乎具有高质量的模型。在继续进行模型的最终部署之前，务必彻底评估模型并审查构建模型所执行的步骤，以确保它能够正确实现业务目标。

一个关键目标是确定是否存在尚未充分考虑的重要业务问题。在本阶段结束时，应就使用数据挖掘结果做出决定。
部署 - 模型的创建通常不是项目的结束。即使模型的目的是增加对数据的了解，获得的知识也需要以对客户有用的方式进行组织和呈现。

根据需求，部署阶段可以像生成报告一样简单，也可以像实现可重复的数据评分（例如，细分分配）或数据挖掘过程一样复杂。

在许多情况下，将是客户而不是数据分析师执行部署步骤。即使分析师部署了模型，客户也需要预先了解需要执行哪些操作才能实际使用创建的模型。

SEMMA方法论

SEMMA是SAS为数据挖掘建模开发的另一种方法论。它代表Sample（抽样）、Explore（探索）、Modify（修改）、Model（建模）和Asses（评估）。以下是其各个阶段的简要说明：

抽样 - 该过程从数据抽样开始，例如，选择用于建模的数据集。数据集应足够大，包含足够的信息来检索，但又足够小以高效使用。此阶段还处理数据分区。
探索 - 此阶段涵盖通过发现变量之间预期和意外的关系以及异常值（借助数据可视化）来理解数据。
修改 - 修改阶段包含用于选择、创建和转换变量以准备数据建模的方法。
建模 - 在建模阶段，重点是将各种建模（数据挖掘）技术应用于准备好的变量，以便创建可能提供所需结果的模型。
评估 - 对建模结果的评估显示了创建模型的可靠性和实用性。

CRISM-DM和SEMMA的主要区别在于，SEMMA侧重于建模方面，而CRISP-DM更重视建模之前的周期阶段，例如理解要解决的业务问题，理解和预处理将用作输入的数据（例如，机器学习算法）。

大数据生命周期

在当今大数据环境中，以前的方法要么不完整，要么次优。例如，SEMMA方法论完全忽略了不同数据源的数据收集和预处理。这些阶段通常构成成功大数据项目的大部分工作。

大数据分析周期可以通过以下阶段来描述：

业务问题定义
研究
人力资源评估
数据采集
数据清洗
数据存储
探索性数据分析
数据准备建模和评估
建模
实施

在本节中，我们将对大数据生命周期的这些阶段进行一些阐述。

业务问题定义

这是传统BI和大数据分析生命周期中的一个共同点。通常，定义问题并正确评估它可能为组织带来的潜在收益，是大数据项目的一个非平凡阶段。提及这一点似乎显而易见，但必须评估项目的预期收益和成本。

研究

分析其他公司在相同情况下所做的事情。这包括寻找对您的公司来说合理的解决方案，即使这意味着将其他解决方案适应您的公司拥有的资源和要求。在此阶段，应定义未来阶段的方法论。

人力资源评估

一旦定义了问题，接下来就合理地分析当前员工是否能够成功完成项目。传统的BI团队可能无法为所有阶段提供最佳解决方案，因此在启动项目之前，应该考虑是否有必要外包项目的一部分或雇佣更多人员。

数据采集

本节是大数据生命周期中的关键；它定义了哪些类型的配置文件将需要交付最终数据产品。数据收集是该过程中的一个非平凡步骤；它通常涉及从不同来源收集非结构化数据。例如，它可能涉及编写爬虫程序从网站检索评论。这涉及处理文本，可能使用不同的语言，通常需要大量时间才能完成。

数据清洗

一旦数据被检索（例如，从网络），就需要将其存储在易于使用的格式中。继续以评论为例，假设数据是从不同的站点检索的，每个站点的数据显示方式都不同。

假设一个数据源以星级评定的形式提供评论，因此可以将其读取为响应变量y ∈ {1, 2, 3, 4, 5}的映射。另一个数据源使用两个箭头系统提供评论，一个用于向上投票，另一个用于向下投票。这意味着响应变量的形式为y ∈ {正面，负面}。

为了组合这两个数据源，必须做出决定以使这两个响应表示等效。这可能涉及将第一个数据源的响应表示转换为第二种形式，将一颗星视为负面，五颗星视为正面。此过程通常需要大量时间才能高质量地交付。

数据存储

数据处理完成后，有时需要将其存储在数据库中。大数据技术为此提供了多种选择。最常见的选择是使用Hadoop文件系统进行存储，它为用户提供了一个有限版本的SQL，称为Hive查询语言。从用户的角度来看，这允许大多数分析任务以与传统BI数据仓库中类似的方式完成。其他可考虑的存储选项包括MongoDB、Redis和SPARK。

这个阶段与人力资源的知识有关，具体来说是他们实施不同架构的能力。传统数据仓库的修改版本仍在大型应用中使用。例如，Teradata和IBM提供可以处理TB级数据的SQL数据库；PostgreSQL和MySQL等开源解决方案仍在大型应用中使用。

尽管不同存储在后台的工作方式有所不同，但从客户端来看，大多数解决方案都提供SQL API。因此，深入了解SQL仍然是大数据分析的关键技能。

这个阶段事先似乎是最重要的主题，但实际上并非如此，它甚至不是必需的阶段。可以实现一个使用实时数据的大数据解决方案，在这种情况下，我们只需要收集数据来开发模型，然后实时实施它。因此，根本不需要正式存储数据。

探索性数据分析

一旦数据被清洗并以可以从中检索洞察信息的方式存储，数据探索阶段就是强制性的。此阶段的目标是理解数据，这通常是通过统计技术和数据绘图完成的。这是一个评估问题定义是否有意义或是否可行的良好阶段。

数据准备建模和评估

此阶段涉及重塑先前检索到的清洗数据，并使用统计预处理进行缺失值插补、异常值检测、归一化、特征提取和特征选择。

建模

之前的阶段应该已经产生了多个用于训练和测试的数据集，例如，一个预测模型。此阶段涉及尝试不同的模型，并期待解决手头的业务问题。实际上，通常希望模型能够对业务提供一些见解。最后，选择最佳模型或模型组合，评估其在留出数据集上的性能。

实施

在此阶段，开发的数据产品将被部署到公司的数管道中。这涉及在数据产品运行时设置验证方案，以跟踪其性能。例如，在实施预测模型的情况下，此阶段将涉及将模型应用于新数据，并在响应可用后评估模型。

大数据分析 - 方法论

在方法论方面，大数据分析与传统的实验设计统计方法有很大不同。分析始于数据。我们通常以某种方式对数据建模以解释响应。这种方法的目标是预测响应行为或理解输入变量如何与响应相关。在统计实验设计中，通常会开发一个实验，并由此检索数据。这允许以统计模型可以使用的方式生成数据，其中某些假设成立，例如独立性、正态性和随机性。

在大数据分析中，我们得到了数据。我们无法设计一个满足我们最喜欢的统计模型的实验。在大规模应用分析中，大量工作（通常占80%的努力）仅用于数据清洗，以便机器学习模型可以使用它。

在实际的大规模应用中，我们没有一个独特的方法论可以遵循。通常，一旦定义了业务问题，就需要一个研究阶段来设计使用方法论。但是，一些通用指南是相关的，并且几乎适用于所有问题。

大数据分析中最重要的任务之一是**统计建模**，这意味着监督和无监督分类或回归问题。一旦数据被清洗和预处理，就可以用于建模，应注意使用合理的损失度量评估不同的模型，然后一旦模型被实施，就应该进一步评估和报告结果。预测建模中一个常见的陷阱是仅仅实施模型而从不衡量其性能。

大数据分析 - 核心交付成果

如大数据生命周期中所述，开发大数据产品产生的数据产品在大多数情况下是以下一些——

**机器学习实现**——这可能是一个分类算法、一个回归模型或一个细分模型。
**推荐系统**——目标是开发一个基于用户行为推荐选择的系统。**Netflix**是这种数据产品的典型示例，它根据用户的评分推荐其他电影。
**仪表盘**——企业通常需要工具来可视化聚合数据。仪表盘是一种图形机制，使这些数据可访问。
**临时分析**——通常业务领域有一些问题、假设或谬论可以通过对数据进行临时分析来解答。

大数据分析 - 主要利益相关者

在大公司中，为了成功地开发大数据项目，需要管理层支持该项目。这通常涉及找到一种方法来展示该项目的业务优势。我们没有一个独特的解决方案来寻找项目的赞助商，但下面给出了一些指导方针——

检查谁以及在哪里是与您感兴趣的项目类似的其他项目的赞助商。
在关键管理职位上有私人联系会有帮助，因此如果项目有前景，可以启动任何联系。
谁会从您的项目中受益？项目启动后，您的客户是谁？
制定一个简单、清晰且令人兴奋的提案，并与您组织中的关键参与者分享。

寻找项目赞助商的最佳方法是理解问题以及实施后将产生的数据产品是什么。这种理解将有助于说服管理层大数据项目的重要性。

大数据分析 - 数据分析师

数据分析师拥有面向报告的个人资料，具有使用SQL从传统数据仓库中提取和分析数据的经验。他们的任务通常是在数据存储方面或报告一般业务结果方面。

许多组织都在努力寻找市场上合格的数据科学家。然而，选择未来的数据分析师并教授他们成为数据科学家的相关技能是一个好主意。这绝非易事，通常需要这个人获得定量领域的硕士学位，但这绝对是一个可行的选择。合格的数据分析师必须具备的基本技能如下：

商业理解
SQL编程
报表设计和实施
仪表盘开发

大数据分析 - 数据科学家

数据科学家的角色通常与诸如预测建模、开发细分算法、推荐系统、A/B 测试框架以及经常处理原始非结构化数据等任务相关联。

他们工作的性质需要对数学、应用统计和编程有深入的理解。数据分析师和数据科学家之间有一些共同的技能，例如查询数据库的能力。两者都分析数据，但数据科学家的决策可能会对组织产生更大的影响。

以下是数据科学家通常需要具备的一些技能：

使用统计软件包进行编程，例如：R、Python、SAS、SPSS或Julia
能够从不同来源清洗、提取和探索数据
统计模型的研究、设计和实施
深厚的统计、数学和计算机科学知识

在大数据分析中，人们通常会混淆数据科学家和数据架构师的角色。实际上，区别很简单。数据架构师定义数据将存储在其中的工具和架构，而数据科学家则使用此架构。当然，如果需要进行临时项目，数据科学家应该能够设置新的工具，但基础设施的定义和设计不应属于他的任务。

大数据分析 - 问题定义

在本教程中，我们将开发一个项目。本教程中的每个后续章节都处理迷你项目部分中较大项目的一部分。这被认为是一个应用教程部分，将提供对现实世界问题的了解。在这种情况下，我们将从项目的问题定义开始。

项目描述

本项目的目标是开发一个机器学习模型，使用求职者的简历文本作为输入来预测他们的每小时薪水。

使用上面定义的框架，很容易定义这个问题。我们可以将 *X = {x₁, x₂, …, x_n}* 定义为用户的简历，其中每个特征都可以以最简单的方式表示该单词出现的次数。然后响应是实值，我们试图预测个人的每小时美元薪水。

这两个考虑足以得出结论：提出的问题可以用监督回归算法来解决。

问题定义

**问题定义**可能是大数据分析管道中最复杂和最被忽视的阶段之一。为了定义数据产品将解决的问题，经验是强制性的。大多数 aspiring 数据科学家在这个阶段几乎没有经验。

大多数大数据问题可以分为以下几种类型：

监督分类
监督回归
无监督学习
学习排序

现在让我们进一步了解这四个概念。

监督分类

给定一个特征矩阵 *X = {x₁, x₂, ..., x_n}*，我们开发一个模型 M 来预测定义为 *y = {c₁, c₂, ..., c_n}* 的不同类别。例如：给定保险公司客户的交易数据，可以开发一个模型来预测客户是否会流失。后者是一个二元分类问题，有两个类别或目标变量：流失和未流失。

其他问题涉及预测多个类别，我们可能对进行数字识别感兴趣，因此响应向量将定义为：*y = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}*，最先进的模型将是卷积神经网络，特征矩阵将定义为图像的像素。

监督回归

在这种情况下，问题定义与前面的示例非常相似；区别在于响应。在回归问题中，响应 y ∈ ℜ，这意味着响应是实值的。例如，我们可以开发一个模型，根据求职者的简历内容预测他们的每小时薪水。

无监督学习

管理层常常渴望获得新的见解。细分模型可以提供这种见解，以便营销部门为不同的细分市场开发产品。开发细分模型的一个好方法，与其考虑算法，不如选择与所需细分相关的特征。

例如，在电信公司中，按客户的手机使用情况细分客户是很有意义的。这将涉及忽略与细分目标无关的特征，而只包含相关的特征。在这种情况下，这将选择诸如每月使用的短信数量、呼入和呼出分钟数等特征。

学习排序

这个问题可以被认为是一个回归问题，但它具有特殊的特点，值得单独处理。这个问题涉及到给定一个文档集合，我们寻求找到给定查询的最相关的排序。为了开发一个监督学习算法，需要标记给定查询的排序的相关性。

值得注意的是，为了开发监督学习算法，需要标记训练数据。这意味着，为了训练一个模型，例如，从图像中识别数字，我们需要手工标记大量的示例。有一些网络服务可以加快这个过程，并且通常用于此任务，例如亚马逊Mechanical Turk。事实证明，当提供更多数据时，学习算法会提高其性能，因此在监督学习中，标记相当数量的示例实际上是强制性的。

大数据分析 - 数据收集

数据收集在大数据周期中扮演着最重要的角色。互联网为各种主题提供了几乎无限的数据来源。该领域的重要性取决于业务类型，但传统行业可以获取各种外部数据来源，并将其与交易数据相结合。

例如，假设我们想构建一个推荐餐厅的系统。第一步是收集数据，在这种情况下，是从不同的网站收集餐厅评论并将它们存储在数据库中。由于我们感兴趣的是原始文本，并将使用它进行分析，因此数据存储位置与开发模型无关。这听起来可能与大数据主要技术相矛盾，但是为了实现大数据应用程序，我们只需要使其实时运行。

Twitter小型项目

一旦问题被定义，接下来的阶段就是收集数据。以下小型项目的想法是致力于从网络上收集数据并将其结构化，以便用于机器学习模型。我们将使用R编程语言从Twitter REST API收集一些推文。

首先创建一个Twitter帐户，然后按照**twitteR**软件包vignette中的说明创建一个Twitter开发者帐户。这是这些说明的摘要：

访问https://twitter.com/apps/new并登录。
填写基本信息后，转到“设置”选项卡，然后选择“读取、写入和访问直接消息”。
确保在执行此操作后单击保存按钮。
在“详细信息”选项卡中，记下您的消费者密钥和消费者密钥。
在您的R会话中，您将使用API密钥和API密钥值。
最后运行以下脚本。这将从其在github上的存储库安装**twitteR**软件包。

install.packages(c("devtools", "rjson", "bit64", "httr"))  

# Make sure to restart your R session at this point 
library(devtools) 
install_github("geoffjentry/twitteR")

我们感兴趣的是获取包含字符串“big mac”的数据，并找出哪些主题在此方面脱颖而出。为此，第一步是从Twitter收集数据。以下是我们收集Twitter所需数据的R脚本。此代码也可在bda/part1/collect_data/collect_data_twitter.R文件中找到。

rm(list = ls(all = TRUE)); gc() # Clears the global environment
library(twitteR)
Sys.setlocale(category = "LC_ALL", locale = "C")

### Replace the xxx’s with the values you got from the previous instructions

# consumer_key = "xxxxxxxxxxxxxxxxxxxx"
# consumer_secret = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"

# access_token = "xxxxxxxxxx-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"

# access_token_secret= "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"

# Connect to twitter rest API
setup_twitter_oauth(consumer_key, consumer_secret, access_token, access_token_secret)

# Get tweets related to big mac
tweets <- searchTwitter(’big mac’, n = 200, lang = ’en’)
df <- twListToDF(tweets)

# Take a look at the data
head(df)

# Check which device is most used
sources <- sapply(tweets, function(x) x$getStatusSource())
sources <- gsub("</a>", "", sources)
sources <- strsplit(sources, ">")
sources <- sapply(sources, function(x) ifelse(length(x) > 1, x[2], x[1]))
source_table = table(sources)
source_table = source_table[source_table > 1]
freq = source_table[order(source_table, decreasing = T)]
as.data.frame(freq)

#                       Frequency
# Twitter for iPhone       71
# Twitter for Android      29
# Twitter Web Client       25
# recognia                 20

大数据分析 - 数据清洗

一旦数据被收集，我们通常会有来自不同来源、具有不同特征的数据。最直接的步骤是使这些数据源同质化，并继续开发我们的数据产品。然而，这取决于数据的类型。我们应该问问自己，使数据同质化是否实际。

也许数据来源完全不同，如果对来源进行同质化，信息丢失将会很大。在这种情况下，我们可以考虑替代方案。一个数据源能否帮助我构建回归模型，而另一个数据源构建分类模型？是否可以利用异质性而不是仅仅丢失信息？做出这些决定是使分析变得有趣和具有挑战性的原因。

对于评论来说，每个数据源都可能有一种语言。同样，我们有两个选择：

**同质化** - 这涉及将不同的语言翻译成我们拥有更多数据的语言。翻译服务的质量是可以接受的，但是如果我们想使用API翻译大量数据，成本将非常高。有一些软件工具可以完成这项任务，但这也会很昂贵。
**异质化** - 是否可以为每种语言开发一种解决方案？由于很容易检测语料库的语言，我们可以为每种语言开发一个推荐器。这将涉及更多工作，需要根据可用语言的数量调整每个推荐器，但如果我们只有几种语言可用，这绝对是一个可行的选择。

Twitter小型项目

在本例中，我们需要首先清理非结构化数据，然后将其转换为数据矩阵，以便对其应用主题建模。一般来说，从Twitter获取数据时，有一些字符我们不感兴趣，至少在数据清洗过程的第一阶段是这样。

例如，在获取推文后，我们得到这些奇怪的字符：“<ed><U+00A0><U+00BD><ed><U+00B8><U+008B>”。这些可能是表情符号，因此为了清理数据，我们将使用以下脚本将其删除。此代码也可在bda/part1/collect_data/cleaning_data.R文件中找到。

rm(list = ls(all = TRUE)); gc() # Clears the global environment
source('collect_data_twitter.R')
# Some tweets
head(df$text)

[1] "I’m not a big fan of turkey but baked Mac &
cheese <ed><U+00A0><U+00BD><ed><U+00B8><U+008B>"
[2] "@Jayoh30 Like no special sauce on a big mac. HOW"
### We are interested in the text - Let’s clean it!

# We first convert the encoding of the text from latin1 to ASCII
df$text <- sapply(df$text,function(row) iconv(row, "latin1", "ASCII", sub = ""))

# Create a function to clean tweets
clean.text <- function(tx) {
  tx <- gsub("htt.{1,20}", " ", tx, ignore.case = TRUE)
  tx = gsub("[^#[:^punct:]]|@|RT", " ", tx, perl = TRUE, ignore.case = TRUE)
  tx = gsub("[[:digit:]]", " ", tx, ignore.case = TRUE)
  tx = gsub(" {1,}", " ", tx, ignore.case = TRUE)
  tx = gsub("^\\s+|\\s+$", " ", tx, ignore.case = TRUE)
  return(tx)
}  

clean_tweets <- lapply(df$text, clean.text)

# Cleaned tweets
head(clean_tweets)
[1] " WeNeedFeminlsm MAC s new make up line features men woc and big girls "
[1] " TravelsPhoto What Happens To Your Body One Hour After A Big Mac "

数据清洗小型项目的最后一步是拥有经过清理的文本，我们可以将其转换为矩阵并应用算法。从存储在**clean_tweets**向量中的文本中，我们可以轻松地将其转换为词袋矩阵并应用无监督学习算法。

大数据分析 - 数据汇总

报告在大数据分析中非常重要。每个组织都必须定期提供信息来支持其决策过程。这项任务通常由具有SQL和ETL（提取、转换和加载）经验的数据分析师处理。

负责这项任务的团队有责任将大数据分析部门产生的信息传播到组织的不同领域。

以下示例演示了数据汇总的含义。导航到文件夹**bda/part1/summarize_data**，并在文件夹内双击打开**summarize_data.Rproj**文件。然后，打开**summarize_data.R**脚本并查看代码，并按照提供的说明进行操作。

# Install the following packages by running the following code in R. 
pkgs = c('data.table', 'ggplot2', 'nycflights13', 'reshape2') 
install.packages(pkgs)

**ggplot2**软件包非常适合数据可视化。**data.table**软件包是在**R**中进行快速且内存高效的汇总的一个不错的选择。最近的基准测试表明，它甚至比**pandas**（用于执行类似任务的Python库）更快。

使用以下代码查看数据。此代码也可在**bda/part1/summarize_data/summarize_data.Rproj**文件中找到。

library(nycflights13) 
library(ggplot2) 
library(data.table) 
library(reshape2)  

# Convert the flights data.frame to a data.table object and call it DT 
DT <- as.data.table(flights)  

# The data has 336776 rows and 16 columns 
dim(DT)  

# Take a look at the first rows 
head(DT) 

#   year    month  day   dep_time  dep_delay  arr_time  arr_delay  carrier 
# 1: 2013     1     1      517       2         830         11       UA 
# 2: 2013     1     1      533       4         850         20       UA 
# 3: 2013     1     1      542       2         923         33       AA 
# 4: 2013     1     1      544      -1        1004        -18       B6 
# 5: 2013     1     1      554      -6         812        -25       DL 
# 6: 2013     1     1      554      -4         740         12       UA  

#     tailnum  flight  origin   dest    air_time   distance    hour   minute 
# 1:  N14228   1545     EWR      IAH      227        1400       5       17 
# 2:  N24211   1714     LGA      IAH      227        1416       5       33 
# 3:  N619AA   1141     JFK      MIA      160        1089       5       42 
# 4:  N804JB    725     JFK      BQN      183        1576       5       44 
# 5:  N668DN    461     LGA      ATL      116        762        5       54 
# 6:  N39463   1696     EWR      ORD      150        719        5       54

以下代码包含数据汇总示例。

### Data Summarization
# Compute the mean arrival delay  
DT[, list(mean_arrival_delay = mean(arr_delay, na.rm = TRUE))] 
#        mean_arrival_delay 
# 1:           6.895377  
# Now, we compute the same value but for each carrier 
mean1 = DT[, list(mean_arrival_delay = mean(arr_delay, na.rm = TRUE)), 
   by = carrier] 
print(mean1) 
#      carrier    mean_arrival_delay 
# 1:      UA          3.5580111 
# 2:      AA          0.3642909 
# 3:      B6          9.4579733 
# 4:      DL          1.6443409 
# 5:      EV         15.7964311 
# 6:      MQ         10.7747334 
# 7:      US          2.1295951 
# 8:      WN          9.6491199 
# 9:      VX          1.7644644 
# 10:     FL         20.1159055 
# 11:     AS         -9.9308886 
# 12:     9E          7.3796692
# 13:     F9         21.9207048 
# 14:     HA         -6.9152047 
# 15:     YV         15.5569853 
# 16:     OO         11.9310345

# Now let’s compute to means in the same line of code 
mean2 = DT[, list(mean_departure_delay = mean(dep_delay, na.rm = TRUE), 
   mean_arrival_delay = mean(arr_delay, na.rm = TRUE)), 
   by = carrier] 
print(mean2) 

#       carrier    mean_departure_delay   mean_arrival_delay 
# 1:      UA            12.106073          3.5580111 
# 2:      AA             8.586016          0.3642909 
# 3:      B6            13.022522          9.4579733 
# 4:      DL             9.264505          1.6443409 
# 5:      EV            19.955390         15.7964311 
# 6:      MQ            10.552041         10.7747334 
# 7:      US             3.782418          2.1295951 
# 8:      WN            17.711744          9.6491199 
# 9:      VX            12.869421          1.7644644 
# 10:     FL            18.726075         20.1159055 
# 11:     AS             5.804775         -9.9308886 
# 12:     9E            16.725769          7.3796692 
# 13:     F9            20.215543         21.9207048 
# 14:     HA             4.900585         -6.9152047 
# 15:     YV            18.996330         15.5569853 
# 16:     OO            12.586207         11.9310345

### Create a new variable called gain 
# this is the difference between arrival delay and departure delay 
DT[, gain:= arr_delay - dep_delay]  

# Compute the median gain per carrier 
median_gain = DT[, median(gain, na.rm = TRUE), by = carrier] 
print(median_gain)

大数据分析 - 数据探索

**探索性数据分析**是由John Tuckey（1977）提出的一个概念，它代表了统计学的一种新视角。Tuckey的想法是，在传统的统计学中，数据没有以图形方式进行探索，它只是用于检验假设。第一次尝试开发工具是在斯坦福大学进行的，该项目名为prim9。该工具能够将数据可视化为九个维度，因此能够提供数据的多元视角。

近年来，探索性数据分析是必不可少的，并且已被纳入大数据分析生命周期。能够发现见解并能够在组织中有效地传达见解的能力，是强大的EDA能力的推动。

基于Tuckey的想法，贝尔实验室开发了**S编程语言**，以便提供一个交互式界面来进行统计分析。S的想法是提供具有易于使用的语言的广泛图形功能。在当今的大数据背景下，基于**S**编程语言的**R**是最流行的分析软件。

以下程序演示了探索性数据分析的使用。

以下是探索性数据分析的示例。此代码也可在**part1/eda/exploratory_data_analysis.R**文件中找到。

library(nycflights13) 
library(ggplot2) 
library(data.table) 
library(reshape2)  

# Using the code from the previous section 
# This computes the mean arrival and departure delays by carrier. 
DT <- as.data.table(flights) 
mean2 = DT[, list(mean_departure_delay = mean(dep_delay, na.rm = TRUE), 
   mean_arrival_delay = mean(arr_delay, na.rm = TRUE)), 
   by = carrier]  

# In order to plot data in R usign ggplot, it is normally needed to reshape the data 
# We want to have the data in long format for plotting with ggplot 
dt = melt(mean2, id.vars = ’carrier’)  

# Take a look at the first rows 
print(head(dt))  

# Take a look at the help for ?geom_point and geom_line to find similar examples 
# Here we take the carrier code as the x axis 
# the value from the dt data.table goes in the y axis 

# The variable column represents the color 
p = ggplot(dt, aes(x = carrier, y = value, color = variable, group = variable)) +
   geom_point() + # Plots points 
   geom_line() + # Plots lines 
   theme_bw() + # Uses a white background 
   labs(list(title = 'Mean arrival and departure delay by carrier', 
      x = 'Carrier', y = 'Mean delay')) 
print(p)  

# Save the plot to disk 
ggsave('mean_delay_by_carrier.png', p,  
   width = 10.4, height = 5.07)

该代码应该生成如下所示的图像：

大数据分析 - 数据可视化

为了理解数据，将其可视化通常很有用。通常在大数据应用程序中，兴趣在于发现见解，而不仅仅是制作漂亮的图表。以下是使用图表来理解数据的一些不同方法的示例。

要开始分析航班数据，我们可以首先检查数值变量之间是否存在相关性。此代码也可在**bda/part1/data_visualization/data_visualization.R**文件中找到。

# Install the package corrplot by running
install.packages('corrplot')  

# then load the library 
library(corrplot)  

# Load the following libraries  
library(nycflights13) 
library(ggplot2) 
library(data.table) 
library(reshape2)  

# We will continue working with the flights data 
DT <- as.data.table(flights)  
head(DT) # take a look  

# We select the numeric variables after inspecting the first rows. 
numeric_variables = c('dep_time', 'dep_delay',  
   'arr_time', 'arr_delay', 'air_time', 'distance')

# Select numeric variables from the DT data.table 
dt_num = DT[, numeric_variables, with = FALSE]  

# Compute the correlation matrix of dt_num 
cor_mat = cor(dt_num, use = "complete.obs")  

print(cor_mat) 
### Here is the correlation matrix 
#              dep_time   dep_delay   arr_time   arr_delay    air_time    distance 
# dep_time   1.00000000  0.25961272 0.66250900  0.23230573 -0.01461948 -0.01413373 
# dep_delay  0.25961272  1.00000000 0.02942101  0.91480276 -0.02240508 -0.02168090 
# arr_time   0.66250900  0.02942101 1.00000000  0.02448214  0.05429603  0.04718917 
# arr_delay  0.23230573  0.91480276 0.02448214  1.00000000 -0.03529709 -0.06186776 
# air_time  -0.01461948 -0.02240508 0.05429603 -0.03529709  1.00000000  0.99064965 
# distance  -0.01413373 -0.02168090 0.04718917 -0.06186776  0.99064965  1.00000000  

# We can display it visually to get a better understanding of the data 
corrplot.mixed(cor_mat, lower = "circle", upper = "ellipse")  

# save it to disk 
png('corrplot.png') 
print(corrplot.mixed(cor_mat, lower = "circle", upper = "ellipse")) 
dev.off()

此代码生成以下相关矩阵可视化：

我们可以从图中看到，数据集中的某些变量之间存在很强的相关性。例如，到达延误和出发延误似乎高度相关。我们可以看到这一点，因为椭圆显示这两个变量之间几乎呈线性关系，但是，从这个结果中找到因果关系并不简单。

我们不能说，因为两个变量是相关的，所以一个变量会影响另一个变量。我们还在图中发现飞行时间和距离之间存在很强的相关性，这在预期中是相当合理的，因为距离越远，飞行时间应该越长。

我们还可以对数据进行单变量分析。可视化分布的一种简单有效的方法是**箱线图**。以下代码演示了如何使用ggplot2库生成箱线图和格子图。此代码也可在**bda/part1/data_visualization/boxplots.R**文件中找到。

source('data_visualization.R') 
### Analyzing Distributions using box-plots  
# The following shows the distance as a function of the carrier 

p = ggplot(DT, aes(x = carrier, y = distance, fill = carrier)) + # Define the carrier 
   in the x axis and distance in the y axis 
   geom_box-plot() + # Use the box-plot geom 
   theme_bw() + # Leave a white background - More in line with tufte's 
      principles than the default 
   guides(fill = FALSE) + # Remove legend 
   labs(list(title = 'Distance as a function of carrier', # Add labels 
      x = 'Carrier', y = 'Distance')) 
p   
# Save to disk 
png(‘boxplot_carrier.png’) 
print(p) 
dev.off()   

# Let's add now another variable, the month of each flight 
# We will be using facet_wrap for this 
p = ggplot(DT, aes(carrier, distance, fill = carrier)) + 
   geom_box-plot() + 
   theme_bw() + 
   guides(fill = FALSE) +  
   facet_wrap(~month) + # This creates the trellis plot with the by month variable
   labs(list(title = 'Distance as a function of carrier by month', 
      x = 'Carrier', y = 'Distance')) 
p   
# The plot shows there aren't clear differences between distance in different months  

# Save to disk 
png('boxplot_carrier_by_month.png') 
print(p) 
dev.off()

大数据分析 - R语言入门

本节旨在向用户介绍R编程语言。R可以从cran网站下载。对于Windows用户，安装rtools和rstudio IDE很有用。

**R**背后的总体概念是充当用编译语言（如C、C++和Fortran）开发的其他软件的接口，并为用户提供一个交互式工具来分析数据。

导航到书籍压缩文件bda/part2/R_introduction的文件夹，并打开R_introduction.Rproj文件。这将打开一个RStudio会话。然后打开01_vectors.R文件。逐行运行脚本并遵循代码中的注释。另一个有用的学习方法是直接键入代码，这将帮助你熟悉R语法。在R中，注释用#符号编写。

为了显示本书中运行R代码的结果，在代码评估后，R返回的结果将被注释。这样，你就可以复制粘贴书中的代码，并在R中直接尝试其中的部分代码。

# Create a vector of numbers 
numbers = c(1, 2, 3, 4, 5) 
print(numbers) 

# [1] 1 2 3 4 5  
# Create a vector of letters 
ltrs = c('a', 'b', 'c', 'd', 'e') 
# [1] "a" "b" "c" "d" "e"  

# Concatenate both  
mixed_vec = c(numbers, ltrs) 
print(mixed_vec) 
# [1] "1" "2" "3" "4" "5" "a" "b" "c" "d" "e"

让我们分析一下前面代码中发生的情况。我们可以看到，可以用数字和字母创建向量。我们不需要事先告诉R我们想要哪种数据类型。最后，我们能够创建一个包含数字和字母的向量。向量mixed_vec已将数字强制转换为字符，我们可以通过查看值如何在引号内打印来看到这一点。

下面的代码显示了由class函数返回的不同向量的类型。通常使用class函数来“询问”对象，询问它的类是什么。

### Evaluate the data types using class

### One dimensional objects 
# Integer vector 
num = 1:10 
class(num) 
# [1] "integer"  

# Numeric vector, it has a float, 10.5 
num = c(1:10, 10.5) 
class(num) 
# [1] "numeric"  

# Character vector 
ltrs = letters[1:10] 
class(ltrs) 
# [1] "character"  

# Factor vector 
fac = as.factor(ltrs) 
class(fac) 
# [1] "factor"

R也支持二维对象。在下面的代码中，有一些在R中最常用的两种数据结构的示例：矩阵和data.frame。

# Matrix
M = matrix(1:12, ncol = 4) 
#      [,1] [,2] [,3] [,4] 
# [1,]    1    4    7   10 
# [2,]    2    5    8   11 
# [3,]    3    6    9   12 
lM = matrix(letters[1:12], ncol = 4) 
#     [,1] [,2] [,3] [,4] 
# [1,] "a"  "d"  "g"  "j"  
# [2,] "b"  "e"  "h"  "k"  
# [3,] "c"  "f"  "i"  "l"   

# Coerces the numbers to character 
# cbind concatenates two matrices (or vectors) in one matrix 
cbind(M, lM) 
#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] 
# [1,] "1"  "4"  "7"  "10" "a"  "d"  "g"  "j"  
# [2,] "2"  "5"  "8"  "11" "b"  "e"  "h"  "k"  
# [3,] "3"  "6"  "9"  "12" "c"  "f"  "i"  "l"   

class(M) 
# [1] "matrix" 
class(lM) 
# [1] "matrix"  

# data.frame 
# One of the main objects of R, handles different data types in the same object.  
# It is possible to have numeric, character and factor vectors in the same data.frame  

df = data.frame(n = 1:5, l = letters[1:5]) 
df 
#   n l 
# 1 1 a 
# 2 2 b 
# 3 3 c 
# 4 4 d 
# 5 5 e

如前面的示例所示，可以在同一个对象中使用不同的数据类型。通常情况下，数据就是这样在数据库、API中呈现的，部分数据是文本或字符向量，其他是数值型。分析师的工作是确定要分配哪种统计数据类型，然后为其使用正确的R数据类型。在统计学中，我们通常认为变量具有以下类型：

数值型
名义型或分类型
有序型

在R中，向量可以具有以下类别：

数值型 - 整数型
因子型
有序因子型

R为每种统计类型的变量提供了一种数据类型。然而，有序因子很少使用，但可以通过函数factor或ordered创建。

下一节讨论索引的概念。这是一个相当常见的操作，它处理选择对象的部分并对其进行转换的问题。

# Let's create a data.frame
df = data.frame(numbers = 1:26, letters) 
head(df) 
#      numbers  letters 
# 1       1       a 
# 2       2       b 
# 3       3       c 
# 4       4       d 
# 5       5       e 
# 6       6       f 

# str gives the structure of a data.frame, it’s a good summary to inspect an object 
str(df) 
#   'data.frame': 26 obs. of  2 variables: 
#   $ numbers: int  1 2 3 4 5 6 7 8 9 10 ... 
#   $ letters: Factor w/ 26 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ...  

# The latter shows the letters character vector was coerced as a factor. 
# This can be explained by the stringsAsFactors = TRUE argumnet in data.frame 
# read ?data.frame for more information  

class(df) 
# [1] "data.frame"  

### Indexing
# Get the first row 
df[1, ] 
#     numbers  letters 
# 1       1       a  

# Used for programming normally - returns the output as a list 
df[1, , drop = TRUE] 
# $numbers 
# [1] 1 
#  
# $letters 
# [1] a 
# Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z  

# Get several rows of the data.frame 
df[5:7, ] 
#      numbers  letters 
# 5       5       e 
# 6       6       f 
# 7       7       g  

### Add one column that mixes the numeric column with the factor column 
df$mixed = paste(df$numbers, df$letters, sep = ’’)  

str(df) 
# 'data.frame': 26 obs. of  3 variables: 
# $ numbers: int  1 2 3 4 5 6 7 8 9 10 ...
# $ letters: Factor w/ 26 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ... 
# $ mixed  : chr  "1a" "2b" "3c" "4d" ...  

### Get columns 
# Get the first column 
df[, 1]  
# It returns a one dimensional vector with that column  

# Get two columns 
df2 = df[, 1:2] 
head(df2)  

#      numbers  letters 
# 1       1       a 
# 2       2       b 
# 3       3       c 
# 4       4       d 
# 5       5       e 
# 6       6       f  

# Get the first and third columns 
df3 = df[, c(1, 3)] 
df3[1:3, ]  

#      numbers  mixed 
# 1       1     1a
# 2       2     2b 
# 3       3     3c  

### Index columns from their names 
names(df) 
# [1] "numbers" "letters" "mixed"   
# This is the best practice in programming, as many times indeces change, but 
variable names don’t 
# We create a variable with the names we want to subset 
keep_vars = c("numbers", "mixed") 
df4 = df[, keep_vars]  

head(df4) 
#      numbers  mixed 
# 1       1     1a 
# 2       2     2b 
# 3       3     3c 
# 4       4     4d 
# 5       5     5e 
# 6       6     6f  

### subset rows and columns 
# Keep the first five rows 
df5 = df[1:5, keep_vars] 
df5 

#      numbers  mixed 
# 1       1     1a 
# 2       2     2b
# 3       3     3c 
# 4       4     4d 
# 5       5     5e  

# subset rows using a logical condition 
df6 = df[df$numbers < 10, keep_vars] 
df6 

#      numbers  mixed 
# 1       1     1a 
# 2       2     2b 
# 3       3     3c 
# 4       4     4d 
# 5       5     5e 
# 6       6     6f 
# 7       7     7g 
# 8       8     8h 
# 9       9     9i

大数据分析 - SQL入门

SQL代表结构化查询语言。它是从传统数据仓库和大数据技术中的数据库中提取数据的最广泛使用的语言之一。为了演示SQL的基础知识，我们将使用示例。为了专注于语言本身，我们将在R内部使用SQL。在编写SQL代码方面，这与在数据库中进行的操作完全一样。

SQL的核心是三个语句：SELECT、FROM和WHERE。以下示例使用了SQL最常见的用例。导航到文件夹bda/part2/SQL_introduction并打开SQL_introduction.Rproj文件。然后打开01_select.R脚本。为了在R中编写SQL代码，我们需要安装sqldf包，如下面的代码所示。

# Install the sqldf package
install.packages('sqldf')  

# load the library 
library('sqldf') 
library(nycflights13)  

# We will be working with the fligths dataset in order to introduce SQL  

# Let’s take a look at the table 
str(flights) 
# Classes 'tbl_d', 'tbl' and 'data.frame': 336776 obs. of  16 variables: 

# $ year     : int  2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ... 
# $ month    : int  1 1 1 1 1 1 1 1 1 1 ... 
# $ day      : int  1 1 1 1 1 1 1 1 1 1 ... 
# $ dep_time : int  517 533 542 544 554 554 555 557 557 558 ... 
# $ dep_delay: num  2 4 2 -1 -6 -4 -5 -3 -3 -2 ... 
# $ arr_time : int  830 850 923 1004 812 740 913 709 838 753 ... 
# $ arr_delay: num  11 20 33 -18 -25 12 19 -14 -8 8 ...
# $ carrier  : chr  "UA" "UA" "AA" "B6" ... 

# $ tailnum  : chr  "N14228" "N24211" "N619AA" "N804JB" ... 
# $ flight   : int  1545 1714 1141 725 461 1696 507 5708 79 301 ... 
# $ origin   : chr  "EWR" "LGA" "JFK" "JFK" ... 
# $ dest     : chr  "IAH" "IAH" "MIA" "BQN" ... 
# $ air_time : num  227 227 160 183 116 150 158 53 140 138 ... 
# $ distance : num  1400 1416 1089 1576 762 ... 
# $ hour     : num  5 5 5 5 5 5 5 5 5 5 ... 
# $ minute   : num  17 33 42 44 54 54 55 57 57 58 ...

select语句用于从表中检索列并在其上进行计算。最简单的SELECT语句在ej1中演示。我们还可以创建新的变量，如ej2所示。

### SELECT statement
ej1 = sqldf(" 
   SELECT  
   dep_time 
   ,dep_delay 
   ,arr_time 
   ,carrier 
   ,tailnum 
   FROM 
   flights
")  

head(ej1) 
#    dep_time   dep_delay  arr_time  carrier  tailnum 
# 1      517         2      830      UA       N14228 
# 2      533         4      850      UA       N24211 
# 3      542         2      923      AA       N619AA 
# 4      544        -1     1004      B6       N804JB 
# 5      554        -6      812      DL       N668DN 
# 6      554        -4      740      UA       N39463  

# In R we can use SQL with the sqldf function. It works exactly the same as in 
a database 

# The data.frame (in this case flights) represents the table we are querying 
and goes in the FROM statement  
# We can also compute new variables in the select statement using the syntax: 

# old_variables as new_variable 
ej2 = sqldf(" 
   SELECT 
   arr_delay - dep_delay as gain, 
   carrier 
   FROM 
   flights
")  

ej2[1:5, ] 
#    gain   carrier 
# 1    9      UA 
# 2   16      UA 
# 3   31      AA 
# 4  -17      B6 
# 5  -19      DL

SQL最常用的功能之一是group by语句。这允许计算另一个变量的不同组的数值。打开脚本02_group_by.R。

### GROUP BY      

# Computing the average 
ej3 = sqldf(" 
  SELECT 
   avg(arr_delay) as mean_arr_delay, 
   avg(dep_delay) as mean_dep_delay, 
   carrier 
   FROM 
   flights 
   GROUP BY 
   carrier 
")  

#    mean_arr_delay   mean_dep_delay carrier 
# 1       7.3796692      16.725769      9E 
# 2       0.3642909       8.586016      AA 
# 3      -9.9308886       5.804775      AS 
# 4       9.4579733      13.022522      B6 
# 5       1.6443409       9.264505      DL 
# 6      15.7964311      19.955390      EV 
# 7      21.9207048      20.215543      F9 
# 8      20.1159055      18.726075      FL 
# 9      -6.9152047       4.900585      HA 
# 10     10.7747334      10.552041      MQ
# 11     11.9310345      12.586207      OO 
# 12      3.5580111      12.106073      UA 
# 13      2.1295951       3.782418      US 
# 14      1.7644644      12.869421      VX 
# 15      9.6491199      17.711744      WN 
# 16     15.5569853      18.996330      YV  

# Other aggregations 
ej4 = sqldf(" 
   SELECT 
   avg(arr_delay) as mean_arr_delay, 
   min(dep_delay) as min_dep_delay, 
   max(dep_delay) as max_dep_delay, 
   carrier 
   FROM  
   flights 
   GROUP BY 
   carrier 
")  

# We can compute the minimun, mean, and maximum values of a numeric value 
ej4 
#      mean_arr_delay    min_dep_delay   max_dep_delay   carrier 
# 1       7.3796692           -24           747          9E 
# 2       0.3642909           -24          1014          AA 
# 3      -9.9308886           -21           225          AS 
# 4       9.4579733           -43           502          B6
# 5       1.6443409           -33           960         DL 
# 6      15.7964311           -32           548         EV 
# 7      21.9207048           -27           853         F9 
# 8      20.1159055           -22           602         FL 
# 9      -6.9152047           -16          1301         HA 
# 10     10.7747334           -26          1137         MQ 
# 11     11.9310345           -14           154         OO 
# 12      3.5580111           -20           483         UA 
# 13      2.1295951           -19           500         US 
# 14      1.7644644           -20           653         VX 
# 15      9.6491199           -13           471         WN 
# 16     15.5569853           -16           387         YV  

### We could be also interested in knowing how many observations each carrier has  
ej5 = sqldf(" 
   SELECT 
   carrier, count(*) as count 
   FROM  
   flights 
   GROUP BY 
   carrier 
")  

ej5 
#      carrier  count 
# 1       9E    18460
# 2       AA   32729 
# 3       AS   714 
# 4       B6   54635 
# 5       DL   48110 
# 6       EV   54173 
# 7       F9   685 
# 8       FL   3260 
# 9       HA   342 
# 10      MQ   26397 
# 11      OO   32 
# 12      UA   58665 
# 13      US   20536 
# 14      VX   5162 
# 15      WN   12275 
# 16      YV   601

SQL最有用的功能是连接。连接意味着我们想要使用一列来匹配两个表的数值，将表A和表B组合成一个表。连接有不同的类型，实际上，为了入门，这些将是最有用的：内部连接和左外部连接。

# Let’s create two tables: A and B to demonstrate joins.
A = data.frame(c1 = 1:4, c2 = letters[1:4]) 
B = data.frame(c1 = c(2,4,5,6), c2 = letters[c(2:5)])  

A 
# c1 c2 
# 1  a 
# 2  b 
# 3  c 
# 4  d  

B 
# c1 c2 
# 2  b 
# 4  c 
# 5  d 
# 6  e  

### INNER JOIN 
# This means to match the observations of the column we would join the tables by.   
inner = sqldf(" 
   SELECT 
   A.c1, B.c2 
   FROM 
   A INNER JOIN B 
   ON A.c1 = B.c1 
")  

# Only the rows that match c1 in both A and B are returned 
inner 
# c1 c2 
#  2  b 
#  4  c  

### LEFT OUTER JOIN
# the left outer join, sometimes just called left join will return the  
# first all the values of the column used from the A table  
left = sqldf(" 
  SELECT 
   A.c1, B.c2 
  FROM 
   A LEFT OUTER JOIN B 
   ON A.c1 = B.c1 
")  

# Only the rows that match c1 in both A and B are returned 
left 
#   c1    c2 
#    1  <NA> 
#    2    b 
#    3  <NA> 
#    4    c

大数据分析 - 图表

分析数据的首要方法是直观地分析它。这样做的目标通常是找到变量之间的关系和变量的单变量描述。我们可以将这些策略分为：

单变量分析
多变量分析

单变量图形方法

单变量是一个统计术语。实际上，这意味着我们想要独立于其余数据分析一个变量。能够有效地做到这一点的图是：

箱线图

箱线图通常用于比较分布。这是直观检查不同分布之间是否存在差异的好方法。我们可以看看不同切割方式下钻石价格是否存在差异。

# We will be using the ggplot2 library for plotting
library(ggplot2)  
data("diamonds")  

# We will be using the diamonds dataset to analyze distributions of numeric variables 
head(diamonds) 

#    carat   cut       color  clarity  depth  table   price    x     y     z 
# 1  0.23    Ideal       E      SI2    61.5    55     326     3.95  3.98  2.43 
# 2  0.21    Premium     E      SI1    59.8    61     326     3.89  3.84  2.31 
# 3  0.23    Good        E      VS1    56.9    65     327     4.05  4.07  2.31 
# 4  0.29    Premium     I      VS2    62.4    58     334     4.20  4.23  2.63 
# 5  0.31    Good        J      SI2    63.3    58     335     4.34  4.35  2.75 
# 6  0.24    Very Good   J      VVS2   62.8    57     336     3.94  3.96  2.48 

### Box-Plots
p = ggplot(diamonds, aes(x = cut, y = price, fill = cut)) + 
   geom_box-plot() + 
   theme_bw() 
print(p)

我们可以看到，在图中，不同切割类型的钻石价格分布存在差异。

直方图

source('01_box_plots.R')

# We can plot histograms for each level of the cut factor variable using 
facet_grid 
p = ggplot(diamonds, aes(x = price, fill = cut)) + 
   geom_histogram() + 
   facet_grid(cut ~ .) + 
   theme_bw() 

p  
# the previous plot doesn’t allow to visuallize correctly the data because of 
the differences in scale 
# we can turn this off using the scales argument of facet_grid  

p = ggplot(diamonds, aes(x = price, fill = cut)) + 
   geom_histogram() + 
   facet_grid(cut ~ ., scales = 'free') + 
   theme_bw() 
p  

png('02_histogram_diamonds_cut.png') 
print(p) 
dev.off()

上述代码的输出如下：

多变量图形方法

探索性数据分析中的多变量图形方法的目标是寻找不同变量之间的关系。通常有两种方法可以做到这一点：绘制数值变量的相关矩阵或简单地将原始数据绘制为散点图矩阵。

为了演示这一点，我们将使用diamonds数据集。要遵循代码，请打开脚本bda/part2/charts/03_multivariate_analysis.R。

library(ggplot2)
data(diamonds) 

# Correlation matrix plots  
keep_vars = c('carat', 'depth', 'price', 'table') 
df = diamonds[, keep_vars]  
# compute the correlation matrix 
M_cor = cor(df) 

#          carat       depth      price      table 
# carat 1.00000000  0.02822431  0.9215913  0.1816175 
# depth 0.02822431  1.00000000 -0.0106474 -0.2957785 
# price 0.92159130 -0.01064740  1.0000000  0.1271339 
# table 0.18161755 -0.29577852  0.1271339  1.0000000  

# plots 
heat-map(M_cor)

代码将产生以下输出：

这是一个摘要，它告诉我们价格和克拉之间存在很强的相关性，而其他变量之间则没有太多相关性。

当我们有很多变量时，相关矩阵很有用，在这种情况下，绘制原始数据是不切实际的。如前所述，也可以显示原始数据：

library(GGally)
ggpairs(df)

我们可以看到，在热力图中显示的结果得到了证实，价格和克拉变量之间存在0.922的相关性。

可以在散点图矩阵的(3, 1)索引中找到的价格-克拉散点图中直观地看到这种关系。

大数据分析 - 数据分析工具

各种工具允许数据科学家有效地分析数据。通常，数据分析的工程方面侧重于数据库，数据科学家则侧重于可以实现数据产品的工具。下一节将讨论不同工具的优势，重点是数据科学家在实践中最常用的统计软件包。

R编程语言

R是一种开源编程语言，专注于统计分析。就统计能力而言，它与SAS、SPSS等商业工具具有竞争力。它被认为是其他编程语言（如C、C++或Fortran）的接口。

R的另一个优势是大量的开源库。在CRAN中，有6000多个可以免费下载的包，在Github中，有各种各样的R包可用。

在性能方面，R对于密集型操作来说速度很慢，鉴于大量的库可用，代码的慢速部分是用编译语言编写的。但是，如果你打算进行需要编写深度for循环的操作，那么R就不是你的最佳选择。对于数据分析目的，有一些很好的库，例如data.table, glmnet, ranger, xgboost, ggplot2, caret，它们允许将R用作更快速编程语言的接口。

用于数据分析的Python

Python是一种通用编程语言，它包含大量用于数据分析的库，例如pandas, scikit-learn, theano, numpy和scipy。

R中提供的大多数功能也可以在Python中完成，但我们发现R更易于使用。如果你正在处理大型数据集，通常Python比R更好。Python可以非常有效地逐行清理和处理数据。R也可以做到这一点，但对于脚本任务来说，它不如Python高效。

对于机器学习，scikit-learn是一个不错的环境，它提供了大量可以轻松处理中等大小数据集的算法。与R的等效库(caret)相比，scikit-learn具有更清晰、更一致的API。

Julia

Julia是一种用于技术计算的高级、高性能动态编程语言。它的语法与R或Python非常相似，因此，如果你已经在使用R或Python，那么用Julia编写相同的代码应该非常简单。这门语言相当新，并且在过去几年里发展迅速，所以它绝对是一个不错的选择。

我们建议将Julia用于原型设计计算密集型算法，例如神经网络。它是进行研究的一个好工具。在生产中实现模型方面，Python可能会有更好的选择。然而，随着提供R、Python和Julia模型实现工程的网络服务的出现，这个问题正变得越来越不重要。

SAS

SAS是一种商业语言，仍在用于商业智能。它有一个基础语言，允许用户编程各种各样的应用程序。它包含相当多的商业产品，使非专业用户能够使用复杂工具，例如神经网络库，而无需编程。

除了商业工具的明显缺点之外，SAS无法很好地扩展到大型数据集。即使是中等大小的数据集也会让SAS出现问题，并导致服务器崩溃。只有当你处理小型数据集并且用户不是专业数据科学家时，才推荐使用SAS。对于高级用户来说，R和Python提供了更高效的环境。

SPSS

SPSS是IBM目前用于统计分析的产品。它主要用于分析调查数据，对于无法编程的用户来说，它是一个不错的选择。它可能与SAS一样易于使用，但在实现模型方面，它更简单，因为它提供SQL代码来评分模型。这段代码通常效率不高，但这是一个开始，而SAS则分别为每个数据库销售评分模型的产品。对于小型数据和缺乏经验的团队，SPSS是一个与SAS一样好的选择。

然而，该软件的功能相当有限，经验丰富的用户使用R或Python的效率要高得多。

Matlab，Octave

还有一些其他可用的工具，例如Matlab或其开源版本（Octave）。这些工具主要用于研究。就功能而言，R或Python可以实现Matlab或Octave的所有功能。只有当您对它们提供的支持感兴趣时，购买产品的许可证才有意义。

大数据分析 - 统计方法

在分析数据时，可以使用统计方法。执行基本分析所需的工具包括：

相关性分析
方差分析
假设检验

处理大型数据集时，这不会造成问题，因为这些方法在计算上并不密集，相关性分析除外。在这种情况下，始终可以抽取样本，结果应该是稳健的。

卡方检验

卡方检验允许我们检验两个随机变量是否独立。这意味着每个变量的概率分布不会影响另一个变量。为了在R中评估该检验，我们首先需要创建一个列联表，然后将该表传递给**chisq.test R**函数。

例如，让我们检查diamonds数据集中的cut和color变量之间是否存在关联。该检验正式定义如下：

H0：变量cut和color是独立的
H1：变量cut和color不是独立的

从变量名称来看，我们假设这两个变量之间存在关系，但是该检验可以提供一个客观的“规则”，说明该结果的显著性如何。

在下面的代码片段中，我们发现该检验的p值为2.2e-16，这在实践中几乎为零。然后在进行**蒙特卡罗模拟**后运行该检验，我们发现p值为0.0004998，这仍然远低于阈值0.05。此结果意味着我们拒绝零假设（H0），因此我们认为变量**cut**和**color**不是独立的。

library(ggplot2)

# Use the table function to compute the contingency table 
tbl = table(diamonds$cut, diamonds$color) 
tbl  

#              D    E    F    G    H    I    J 
# Fair       163  224  312  314  303  175  119 
# Good       662  933  909  871  702  522  307 
# Very Good 1513 2400 2164 2299 1824 1204  678 
# Premium   1603 2337 2331 2924 2360 1428  808 
# Ideal     2834 3903 3826 4884 3115 2093  896  

# In order to run the test we just use the chisq.test function. 
chisq.test(tbl)  

# Pearson’s Chi-squared test 
# data:  tbl 
# X-squared = 310.32, df = 24, p-value < 2.2e-16
# It is also possible to compute the p-values using a monte-carlo simulation 
# It's needed to add the simulate.p.value = TRUE flag and the amount of 
simulations 
chisq.test(tbl, simulate.p.value = TRUE, B = 2000)  

# Pearson’s Chi-squared test with simulated p-value (based on 2000 replicates) 
# data:  tbl 
# X-squared = 310.32, df = NA, p-value = 0.0004998

t检验

**t检验**的目的是评估名义变量的不同组之间数值变量的分布是否存在差异。为了证明这一点，我将选择因子变量cut的Fair和Ideal水平，然后我们将比较这两个组中数值变量的值。

data = diamonds[diamonds$cut %in% c('Fair', 'Ideal'), ]

data$cut = droplevels.factor(data$cut) # Drop levels that aren’t used from the 
cut variable 
df1 = data[, c('cut', 'price')]  

# We can see the price means are different for each group 
tapply(df1$price, df1$cut, mean) 
# Fair    Ideal  
# 4358.758 3457.542

t检验在R中使用**t.test**函数实现。t.test的公式接口是最简单的使用方法，其思想是用一个组变量来解释一个数值变量。

例如：**t.test(numeric_variable ~ group_variable, data = data)**。在前面的示例中，**numeric_variable**是**price**，**group_variable**是**cut**。

从统计学的角度来看，我们正在检验数值变量在两组之间的分布是否存在差异。正式的假设检验是用零假设（H0）和备择假设（H1）来描述的。

H0：Fair和Ideal组之间价格变量的分布没有差异
H1：Fair和Ideal组之间价格变量的分布存在差异

以下可以在R中使用以下代码实现：

t.test(price ~ cut, data = data)

# Welch Two Sample t-test 
#  
# data:  price by cut 
# t = 9.7484, df = 1894.8, p-value < 2.2e-16 
# alternative hypothesis: true difference in means is not equal to 0 
# 95 percent confidence interval: 
#   719.9065 1082.5251 
# sample estimates: 
#   mean in group Fair mean in group Ideal  
#   4358.758            3457.542   

# Another way to validate the previous results is to just plot the 
distributions using a box-plot 
plot(price ~ cut, data = data, ylim = c(0,12000),  
   col = 'deepskyblue3')

我们可以通过检查p值是否低于0.05来分析检验结果。如果是这种情况，我们将保留备择假设。这意味着我们在cut因子的两个水平之间发现了价格差异。从水平名称来看，我们本应该预料到这个结果，但是我们没有预料到Fail组的平均价格会高于Ideal组。我们可以通过比较每个因子的均值来观察这一点。

**plot**命令生成一个图形，显示价格和cut变量之间的关系。这是一个箱线图；我们在16.0.1节中介绍过这种图，但它基本上显示了我们正在分析的cut的两个水平的价格变量的分布。

方差分析

方差分析（ANOVA）是一种统计模型，用于通过比较每个组的均值和方差来分析组分布之间的差异，该模型由罗纳德·费舍尔开发。ANOVA提供了一个统计检验，用于检验多个组的均值是否相等，因此它将t检验推广到两个以上组。

ANOVA对于比较三个或更多组的统计显著性非常有用，因为进行多个双样本t检验会导致犯I类统计错误的可能性增加。

就提供数学解释而言，需要理解以下内容才能理解该检验。

x_ij = x + (x_i − x) + (x_ij − x)

这导致以下模型：

x_ij = μ + α_i + ∈_ij

其中μ是总均值，α_i是第i个组的均值。误差项∈_ij假定是从正态分布中独立同分布的。该检验的零假设是：

α₁ = α₂ = … = α_k

就计算检验统计量而言，我们需要计算两个值：

组间差异的平方和：

$$SSD_B = \sum_{i}^{k} \sum_{j}^{n}(\bar{x_{\bar{i}}} - \bar{x})^2$$

组内平方和

$$SSD_W = \sum_{i}^{k} \sum_{j}^{n}(\bar{x_{\bar{ij}}} - \bar{x_{\bar{i}}})^2$$

其中SSD_B的自由度为k−1，SSD_W的自由度为N−k。然后我们可以定义每个度量的均方差。

MS_B = SSD_B / (k - 1)

MS_w = SSD_w / (N - k)

最后，ANOVA中的检验统计量定义为上述两个量的比率

F = MS_B / MS_w

它服从自由度为k−1和N−k的F分布。如果零假设为真，F值可能接近1。否则，组间均方MSB可能较大，导致F值较大。

基本上，ANOVA检查总方差的两个来源，并查看哪个部分贡献更大。这就是为什么它被称为方差分析，尽管目的是比较组均值。

就计算统计量而言，在R中实际上相当简单。以下示例将演示如何执行此操作以及如何绘制结果。

library(ggplot2)
# We will be using the mtcars dataset 

head(mtcars) 
#                    mpg  cyl disp  hp drat  wt  qsec   vs am  gear carb 
# Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4 
# Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4 
# Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1 
# Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1 
# Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2 
# Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1  

# Let's see if there are differences between the groups of cyl in the mpg variable. 
data = mtcars[, c('mpg', 'cyl')]  
fit = lm(mpg ~ cyl, data = mtcars) 
anova(fit)  

# Analysis of Variance Table 
# Response: mpg 
#           Df Sum Sq Mean Sq F value    Pr(>F)     
# cyl        1 817.71  817.71  79.561 6.113e-10 *** 
# Residuals 30 308.33   10.28 
# Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 
# Plot the distribution 
plot(mpg ~ as.factor(cyl), data = mtcars, col = 'deepskyblue3')

代码将产生以下输出：

我们在示例中获得的p值明显小于0.05，因此R返回符号'***'来表示这一点。这意味着我们拒绝零假设，并且我们在cyl变量的不同组之间发现了mpg均值的差异。

机器学习用于数据分析

机器学习是计算机科学的一个子领域，它处理模式识别、计算机视觉、语音识别、文本分析等任务，并且与统计学和数学优化有着密切的联系。应用包括搜索引擎的开发、垃圾邮件过滤、光学字符识别（OCR）等。数据挖掘、模式识别和统计学习领域之间的界限并不清晰，基本上都指类似的问题。

机器学习可以分为两种类型的任务：

监督学习
无监督学习

监督学习

监督学习指的是一种问题类型，其中输入数据定义为矩阵X，我们感兴趣的是预测响应y。其中X = {x₁, x₂, …, x_n}具有n个预测变量，并且具有两个值y = {c₁, c₂}。

一个示例应用程序是使用人口统计特征作为预测变量来预测网页用户点击广告的概率。这通常被称为预测点击率（CTR）。然后y = {点击，不点击}，预测变量可以是使用的IP地址、用户进入网站的日期、用户的城市、国家以及其他可能可用的特征。

无监督学习

无监督学习处理在没有要学习的类别的情况下查找彼此相似的组的问题。有多种方法可以将预测变量映射到查找在每个组中共享相似实例并在彼此之间不同的组。

无监督学习的一个示例应用程序是客户细分。例如，在电信行业中，一个常见的任务是根据用户对电话的使用情况对用户进行细分。这将允许营销部门针对每个组使用不同的产品。

大数据分析 - 朴素贝叶斯分类器

朴素贝叶斯是一种用于构建分类器的概率技术。朴素贝叶斯分类器的特征假设是认为，给定类别变量，特定特征的值独立于任何其他特征的值。

尽管前面提到了过于简化的假设，但朴素贝叶斯分类器在复杂的现实世界情况中具有良好的结果。朴素贝叶斯的优点是它只需要少量训练数据来估计分类所需的必要参数，并且可以增量地训练分类器。

朴素贝叶斯是一个条件概率模型：给定一个要分类的问题实例，该实例由表示一些n个特征（自变量）的向量x = (x₁, …, x_n)表示，它为每个K个可能的输出或类别分配此实例的概率。

$$p(C_k|x_1,....., x_n)$$

上述公式的问题在于，如果特征数量 n 很大，或者某个特征可以取很多值，那么基于概率表构建这样的模型是不可行的。因此，我们重新构建模型使其更简单。利用贝叶斯定理，条件概率可以分解为：

$$p(C_k|x) = \frac{p(C_k)p(x|C_k)}{p(x)}$$

这意味着在上述独立性假设下，类变量 C 的条件分布为：

$$p(C_k|x_1,....., x_n)\: = \: \frac{1}{Z}p(C_k)\prod_{i = 1}^{n}p(x_i|C_k)$$

其中证据 Z = p(x) 是一个仅依赖于 x₁, …, x_n 的比例因子，如果特征变量的值已知，则为常数。一个常见的规则是选择最可能的假设；这被称为最大后验概率或 MAP 决策规则。相应的分类器，即贝叶斯分类器，是一个函数，它如下分配类标签 $\hat{y} = C_k$ (其中 k 为某值)：

$$\hat{y} = argmax\: p(C_k)\prod_{i = 1}^{n}p(x_i|C_k)$$

在 R 中实现该算法是一个简单的过程。以下示例演示了如何训练朴素贝叶斯分类器并将其用于垃圾邮件过滤问题中的预测。

以下脚本位于 bda/part3/naive_bayes/naive_bayes.R 文件中。

# Install these packages 
pkgs = c("klaR", "caret", "ElemStatLearn") 
install.packages(pkgs)  
library('ElemStatLearn') 
library("klaR") 
library("caret")  

# Split the data in training and testing 
inx = sample(nrow(spam), round(nrow(spam) * 0.9)) 
train = spam[inx,] 
test = spam[-inx,]  

# Define a matrix with features, X_train 
# And a vector with class labels, y_train 
X_train = train[,-58] 
y_train = train$spam  
X_test = test[,-58] 
y_test = test$spam  
# Train the model 
nb_model = train(X_train, y_train, method = 'nb',  
   trControl = trainControl(method = 'cv', number = 3)) 

# Compute  
preds = predict(nb_model$finalModel, X_test)$class 
tbl = table(y_test, yhat = preds) 
sum(diag(tbl)) / sum(tbl) 
# 0.7217391

从结果可以看出，朴素贝叶斯模型的准确率为 72%。这意味着该模型正确分类了 72% 的实例。

大数据分析 - K 均值聚类

K 均值聚类旨在将 n 个观测值划分为 k 个簇，其中每个观测值都属于具有最近均值的簇，该均值作为簇的原型。这导致数据空间被划分为 Voronoi 单元。

给定一组观测值 (x₁, x₂, …, x_n)，其中每个观测值都是一个 d 维实数向量，K 均值聚类旨在将 n 个观测值划分为 k 个组 G = {G₁, G₂, …, G_k}，以最小化如下定义的簇内平方和 (WCSS)：

$$argmin \: \sum_{i = 1}^{k} \sum_{x \in S_{i}}\parallel x - \mu_{i}\parallel ^2$$

后面的公式显示了为了找到 K 均值聚类中的最佳原型而最小化的目标函数。该公式的直觉是，我们希望找到彼此不同的组，并且每个组的每个成员都应该与该簇的其他成员相似。

以下示例演示了如何在 R 中运行 K 均值聚类算法。

library(ggplot2)
# Prepare Data 
data = mtcars  

# We need to scale the data to have zero mean and unit variance 
data <- scale(data)  

# Determine number of clusters 
wss <- (nrow(data)-1)*sum(apply(data,2,var)) 
for (i in 2:dim(data)[2]) { 
   wss[i] <- sum(kmeans(data, centers = i)$withinss) 
}  

# Plot the clusters 
plot(1:dim(data)[2], wss, type = "b", xlab = "Number of Clusters", 
   ylab = "Within groups sum of squares")

为了找到 K 的一个好值，我们可以绘制不同 K 值的组内平方和。此度量通常随着添加更多组而减少，我们希望找到组内平方和的减少开始缓慢减少的点。在图中，此值最好由 K = 6 表示。

现在 K 值已定义，需要使用该值运行算法。

# K-Means Cluster Analysis
fit <- kmeans(data, 5) # 5 cluster solution 

# get cluster means  
aggregate(data,by = list(fit$cluster),FUN = mean) 

# append cluster assignment 
data <- data.frame(data, fit$cluster)

大数据分析 - 关联规则

设 I = i₁, i₂, ..., i_n 为一组 n 个二元属性，称为项目。设 D = t₁, t₂, ..., t_m 为一组事务，称为数据库。D 中的每个事务都有一个唯一的事务 ID，并包含 I 中项目的子集。规则定义为 X ⇒ Y 形式的蕴涵，其中 X, Y ⊆ I 且 X ∩ Y = ∅。

项目集（简称项目集）X 和 Y 分别称为规则的前件（左侧或 LHS）和后件（右侧或 RHS）。

为了说明这些概念，我们使用超市领域的一个小例子。项目集为 I = {牛奶、面包、黄油、啤酒}，包含这些项目的小型数据库如下表所示。

事务 ID	项目
1	牛奶、面包
2	面包、黄油
3	啤酒
4	牛奶、面包、黄油
5	面包、黄油

超市的一个示例规则可能是 {牛奶、面包} ⇒ {黄油}，这意味着如果购买了牛奶和面包，顾客也会购买黄油。为了从所有可能的规则集中选择有趣的规则，可以使用对各种显著性和兴趣度量的约束。最著名的约束是支持和置信度的最小阈值。

项目集 X 的支持 supp(X) 定义为数据集中包含该项目集的事务的比例。在表 1 的示例数据库中，项目集 {牛奶、面包} 的支持度为 2/5 = 0.4，因为它出现在 40% 的所有事务中（5 个事务中的 2 个）。查找频繁项目集可以看作是无监督学习问题的简化。

规则的置信度定义为 conf(X ⇒ Y ) = supp(X ∪ Y )/supp(X)。例如，规则 {牛奶、面包} ⇒ {黄油} 在表 1 的数据库中的置信度为 0.2/0.4 = 0.5，这意味着对于包含牛奶和面包的事务的 50%，该规则是正确的。置信度可以解释为概率 P(Y|X) 的估计值，即在这些事务也包含 LHS 的条件下，在事务中找到规则 RHS 的概率。

在位于 bda/part3/apriori.R 中的脚本中，可以找到实现 Apriori 算法 的代码。

# Load the library for doing association rules
# install.packages(’arules’) 
library(arules)  

# Data preprocessing 
data("AdultUCI") 
AdultUCI[1:2,]  
AdultUCI[["fnlwgt"]] <- NULL 
AdultUCI[["education-num"]] <- NULL  

AdultUCI[[ "age"]] <- ordered(cut(AdultUCI[[ "age"]], c(15,25,45,65,100)), 
   labels = c("Young", "Middle-aged", "Senior", "Old")) 
AdultUCI[[ "hours-per-week"]] <- ordered(cut(AdultUCI[[ "hours-per-week"]], 
   c(0,25,40,60,168)), labels = c("Part-time", "Full-time", "Over-time", "Workaholic")) 
AdultUCI[[ "capital-gain"]] <- ordered(cut(AdultUCI[[ "capital-gain"]], 
   c(-Inf,0,median(AdultUCI[[ "capital-gain"]][AdultUCI[[ "capitalgain"]]>0]),Inf)), 
   labels = c("None", "Low", "High")) 
AdultUCI[[ "capital-loss"]] <- ordered(cut(AdultUCI[[ "capital-loss"]], 
   c(-Inf,0, median(AdultUCI[[ "capital-loss"]][AdultUCI[[ "capitalloss"]]>0]),Inf)), 
   labels = c("none", "low", "high"))

为了使用 Apriori 算法生成规则，我们需要创建一个事务矩阵。以下代码显示了如何在 R 中执行此操作。

# Convert the data into a transactions format
Adult <- as(AdultUCI, "transactions") 
Adult 
# transactions in sparse format with 
# 48842 transactions (rows) and 
# 115 items (columns)  

summary(Adult)  
# Plot frequent item-sets 
itemFrequencyPlot(Adult, support = 0.1, cex.names = 0.8)  

# generate rules 
min_support = 0.01 
confidence = 0.6 
rules <- apriori(Adult, parameter = list(support = min_support, confidence = confidence))

rules 
inspect(rules[100:110, ]) 
# lhs                             rhs                      support     confidence  lift
# {occupation = Farming-fishing} => {sex = Male}        0.02856148  0.9362416   1.4005486
# {occupation = Farming-fishing} => {race = White}      0.02831579  0.9281879   1.0855456
# {occupation = Farming-fishing} => {native-country     0.02671881  0.8758389   0.9759474
                                       = United-States}

大数据分析 - 决策树

决策树是一种用于监督学习问题的算法，例如分类或回归。决策树或分类树是一棵树，其中每个内部（非叶）节点都用输入特征标记。来自用特征标记的节点的弧线用特征的每个可能值标记。树的每个叶子都用一个类或一个类上的概率分布标记。

可以通过基于属性值测试将源集拆分为子集来“学习”树。此过程以递归方式重复应用于每个派生子集，称为递归划分。当节点处的子集具有目标变量的相同值，或者拆分不再增加预测的价值时，递归完成。这种自上而下归纳决策树的过程是贪婪算法的一个例子，也是学习决策树最常见的策略。

数据挖掘中使用的决策树主要有两种：

分类树 - 当响应是名义变量时，例如电子邮件是否是垃圾邮件。
回归树 - 当预测结果可以被认为是一个实数时（例如，工人的工资）。

决策树是一种简单的方法，因此存在一些问题。其中一个问题是决策树产生的结果模型的高方差。为了减轻这个问题，开发了决策树的集成方法。目前广泛使用的集成方法有两类：

Bagging 决策树 - 通过重复地用替换方法重新采样训练数据，并对树进行投票以获得一致性预测，来构建多个决策树。该算法被称为随机森林。
Boosting 决策树 - 梯度提升将弱学习器组合；在这种情况下，将决策树以迭代的方式组合成一个强大的学习器。它将一个弱树拟合到数据，并迭代地继续拟合弱学习器以纠正先前模型的错误。

# Install the party package
# install.packages('party') 
library(party) 
library(ggplot2)  

head(diamonds) 
# We will predict the cut of diamonds using the features available in the 
diamonds dataset. 
ct = ctree(cut ~ ., data = diamonds) 

# plot(ct, main="Conditional Inference Tree") 
# Example output 
# Response:  cut  
# Inputs:  carat, color, clarity, depth, table, price, x, y, z  

# Number of observations:  53940  
#  
# 1) table <= 57; criterion = 1, statistic = 10131.878 
#   2) depth <= 63; criterion = 1, statistic = 8377.279 
#     3) table <= 56.4; criterion = 1, statistic = 226.423 
#       4) z <= 2.64; criterion = 1, statistic = 70.393 
#         5) clarity <= VS1; criterion = 0.989, statistic = 10.48 
#           6) color <= E; criterion = 0.997, statistic = 12.829 
#             7)*  weights = 82  
#           6) color > E  

#Table of prediction errors 
table(predict(ct), diamonds$cut) 
#            Fair  Good Very Good Premium Ideal 
# Fair       1388   171        17       0    14 
# Good        102  2912       499      26    27 
# Very Good    54   998      3334     249   355 
# Premium      44   711      5054   11915  1167 
# Ideal        22   114      3178    1601 19988 
# Estimated class probabilities 
probs = predict(ct, newdata = diamonds, type = "prob") 
probs = do.call(rbind, probs) 
head(probs)

大数据分析 - 逻辑回归

逻辑回归是一种分类模型，其中响应变量是分类的。这是一种来自统计学的算法，用于监督分类问题。在逻辑回归中，我们试图找到以下等式中参数向量 β，以最小化成本函数。

$$logit(p_i) = ln \left ( \frac{p_i}{1 - p_i} \right ) = \beta_0 + \beta_1x_{1,i} + ... + \beta_kx_{k,i}$$

以下代码演示了如何在 R 中拟合逻辑回归模型。我们将在这里使用与朴素贝叶斯相同的垃圾邮件数据集来演示逻辑回归。

从准确性方面的预测结果来看，我们发现回归模型在测试集中的准确率达到了 92.5%，而朴素贝叶斯分类器的准确率为 72%。

library(ElemStatLearn)
head(spam) 

# Split dataset in training and testing 
inx = sample(nrow(spam), round(nrow(spam) * 0.8)) 
train = spam[inx,] 
test = spam[-inx,]  

# Fit regression model 
fit = glm(spam ~ ., data = train, family = binomial()) 
summary(fit)  

# Call: 
#   glm(formula = spam ~ ., family = binomial(), data = train) 
#  

# Deviance Residuals:  
#   Min       1Q   Median       3Q      Max   
# -4.5172  -0.2039   0.0000   0.1111   5.4944
# Coefficients: 
# Estimate Std. Error z value Pr(>|z|)     
# (Intercept) -1.511e+00  1.546e-01  -9.772  < 2e-16 *** 
# A.1         -4.546e-01  2.560e-01  -1.776 0.075720 .   
# A.2         -1.630e-01  7.731e-02  -2.108 0.035043 *   
# A.3          1.487e-01  1.261e-01   1.179 0.238591     
# A.4          2.055e+00  1.467e+00   1.401 0.161153     
# A.5          6.165e-01  1.191e-01   5.177 2.25e-07 *** 
# A.6          7.156e-01  2.768e-01   2.585 0.009747 **  
# A.7          2.606e+00  3.917e-01   6.652 2.88e-11 *** 
# A.8          6.750e-01  2.284e-01   2.955 0.003127 **  
# A.9          1.197e+00  3.362e-01   3.559 0.000373 *** 
# Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1  1  

### Make predictions 
preds = predict(fit, test, type = ’response’) 
preds = ifelse(preds > 0.5, 1, 0) 
tbl = table(target = test$spam, preds) 
tbl 

#         preds 
# target    0   1 
# email   535  23 
# spam     46 316 
sum(diag(tbl)) / sum(tbl) 
# 0.925

大数据分析 - 时间序列分析

时间序列是以日期或时间戳为索引的分类或数值变量的一系列观测值。股票价格的时间序列是时间序列数据的清晰示例。在下表中，我们可以看到时间序列数据的基本结构。在这种情况下，每小时记录一次观测值。

时间戳	股票价格
2015-10-11 09:00:00	100
2015-10-11 10:00:00	110
2015-10-11 11:00:00	105
2015-10-11 12:00:00	90
2015-10-11 13:00:00	120

通常，时间序列分析的第一步是绘制序列，这通常使用折线图完成。

时间序列分析最常见的应用是使用数据的时态结构预测数值变量的未来值。这意味着，使用可用的观测值来预测未来的值。

数据的时态排序意味着传统的回归方法没有用。为了构建稳健的预测，我们需要考虑数据的时态排序的模型。

时间序列分析中最广泛使用的模型称为自回归移动平均 (ARMA)。该模型由两部分组成：自回归 (AR) 部分和移动平均 (MA) 部分。然后，该模型通常被称为ARMA(p, q) 模型，其中p 是自回归部分的阶数，q 是移动平均部分的阶数。

自回归模型

AR(p) 读作 p 阶自回归模型。数学上写成：

$$X_t = c + \sum_{i = 1}^{P} \phi_i X_{t - i} + \varepsilon_{t}$$

其中 {φ₁, …, φ_p} 是要估计的参数，c 是一个常数，随机变量 ε_t 表示白噪声。参数值需要一些约束，以便模型保持平稳。

移动平均

MA(q) 表示 q 阶移动平均模型：

$$X_t = \mu + \varepsilon_t + \sum_{i = 1}^{q} \theta_i \varepsilon_{t - i}$$

其中 θ₁, ..., θ_q 是模型的参数，μ 是 X_t 的期望值，ε_t, ε_{t − 1}, ... 是白噪声误差项。

自回归移动平均

ARMA(p, q) 模型结合了 p 个自回归项和 q 个移动平均项。数学上，该模型用以下公式表示：

$$X_t = c + \varepsilon_t + \sum_{i = 1}^{P} \phi_iX_{t - 1} + \sum_{i = 1}^{q} \theta_i \varepsilon_{t-i}$$

我们可以看到，ARMA(p, q) 模型是AR(p) 和MA(q) 模型的组合。

为了更好地理解模型，可以考虑方程中的AR部分旨在估计X_{t − i}观测值的参数，以预测变量X_t的值。最终，它是一个过去值的加权平均值。MA部分采用相同的方法，但使用的是先前观测值的误差ε_{t − i}。因此，最终模型的结果是一个加权平均值。

下面的代码片段演示了如何在R中实现ARMA(p, q)。

# install.packages("forecast")
library("forecast")  

# Read the data 
data = scan('fancy.dat') 
ts_data <- ts(data, frequency = 12, start = c(1987,1)) 
ts_data  
plot.ts(ts_data)

绘制数据通常是第一步，以找出数据中是否存在时间结构。从图中可以看出，每年的年底都有很强的峰值。

下面的代码将ARMA模型拟合到数据。它运行多个模型组合，并选择误差较小的模型。

# Fit the ARMA model
fit = auto.arima(ts_data) 
summary(fit) 

# Series: ts_data  
# ARIMA(1,1,1)(0,1,1)[12]                     
#    Coefficients: 
#    ar1     ma1    sma1 
# 0.2401  -0.9013  0.7499 
# s.e.  0.1427   0.0709  0.1790 

#  
# sigma^2 estimated as 15464184:  log likelihood = -693.69 
# AIC = 1395.38   AICc = 1395.98   BIC = 1404.43 

# Training set error measures: 
#                 ME        RMSE      MAE        MPE        MAPE      MASE       ACF1 
# Training set   328.301  3615.374  2171.002  -2.481166  15.97302  0.4905797 -0.02521172

大数据分析 - 文本分析

在本章中，我们将使用本书第一部分中抓取的数据。数据包含描述自由职业者个人资料及其以美元计价的小时费率的文本。下一节的目的是拟合一个模型，根据自由职业者的技能，我们可以预测其每小时的工资。

下面的代码展示了如何将此案例中包含用户技能的原始文本转换为词袋矩阵。为此，我们使用一个名为tm的R库。这意味着，对于语料库中的每个单词，我们都会创建一个变量，其中包含每个变量出现的次数。

library(tm)
library(data.table)  

source('text_analytics/text_analytics_functions.R') 
data = fread('text_analytics/data/profiles.txt') 
rate = as.numeric(data$rate) 
keep = !is.na(rate) 
rate = rate[keep]  

### Make bag of words of title and body 
X_all = bag_words(data$user_skills[keep]) 
X_all = removeSparseTerms(X_all, 0.999) 
X_all 

# <<DocumentTermMatrix (documents: 389, terms: 1422)>> 
#   Non-/sparse entries: 4057/549101 
# Sparsity           : 99% 
# Maximal term length: 80 
# Weighting          : term frequency - inverse document frequency (normalized) (tf-idf) 

### Make a sparse matrix with all the data 
X_all <- as_sparseMatrix(X_all)

现在我们已经将文本表示为稀疏矩阵，我们可以拟合一个能够提供稀疏解的模型。对于这种情况，一个不错的选择是使用LASSO（最小绝对收缩和选择算子）。这是一个能够选择最相关的特征来预测目标的回归模型。

train_inx = 1:200
X_train = X_all[train_inx, ] 
y_train = rate[train_inx]  
X_test = X_all[-train_inx, ] 
y_test = rate[-train_inx]  

# Train a regression model 
library(glmnet) 
fit <- cv.glmnet(x = X_train, y = y_train,  
   family = 'gaussian', alpha = 1,  
   nfolds = 3, type.measure = 'mae') 
plot(fit)  

# Make predictions 
predictions = predict(fit, newx = X_test) 
predictions = as.vector(predictions[,1]) 
head(predictions)  

# 36.23598 36.43046 51.69786 26.06811 35.13185 37.66367 
# We can compute the mean absolute error for the test data 
mean(abs(y_test - predictions)) 
# 15.02175

现在我们有一个模型，可以根据一组技能预测自由职业者的每小时工资。如果收集更多数据，模型的性能将会提高，但是实现此流程的代码将保持不变。

大数据分析 - 在线学习

在线学习是机器学习的一个子领域，它允许将监督学习模型扩展到海量数据集。其基本思想是，我们不需要将所有数据都读入内存中才能拟合模型，只需要一次读取每个实例。

在本例中，我们将展示如何使用逻辑回归实现在线学习算法。与大多数监督学习算法一样，存在一个被最小化的成本函数。在逻辑回归中，成本函数定义为：

$$J(\theta) \: = \: \frac{-1}{m} \left [ \sum_{i = 1}^{m}y^{(i)}log(h_{\theta}(x^{(i)})) + (1 - y^{(i)}) log(1 - h_{\theta}(x^{(i)})) \right ]$$

其中J(θ)表示成本函数，h_θ(x)表示假设。在逻辑回归的情况下，它由以下公式定义：

$$h_\theta(x) = \frac{1}{1 + e^{\theta^T x}}$$

现在我们已经定义了成本函数，我们需要找到一个算法来最小化它。实现此目的最简单的算法称为随机梯度下降。逻辑回归模型权重的算法更新规则定义为：

$$\theta_j : = \theta_j - \alpha(h_\theta(x) - y)x$$

此算法有几种实现方式，但是vowpal wabbit库中的实现是迄今为止最完善的。该库允许训练大规模回归模型，并使用少量RAM。用创建者自己的话说，它被描述为：“Vowpal Wabbit (VW) 项目是一个快速的out-of-core学习系统，由微软研究院和（以前）雅虎研究院赞助”。

我们将使用来自kaggle竞赛的泰坦尼克号数据集。原始数据可以在bda/part3/vw文件夹中找到。这里，我们有两个文件：

我们有训练数据(train_titanic.csv)，以及
未标记数据，以便进行新的预测(test_titanic.csv)。

为了将csv格式转换为vowpal wabbit输入格式，请使用csv_to_vowpal_wabbit.py python脚本。显然，你需要安装python才能做到这一点。导航到bda/part3/vw文件夹，打开终端并执行以下命令：

python csv_to_vowpal_wabbit.py

请注意，对于本节，如果您使用的是Windows，则需要安装Unix命令行，请访问cygwin网站。

打开终端，也在bda/part3/vw文件夹中，并执行以下命令：

vw train_titanic.vw -f model.vw --binary --passes 20 -c -q ff --sgd --l1 
0.00000001 --l2 0.0000001 --learning_rate 0.5 --loss_function logistic

让我们分解一下vw调用的每个参数的含义。

-f model.vw − 表示我们将模型保存在model.vw文件中，以便以后进行预测
--binary − 将损失报告为具有-1,1标签的二元分类
--passes 20 − 数据被使用20次来学习权重
-c − 创建缓存文件
-q ff − 在f命名空间中使用二次特征
--sgd − 使用常规/经典/简单的随机梯度下降更新，即非自适应、非归一化和非不变。
--l1 --l2 − L1和L2范数正则化
--learning_rate 0.5 − 更新规则公式中定义的学习率α

下面的代码显示了在命令行中运行回归模型的结果。在结果中，我们得到了平均对数损失和算法性能的小型报告。

-loss_function logistic
creating quadratic features for pairs: ff  
using l1 regularization = 1e-08 
using l2 regularization = 1e-07 

final_regressor = model.vw 
Num weight bits = 18 
learning rate = 0.5 
initial_t = 1 
power_t = 0.5 
decay_learning_rate = 1 
using cache_file = train_titanic.vw.cache 
ignoring text input in favor of cache input 
num sources = 1 

average    since         example   example  current  current  current 
loss       last          counter   weight    label   predict  features 
0.000000   0.000000          1      1.0    -1.0000   -1.0000       57 
0.500000   1.000000          2      2.0     1.0000   -1.0000       57 
0.250000   0.000000          4      4.0     1.0000    1.0000       57 
0.375000   0.500000          8      8.0    -1.0000   -1.0000       73 
0.625000   0.875000         16     16.0    -1.0000    1.0000       73 
0.468750   0.312500         32     32.0    -1.0000   -1.0000       57 
0.468750   0.468750         64     64.0    -1.0000    1.0000       43 
0.375000   0.281250        128    128.0     1.0000   -1.0000       43 
0.351562   0.328125        256    256.0     1.0000   -1.0000       43 
0.359375   0.367188        512    512.0    -1.0000    1.0000       57 
0.274336   0.274336       1024   1024.0    -1.0000   -1.0000       57 h 
0.281938   0.289474       2048   2048.0    -1.0000   -1.0000       43 h 
0.246696   0.211454       4096   4096.0    -1.0000   -1.0000       43 h 
0.218922   0.191209       8192   8192.0     1.0000    1.0000       43 h 

finished run 
number of examples per pass = 802 
passes used = 11 
weighted example sum = 8822 
weighted label sum = -2288 
average loss = 0.179775 h 
best constant = -0.530826 
best constant’s loss = 0.659128 
total feature number = 427878

现在我们可以使用我们训练的model.vw来使用新数据生成预测。

vw -d test_titanic.vw -t -i model.vw -p predictions.txt

前面命令生成的预测并未归一化到[0, 1]范围。为此，我们使用sigmoid变换。

# Read the predictions
preds = fread('vw/predictions.txt')  

# Define the sigmoid function 
sigmoid = function(x) { 
   1 / (1 + exp(-x)) 
} 
probs = sigmoid(preds[[1]])  

# Generate class labels 
preds = ifelse(probs > 0.5, 1, 0) 
head(preds) 
# [1] 0 1 0 0 1 0

打印页面