Jupyter-数据科学-全-

Jupyter 数据科学（全）

原文：annas-archive.org/md5/0029981e1e9d3df96c10d26f20813ae8

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

Jupyter 是一个正在广泛使用的开放平台，支持多种编程语言编写程序。这些语言中许多都与数据科学相关，例如 R 和 Python。在本书中，我们将探讨如何在 Jupyter 平台上使用多种语言解决数据科学问题。

我们将从 Jupyter 的基础知识开始学习。然后我们将使用 Jupyter 作为我们的数据分析和可视化平台。我们将探讨数据挖掘、数据清理和机器学习，所有这些都将在 Jupyter 框架的支持下进行。

你将学习如何使用 Jupyter，通过一套编程语言解决数据科学问题。

本书所需内容

本书聚焦于使用 Jupyter 作为数据科学的主要平台。它假设你已经对数据科学概念有较好的理解，并希望使用 Jupyter 作为展示平台。

适用对象

本书适用于那些希望在保持研究本质的同时，公开展示自己发现的数据科学从业者。通过 Jupyter，你可以以交互的方式展示你的具体方法论。

约定

在本书中，你会看到多种文本样式，用来区分不同类型的信息。以下是这些样式的一些示例及其含义的解释。文中的代码词、数据库表名、文件夹名称、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 账号如下所示：“类似地，前面的 describe 语句为我们提供了数据框的一些快速统计信息。”

一块代码设置如下：

plt.xlabel("Actual Price") plt.ylabel("Predicted Price") plt.title("Actual Price vs Predicted Price")

新术语 和 重要词汇 以粗体显示。

你在屏幕上看到的文字，例如在菜单或对话框中，出现在文本中如下所示：“‘Running’ 选项卡列出了已经启动的笔记本。”

警告或重要提示如下所示。

提示和技巧如下所示。

读者反馈

我们始终欢迎读者的反馈。请告诉我们你对本书的看法——你喜欢或不喜欢的地方。读者的反馈对我们来说非常重要，它帮助我们开发出真正能让你受益的书籍。如果你有任何一般性反馈，可以通过电子邮件发送至 feedback@packtpub.com，并在邮件主题中提及书名。如果你在某个领域有专长，并且有兴趣撰写或参与编写一本书，请参阅我们的作者指南：www.packtpub.com/authors。

客户支持

现在你已经拥有了一本 Packt 出版的书籍，我们为你准备了一些内容，帮助你从购买中获得更多收益。

下载示例代码

您可以从您的账户在www.packtpub.com下载本书的示例代码文件。如果您在其他地方购买了此书，您可以访问www.packtpub.com/support并注册将文件直接通过电子邮件发送给您。您可以按照以下步骤下载代码文件：

使用您的电子邮件地址和密码登录或注册我们的网站。
将鼠标指针悬停在顶部的 SUPPORT 标签上。
点击代码下载和勘误。
在搜索框中输入书籍名称。
选择您希望下载代码文件的书籍。
从下拉菜单中选择您购买这本书的地方。
点击代码下载。

文件下载后，请确保使用以下最新版本解压或提取文件夹：

WinRAR / 7-Zip for Windows
Zipeg / iZip / UnRarX for Mac
7-Zip / PeaZip for Linux

这本书的代码包也托管在 GitHub 上，网址是github.com/PacktPublishing/Jupyter-for-Data-Science。我们还从我们丰富的图书和视频目录中提供了其他代码包，您可以访问github.com/PacktPublishing/查看！

勘误

尽管我们已经尽力确保内容的准确性，但错误还是可能发生。如果您在我们的书籍中发现任何错误——可能是文本或代码的错误——我们将非常感激您能向我们报告。通过这样做，您可以帮助其他读者避免困扰，并帮助我们改进后续版本的书籍。如果您发现任何勘误，请访问www.packtpub.com/submit-errata报告，选择您的书籍，点击勘误提交表格链接，并输入勘误的详细信息。一旦您的勘误被验证，您的提交将被接受，勘误将上传到我们的网站或添加到该书籍勘误列表中。要查看之前提交的勘误，请访问www.packtpub.com/books/content/support并在搜索框中输入书籍名称。所需信息将出现在勘误部分。

盗版

互联网上的版权盗版问题在所有媒体中都存在。我们在 Packt 非常重视保护我们的版权和许可。如果您在互联网上发现我们作品的任何非法副本，请立即提供地址或网站名称，以便我们采取措施解决问题。请通过电子邮件copyright@packtpub.com并附上涉嫌盗版内容的链接与我们联系。感谢您帮助我们保护作者和我们的内容。

问题

如果你在本书的任何方面遇到问题，可以通过questions@packtpub.com联系我们，我们会尽力解决问题。

第一章：Jupyter 与数据科学

Jupyter 产品源自 IPython 项目。IPython 项目用于提供 Python 的交互式在线访问。随着时间的推移，它变得可以用相同的方式与其他编程语言（如 R）进行交互。由于这种从仅限 Python 的分离，这个工具发展成了当前的 Jupyter。IPython 仍然是一个可以使用的活跃工具。

Jupyter 作为一种 Web 应用程序，适用于各种平台。它还可以在桌面/笔记本电脑上通过多种安装方式使用。在本书中，我们将探索在 Windows PC 上使用 Jupyter，并通过互联网访问其他提供商。

Jupyter 概念

Jupyter 基于几个基本概念组织：

Notebook：一组语句（使用某种语言）。例如，这可以是一个完整的 R 脚本，用于加载数据、分析数据、生成图形，并将结果记录到其他地方。
单元格：Jupyter Notebook 中可以操作的最小粒度的部分：
- 当前单元格：当前正在编辑的单元格或已选择的单元格
内核：每个 notebook 都与特定的语言实现相关联。Jupyter 中处理特定语言的部分称为内核。

初步了解 Jupyter 用户界面

我们可以直接跳进去，看看 Jupyter 能提供什么。一个 Jupyter 屏幕看起来是这样的：

所以，Jupyter 被部署为一个网站，可以在你的机器上访问（也可以像访问任何其他网站一样通过互联网访问）。

我们可以看到页面的 URL，http://localhost:8888/tree。localhost 是指在你机器上运行的 Web 服务器的别名。我们在 Web 服务器上访问的网站以树状显示。这是默认的显示方式。它符合 Jupyter 中项目的显示方式。Jupyter 以树状布局显示对象，就像 Windows 文件资源管理器一样。主页列出了多个项目；每个项目都是一个子目录，并包含进一步的内容划分。根据你启动 Jupyter 的位置，当前目录中的现有内容也会显示在界面中。

详细描述 Jupyter 标签页

在网页上，我们可以看到即将熟悉的 Jupyter 标志和三个标签：

文件
运行
集群

文件标签列出了 Jupyter 可用的对象。Jupyter 使用的文件作为常规文件存储在磁盘上。Jupyter 提供了上下文管理器，能够处理你使用的不同类型的文件和程序。你可以通过 Windows 文件资源管理器查看 Jupyter 文件（它们的文件扩展名为 .ipynb）。你也可以在 Jupyter 窗口中看到非 Jupyter 文件。

运行标签列出了已经启动的 notebooks。Jupyter 会跟踪哪些 notebooks 正在运行。这个标签页允许你控制哪些 notebooks 在任何时候都在运行。

集群标签用于多个机器共同运行 Jupyter 的环境。

Jupyter 的集群实现是一个值得独立讨论的话题。

我可以在 Jupyter 中执行哪些操作？

接下来，我们看到：

提示选择要执行操作的项目
上传按钮
新的下拉菜单和
刷新图标

提示告诉你可以选择多个项目，然后对它们执行相同的操作。大多数以下操作（在菜单中）可以对单个项目或选定的一组项目执行。

上传按钮将提示选择一个文件并上传到 Jupyter。这通常用于将数据文件移动到项目中，以便在 Jupyter 作为网站运行在远程位置时访问，在这种情况下，你不能直接将文件复制到 Jupyter 运行的磁盘中。

新的下拉菜单展示了可用的不同类型的 Jupyter 项目（内核）的选择列表：

我们可以看到 Jupyter 能够创建的对象列表：

文本文件：创建一个用于该文件夹的文本文件。例如，如果笔记本要导入一个文件，你可以通过此功能创建该文件。
文件夹：是的，就像 Windows 文件资源管理器一样。
终端不可用：灰显，该功能只能在 Nix 环境中使用。
笔记本：灰显——这实际上不是一个文件类型，而是指该安装可以创建的不同类型的笔记本的标题。
Julia 0.4.5：创建一个使用 Julia 语言编写代码的 Julia 笔记本。
Python 3：创建一个使用 Python 语言编写代码的笔记本。这是默认选项。
R：创建一个使用 R 语言编写代码的笔记本。
根据你在安装中安装的内核，你可能会看到其他类型的笔记本列出。

Jupyter 可以操作哪些对象？

如果我们启动了其中一个笔记本（它会自动在 Jupyter 对象列表中被选中），然后查看针对所选对象的操作下拉菜单，我们会看到如下显示：

我们看到菜单操作已更改为重命名，因为这通常是对一个文件进行的操作，我们还可以看到一个图标用来删除项目（垃圾桶图标）。

项目计数现在是 1（我们在列表中选中了一个对象），该项目的图标是一个填充的蓝色方块（表示它是一个正在运行的项目），以及一个熟悉的主页图标，带我们返回到前面截图中的 Jupyter 主页显示。

对象菜单中有以下选择项：

文件夹：选择可用的文件夹
所有笔记本：选择 Jupyter 笔记本
正在运行：选择正在运行的 Jupyter 笔记本
文件：选择目录中的文件

如果我们在对象显示区域向下滚动，我们会看到在可用对象的列表中显示了一些不同的信息。每个列出的对象都有一个类型（由相关的图标形状表示）和一个用户在创建时分配的名称。

每个对象都是一个 Jupyter 项目，可以单独访问、共享和移动。每个项目都有一个完整的名称（由创建该项目的用户输入），并且有一个图标表示该条目是一个项目。我们还会看到其他 Jupyter 图标，代表其他项目组件，具体如下：

查看 Jupyter 项目显示

如果我们下拉“新建”菜单并选择 Python 3，Jupyter 将创建一个新的 Python 笔记本并展示其内容。我们将看到如下所示的界面：

我们已经创建了一个新的 Jupyter 笔记本，并且进入了它的显示界面。Logo 已经显示，标题默认是Untitled，我们可以点击它进行修改。还有一个（自动保存）标记，告诉你 Jupyter 已将笔记本自动保存到磁盘，并且在你继续工作时会定期自动保存。

我们现在有了一个菜单栏，并且标明了该笔记本使用的是 Python 3 作为源语言。菜单选项包括：

文件：标准文件操作
编辑：用于编辑单元格内容（更多内容待更新）
视图：改变笔记本的显示方式
插入：在笔记本中插入一个单元格
单元格：改变单元格的格式和使用方式
内核：调整笔记本使用的内核
帮助：调出 Jupyter 的帮助系统

文件菜单

文件菜单包含以下选项：

新建笔记本：类似于主页上的下拉操作
打开...：打开一个笔记本
复制一份...：复制一个笔记本
重命名...：重命名一个笔记本
保存和检查点：在检查点保存当前笔记本。检查点是笔记本历史中的特定时刻，用于在你对最近的更改感到后悔时，能够返回到某个检查点。
打印预览：类似于你曾用过的其他打印预览功能。
另存为：允许你将笔记本存储为各种格式。最常见的格式包括 PDF 或 Excel，这样可以将笔记本分享给那些无法访问 Jupyter 的用户。
可信笔记本：（此功能是灰显的）。当一个用户打开一个笔记本时，服务器会用用户的密钥计算一个签名，并将其与存储在笔记本元数据中的签名进行比较。如果签名匹配，那么笔记本中的 HTML 和 JavaScript 输出将在加载时被信任，否则将被视为不可信。
关闭并停止：关闭当前笔记本并停止它在 Jupyter 系统中的运行

编辑菜单

编辑菜单包含以下选项：

剪切单元格：典型的剪切操作。
复制单元格：假设你已经习惯了将单元格复制到内存缓冲区并稍后粘贴到笔记本的其他位置的图形界面操作。
粘贴上方单元格：如果您已选择一个单元格并且已经复制了一个单元格，则此选项不会被灰显，并且会将缓冲区中的单元格粘贴到当前单元格上方。
粘贴下方单元格：与之前的选项类似。
删除单元格：将删除所选的单元格。
撤销删除单元格。
拆分单元格：这里存在样式问题，关于将多少语句放入一个单元格。许多时候，您会从一个包含多个语句的单元格开始，然后将该单元格多次拆分，将单个语句或语句组分割到各自的单元格中。
合并上方单元格：将当前单元格与上方的单元格合并。
合并下方单元格：与之前的选项类似。
向上移动单元格：将当前单元格移到上方单元格之前。
向下移动单元格。
编辑笔记本元数据：供高级用户修改 Jupyter 为您的笔记本使用的内部编程语言。
查找与替换：在单元格中定位特定文本并可能进行替换。

查看菜单

查看菜单有以下选项：

切换标题：切换 Jupyter 标题的显示
切换工具栏：切换 Jupyter 工具栏的显示
单元格工具栏：更改正在编辑的单元格的显示项：
- 无：不显示单元格工具栏
- 编辑元数据：直接编辑单元格的元数据
- 原始单元格格式：编辑 Jupyter 使用的单元格原始格式
- 幻灯片：以幻灯片方式浏览单元格

插入菜单

插入菜单有以下选项：

插入上方单元格：将复制的缓冲区单元格插入到当前单元格前面
插入下方单元格：与之前的选项相同

单元格菜单

单元格菜单有以下选项：

运行单元格：运行笔记本中的所有单元格
运行单元格并选择下方：运行单元格并选择当前单元格下方的所有单元格
运行单元格并插入下方：运行单元格并添加一个空白单元格
运行所有单元格：运行笔记本中的所有单元格
运行上方所有单元格：运行当前单元格上方的所有单元格
运行下方所有单元格：运行当前单元格下方的所有单元格
单元格类型：将选定的单元格类型更改为：
- 代码：这是默认的——该单元格期望包含语言语句
- Markdown：单元格包含 HTML markdown，通常用于以最佳方式显示笔记本（因为它是一个网站，因此可以使用所有 HTML）
- 原始 NBConvert：这是 Jupyter 的一种内部格式，基本上是纯文本
当前输出：是否清除或继续单元格的输出
所有输出

Kernel 菜单

Kernel 菜单用于控制笔记本使用的底层语言引擎。菜单选项如下。我认为这个菜单中的许多选项使用得非常少：

中断：暂时停止底层语言引擎，然后让它继续
重启：重启底层语言引擎
重启并清除输出
重启并运行全部
重新连接：如果你中断了内核，那么你需要重新连接才能继续运行
更改内核：将本笔记本的语言更改为你安装中可用的语言

帮助菜单

帮助菜单显示 Jupyter 的帮助选项和语言上下文选择。例如，在我们的 Python 笔记本中，我们可以看到用于常用 Python 库的选择：

图标工具栏

正常菜单下方是一个图标工具栏，包含许多常用菜单项，方便快速使用，如图所示：

这些图标对应前面菜单中的选项（按出现顺序列出）：

文件/保存当前笔记本
插入下方单元格
剪切当前单元格
复制当前单元格
将单元格粘贴到下方
向上移动选中的单元格
向下移动选中的单元格
从选中的单元格开始运行
中断内核
重启内核
我们可以应用于当前单元格的格式列表
一个图标用于打开一个带有描述性名称的命令面板
一个图标用于打开单元格工具栏

执行脚本时效果如何？

如果我们为笔记本命名，输入一个简单的 Python 脚本，并执行笔记本单元格，我们将看到如下显示：

脚本是：

name = "Dan Toomey"
state = "MA"
print(name + " lives in " + state)

我们为变量name和state赋值，然后打印它们。

如果你注意到，我将语句放入了两个不同的单元格。这仅仅是为了提高可读性。它们本可以放在同一个单元格或三个不同的单元格中。

每个单元格都有行号。编号总是从第一个单元格开始为 1，然后随着单元格的移动，编号可能会增加（正如你所看到的，显示中的第一个单元格标记为单元格 2）。

在第二个单元格下方，我们有不可编辑的显示结果。Jupyter 总是将单元格的任何对应输出显示在其下方。这可能还包括错误信息。

行业数据科学应用

本书讲解 Jupyter 和数据科学。我们已经介绍了 Jupyter。现在，我们可以看看数据科学的实践，再看看这两个概念是如何协同工作的。

数据科学在许多行业中得到应用。有趣的是，行业中涉及的主要技术和使用的算法。我们可以看到 Jupyter 中可用的相同技术。

一些数据科学的大型用户行业包括：

行业	更大的数据科学应用	技术/算法
金融	对冲基金	Python
赌博	确定赔率	R
保险	风险衡量与定价	Domino (R)
零售银行	风险、客户分析、产品分析	R
采矿	智能勘探、产量优化	Python
消费品	定价与分销	R
医疗保健	药物发现与试验	Python

所有这些数据科学分析可以在 Jupyter 中完成，因为所使用的语言都得到了完全支持。

现实生活中的例子

在这一部分，我们展示了从当前行业焦点中提取的几个例子，并将它们应用于 Jupyter，以确保其实用性。

金融，Python - 欧式看涨期权估值

在www.safaribooksonline.com/library/view/python-for-finance/9781491945360/ch03.html中有一个例子，来自 Yves Hilpisch 的书《Python for Finance》。所使用的模型在金融工作中非常标准。

我们希望得出看涨期权的理论价值。看涨期权是指在特定时间内以特定（执行）价格购买证券（如 IBM 股票）的权利。期权定价基于证券相对于执行价格和当前价格的风险性或波动性。这个例子使用的是欧式期权，只有在到期时才能行使——这简化了问题集。

这个例子使用了 Black-Scholes 模型来估值期权，其中我们有：

初始股票指数水平 S[0] = 100
欧式看涨期权的执行价格 K = 105
到期时间 T = 1 年
常数，无风险短期利率 r = 5%
常数波动率 σ = 20%

这些元素组成了以下公式：

使用的算法如下：

从标准正态分布中绘制 I（伪）随机数。
计算给定 z(i) 下的所有期权到期时的指数水平 S[T](i)。计算期权到期时的所有内在价值 h[T](i) = max(S[T](i) - K,0)。
通过蒙特卡洛估算器估算期权现值，公式如下：

脚本如下。我们使用 numpy 进行复杂的数学计算。其余的代码是典型的：

from numpy import * # set parameters S0 = 100. K = 105. T = 1.0 r = 0.05 sigma = 0.2 # how many samples we are using I = 100000 random.seed(103) z = random.standard_normal(I) ST = S0 * exp((r - 0.5 * sigma ** 2) * T + sigma * sqrt(T) * z) hT = maximum(ST - K, 0) C0 = exp(-r * T) * sum(hT) / I # tell user results print ("Value of the European Call Option %5.3f" % C0)

Jupyter 下的结果如以下截图所示：

8.071 的值与已发布的期望值 8.019 一致，原因是使用的随机数有方差。（我在给随机数生成器设置种子，以确保结果可复现）。

金融，Python - 蒙特卡洛定价

另一个常用的算法是蒙特卡罗模拟。正如赌博度假村的名字所暗示的，在蒙特卡罗模拟中，我们模拟了一种情境，其中我们知道不同结果的百分比，但不确定接下来 N 次机会会发生什么。我们可以在www.codeandfinance.com/pricing-options-monte-carlo.html看到这个模型的应用。在这个例子中，我们再次使用了 Black-Scholes 模型，但采用了另一种直接方法，我们可以看到每个步骤。

编码如下。Jupyter 中的 Python 编码风格与直接在 Python 中使用的略有不同，您可以从代码顶部附近的更改导入看到。与仅从库中提取所需函数不同，您会导入整个库，然后使用所需的编码：

import datetime import random # import gauss import math #import exp, sqrt random.seed(103) def generate_asset_price(S,v,r,T):
 return S * exp((r - 0.5 * v**2) * T + v * sqrt(T) * gauss(0,1.0)) def call_payoff(S_T,K):
 return max(0.0,S_T-K) S = 857.29 # underlying price v = 0.2076 # vol of 20.76% r = 0.0014 # rate of 0.14% T = (datetime.date(2013,9,21) - datetime.date(2013,9,3)).days / 365.0 K = 860. simulations = 90000 payoffs = [] discount_factor = math.exp(-r * T) for i in xrange(simulations):
 S_T = generate_asset_price(S,v,r,T) payoffs.append( call_payoff(S_T, K) ) price = discount_factor * (sum(payoffs) / float(simulations)) print ('Price: %.4f' % price)

在 Jupyter 下的结果如下所示：

结果价格为14.4452，接近发布的值14.5069。

赌博，R - 投注分析

一些赌博游戏实际上是硬币翻转，成功的几率是 50/50。沿着这些线路，我们有从forumserver.twoplustwo.com/25/probability/flipping-coins-getting-3-row-1233506/获取的编码，用于确定硬币翻转中一系列正面或反面的概率，如果您知道硬币/游戏偏向于一个结果，可以使用触发器。

我们有以下脚本：

##############################################
# Biased/unbiased  recursion of heads OR tails
##############################################
import numpy as np
import math

N = 14     # number of flips
m = 3      # length of run (must be  > 1 and <= N/2)
p = 0.5   # P(heads)

prob = np.repeat(0.0,N)
h = np.repeat(0.0,N)
t = np.repeat(0.0,N)

h[m] = math.pow(p,m)
t[m] = math.pow(1-p,m)
prob[m] = h[m] + t[m]

for n in range(m+1,2*m):
 h[n] = (1-p)*math.pow(p,m)
 t[n] = p*math.pow(1-p,m)
 prob[n] = prob[n-1] + h[n] + t[n]

for n in range(2*m,N):
 h[n] = ((1-p) - t[n-m] - prob[n-m-1]*(1-p))*math.pow(p,m)
 t[n] = (p - h[n-m] - prob[n-m-1]*p)*math.pow(1-p,m)
 prob[n] = prob[n-1] + h[n] + t[n]

prob[N-1]

前面的代码在 Jupyter 中生成以下输出：

我们最终得到在一个公平的游戏中连续三次抛硬币的概率。在这种情况下，有 92%的机会（在我们运行的 14 次抛硬币的范围内）。

保险，R - 非人寿保险定价

我们有一个使用 R 为非寿险产品，特别是脚踏车，制定定价的例子，网址为www.cybaea.net/journal/2012/03/13/R-code-for-Chapter-2-of-Non_Life-Insurance-Pricing-with-GLM/。该代码首先创建产品线的统计数据表，然后将定价与实际使用的统计数据进行比较。

积累数据的代码的第一部分如下所示：

con <- url("http://www2.math.su.se/~esbj/GLMbook/moppe.sas") data <- readLines(con, n = 200L, warn = FALSE, encoding = "unknown") close(con) ## Find the data range data.start <- grep("^cards;", data) + 1L data.end   <- grep("^;", data[data.start:999L]) + data.start - 2L table.1.2  <- read.table(text = data[data.start:data.end],
 header = FALSE, sep = "", quote = "", col.names = c("premiekl", "moptva", "zon", "dur",
 "medskad", "antskad", "riskpre", "helpre", "cell"), na.strings = NULL, colClasses = c(rep("factor", 3), "numeric", rep("integer", 4), "NULL"), comment.char = "") rm(con, data, data.start, data.end) # Remainder of Script adds comments/descriptions comment(table.1.2) <-
 c("Title: Partial casco moped insurance from Wasa insurance, 1994--1999", "Source: http://www2.math.su.se/~esbj/GLMbook/moppe.sas", "Copyright: http://www2.math.su.se/~esbj/GLMbook/") ## See the SAS code for this derived field table.1.2$skadfre = with(table.1.2, antskad / dur) ## English language column names as comments: comment(table.1.2$premiekl) <-
 c("Name: Class", "Code: 1=Weight over 60kg and more than 2 gears", "Code: 2=Other") comment(table.1.2$moptva)   <-
 c("Name: Age", "Code: 1=At most 1 year", "Code: 2=2 years or more") comment(table.1.2$zon)      <-
 c("Name: Zone", "Code: 1=Central and semi-central parts of Sweden's three largest cities", "Code: 2=suburbs and middle-sized towns", "Code: 3=Lesser towns, except those in 5 or 7", "Code: 4=Small towns and countryside, except 5--7", "Code: 5=Northern towns", "Code: 6=Northern countryside", "Code: 7=Gotland (Sweden's largest island)") comment(table.1.2$dur)      <-
 c("Name: Duration", "Unit: year") comment(table.1.2$medskad)  <-
 c("Name: Claim severity", "Unit: SEK") comment(table.1.2$antskad)  <- "Name: No. claims" comment(table.1.2$riskpre)  <-
 c("Name: Pure premium", "Unit: SEK") comment(table.1.2$helpre)   <-
 c("Name: Actual premium", "Note: The premium for one year according to the tariff in force 1999", "Unit: SEK") comment(table.1.2$skadfre)  <-
 c("Name: Claim frequency", "Unit: /year") ## Save results for later save(table.1.2, file = "table.1.2.RData") ## Print the table (not as pretty as the book) print(table.1.2)

表的前 10 行结果如下：

       premiekl moptva zon    dur medskad antskad riskpre helpre    skadfre
    1         1      1   1   62.9   18256      17    4936   2049 0.27027027
    2         1      1   2  112.9   13632       7     845   1230 0.06200177
    3         1      1   3  133.1   20877       9    1411    762 0.06761833
    4         1      1   4  376.6   13045       7     242    396 0.01858736
    5         1      1   5    9.4       0       0       0    990 0.00000000
    6         1      1   6   70.8   15000       1     212    594 0.01412429
    7         1      1   7    4.4    8018       1    1829    396 0.22727273
    8         1      2   1  352.1    8232      52    1216   1229 0.14768532
    9         1      2   2  840.1    7418      69     609    738 0.08213308
    10        1      2   3 1378.3    7318      75     398    457 0.05441486

然后，我们逐个产品/统计数据确定产品的定价是否与其他产品一致。请注意，install.packages语句中的repos =子句是 R 的一个相当新的添加：

# make sure the packages we want to use are installed install.packages(c("data.table", "foreach", "ggplot2"), dependencies = TRUE, repos = "http://cran.us.r-project.org") # load the data table we need if (!exists("table.1.2"))
 load("table.1.2.RData") library("foreach") ## We are looking to reproduce table 2.7 which we start building here, ## add columns for our results. table27 <-
 data.frame(rating.factor = c(rep("Vehicle class", nlevels(table.1.2$premiekl)), rep("Vehicle age",   nlevels(table.1.2$moptva)), rep("Zone",          nlevels(table.1.2$zon))), class = c(levels(table.1.2$premiekl), levels(table.1.2$moptva), levels(table.1.2$zon)), stringsAsFactors = FALSE) ## Calculate duration per rating factor level and also set the ## contrasts (using the same idiom as in the code for the previous ## chapter). We use foreach here to execute the loop both for its ## side-effect (setting the contrasts) and to accumulate the sums. # new.cols are set to claims, sums, levels new.cols <-
 foreach (rating.factor = c("premiekl", "moptva", "zon"), .combine = rbind) %do% {
 nclaims <- tapply(table.1.2$antskad, table.1.2[[rating.factor]], sum) sums <- tapply(table.1.2$dur, table.1.2[[rating.factor]], sum) n.levels <- nlevels(table.1.2[[rating.factor]]) contrasts(table.1.2[[rating.factor]]) <- contr.treatment(n.levels)[rank(-sums, ties.method = "first"), ] data.frame(duration = sums, n.claims = nclaims) } table27 <- cbind(table27, new.cols) rm(new.cols) #build frequency distribution model.frequency <-
 glm(antskad ~ premiekl + moptva + zon + offset(log(dur)), data = table.1.2, family = poisson) rels <- coef( model.frequency ) rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] ) table27$rels.frequency <-
 c(c(1, rels[1])[rank(-table27$duration[1:2], ties.method = "first")], c(1, rels[2])[rank(-table27$duration[3:4], ties.method = "first")], c(1, rels[3:8])[rank(-table27$duration[5:11], ties.method = "first")]) # note the severities involved model.severity <-
 glm(medskad ~ premiekl + moptva + zon, data = table.1.2[table.1.2$medskad > 0, ], family = Gamma("log"), weights = antskad) rels <- coef( model.severity ) rels <- exp( rels[1] + rels[-1] ) / exp( rels[1] ) ## Aside: For the canonical link function use ## rels <- rels[1] / (rels[1] + rels[-1]) table27$rels.severity <-
 c(c(1, rels[1])[rank(-table27$duration[1:2], ties.method = "first")], c(1, rels[2])[rank(-table27$duration[3:4], ties.method = "first")], c(1, rels[3:8])[rank(-table27$duration[5:11], ties.method = "first")]) table27$rels.pure.premium <- with(table27, rels.frequency * rels.severity) print(table27, digits = 2)

结果显示如下：

       rating.factor class duration n.claims rels.frequency rels.severity
    1  Vehicle class     1     9833      391           1.00          1.00
    2  Vehicle class     2     8825      395           0.78          0.55
    11   Vehicle age     1     1918      141           1.55          1.79
    21   Vehicle age     2    16740      645           1.00          1.00
    12          Zone     1     1451      206           7.10          1.21
    22          Zone     2     2486      209           4.17          1.07
    3           Zone     3     2889      132           2.23          1.07
    4           Zone     4    10069      207           1.00          1.00
    5           Zone     5      246        6           1.20          1.21
    6           Zone     6     1369       23           0.79          0.98
    7           Zone     7      148        3           1.00          1.20
       rels.pure.premium
    1               1.00
    2               0.42
    11              2.78
    21              1.00
    12              8.62
    22              4.48
    3               2.38
    4               1.00
    5               1.46
    6               0.78
    7               1.20

在这里，我们可以看到一些车辆类别（2，6）的价格非常低，与该车辆的统计数据相比，而其他车辆则被高估（12，22）。

消费品，R - 营销效果

我们以一份我在www.dantoomeysoftware.com/Using_R_for_Marketing_Research.pptx上做的演示为例，分析不同广告系列对葡萄果汁的效果。

代码如下：

#library(s20x)
library(car)

#read the dataset from an existing .csv file
df <- read.csv("C:/Users/Dan/grapeJuice.csv",header=T)

#list the name of each variable (data column) and the first six rows of the dataset
head(df)

# basic statistics of the variables
summary(df)

#set the 1 by 2 layout plot window
par(mfrow = c(1,2))

# boxplot to check if there are outliers
boxplot(df$sales,horizontal = TRUE, xlab="sales")

# histogram to explore the data distribution shape
hist(df$sales,main="",xlab="sales",prob=T)
lines(density(df$sales),lty="dashed",lwd=2.5,col="red")

#divide the dataset into two sub dataset by ad_type
sales_ad_nature = subset(df,ad_type==0)
sales_ad_family = subset(df,ad_type==1)

#calculate the mean of sales with different ad_type
mean(sales_ad_nature$sales)
mean(sales_ad_family$sales)

#set the 1 by 2 layout plot window
par(mfrow = c(1,2))

# histogram to explore the data distribution shapes
hist(sales_ad_nature$sales,main="",xlab="sales with nature production theme ad",prob=T)
lines(density(sales_ad_nature$sales),lty="dashed",lwd=2.5,col="red")

hist(sales_ad_family$sales,main="",xlab="sales with family health caring theme ad",prob=T)
lines(density(sales_ad_family$sales),lty="dashed",lwd=2.5,col="red")

带有输出（多个部分）：

（来自文件的原始数据，前 10 行）：

	销售	价格	广告类型	苹果价格	饼干价格
1	222	9.83	0	7.36	8.8
2	201	9.72	1	7.43	9.62
3	247	10.15	1	7.66	8.9
4	169	10.04	0	7.57	10.26
5	317	8.38	1	7.33	9.54
6	227	9.74	0	7.51	9.49

数据的统计如下：

    sales           price           ad_type     price_apple 
    Min.   :131.0   Min.   : 8.200   Min.   :0.0   Min.   :7.300 
    1st Qu.:182.5   1st Qu.: 9.585   1st Qu.:0.0   1st Qu.:7.438 
    Median :204.5   Median : 9.855   Median :0.5   Median :7.580 
    Mean   :216.7   Mean   : 9.738   Mean   :0.5   Mean   :7.659 
    3rd Qu.:244.2   3rd Qu.:10.268   3rd Qu.:1.0   3rd Qu.:7.805 
    Max.   :335.0   Max.   :10.490   Max.   :1.0   Max.   :8.290 
    price_cookies 
    Min.   : 8.790 
    1st Qu.: 9.190 
    Median : 9.515 
    Mean   : 9.622 
    3rd Qu.:10.140 
    Max.   :10.580

数据显示了每个活动的效果。家庭销售更为有效：

186.666666666667//自然销售平均值
246.666666666667//家庭销售平均值

在直方图显示中，差异更加明显：

在 Jupyter 中使用 Docker

Docker 是一种机制，允许您在一台机器上拥有许多完整的应用程序虚拟实例。许多软件公司使用 Docker 来提供其服务的完全可扩展实现，并支持尽可能多的并发用户。

以前处理多个实例的机制共享公共资源（如磁盘地址空间）。在 Docker 下，每个实例都是一个与其他实例完全独立的完整实体。

在 Docker 环境中实现 Jupyter 允许多个用户访问各自的 Jupyter 实例，而无需担心干扰他人的计算。

Docker 的关键特性是允许在任何时候使用多个实例的笔记本。Docker 控制系统可以设置为为每个访问您笔记本的用户创建新的实例。所有这些都内建在 Docker 中，无需编程；只需使用用户界面来决定如何创建实例。

您可以通过两种方式使用 Docker：

从公共服务
在您的机器上安装 Docker

使用公共 Docker 服务

有几个服务可以使用。我认为它们的工作方式基本相同：注册服务，上传您的笔记本，监控使用情况（Docker 控制程序会自动跟踪使用情况）。例如，如果我们使用 hub.docker.com/，我们实际上是在使用一个笔记本的版本库。版本控制在软件开发中用于跟踪随时间做出的更改。这也允许多个用户的访问权限：

首先，注册。这将为您提供对服务提供商的身份验证。
创建一个仓库——在这里您将保留您的笔记本版本。
您需要在您的机器上安装 Docker，才能从仓库中拉取/推送笔记本。

安装 Docker 取决于操作系统。请访问 www.docker.com/ 首页，按照您机器的说明进行操作。

将您的 Jupyter 镜像上传（推送）到仓库中。
访问仓库中的笔记本。您可以在 Docker 控制下与他人共享笔记本的地址（URL），并为不同用户设置特定的访问权限。
从那时起，它将像在本地运行一样工作。

在您的机器上安装 Docker

在您的本地机器上安装 Docker 只会是发布到公共 Docker 服务的前提，除非您安装 Docker 的机器可以被他人访问。

另一个选择是将 Docker 安装在您的机器上。它与之前的情况完全相同，只是您正在管理 Docker 镜像空间。

如何与他人共享笔记本

有几种方法可以与他人共享 Jupyter 笔记本：

电子邮件
上传到 Google Drive
在 GitHub 上共享
将其存储为 HTML 文件在网络服务器上
在网络服务器上安装 Jupyter

你能通过电子邮件发送一个笔记本吗？

为了通过电子邮件发送笔记本，必须将笔记本转换为纯文本格式，作为附件发送给接收者，然后接收者必须将其转换回“二进制”笔记本格式。

邮件附件通常会被转换为明确的MIME（多用途互联网邮件扩展）格式。有一个程序可以将笔记本格式转换为 MIME 格式，名为nb2mail，它可以将笔记本转换为笔记本 MIME 格式。该程序可以在github.com/nfultz/nb2mail找到。

使用方法如下：

使用pip命令安装nb2mail（见网站）
将选定的笔记本转换为 MIME 格式
发送给接收者
接收者的 MIME 转换过程将以正确的方式存储文件（假设他们也安装了nb2mail）

在 Google Drive 上共享笔记本

Google Drive 可用于存储你的笔记本档案信息。当与之前的笔记本邮件发送功能结合使用时，接收者可以使用 Google Drive 档案，这样没有档案信息的人就无法与笔记本互动。

你可以使用pip安装 python 扩展（来自github.com/jupyter/jupyter-drive），然后使用python -m。从此以后，你可以通过 Google Drive 档案来访问笔记本，命令为ipython notebook -profile <profilename>。

在 GitHub 上共享

GitHub（以及其他平台）允许你将笔记本放在他们的服务器上，一旦放置，便可以通过 nbviewer 直接访问。服务器已安装支持你的笔记本所需的 Python（及其他语言）代码。nbviewer 是笔记本的只读版本，不支持互动。

nbviewer 可以在github.com/jupyter/nbviewer找到。该站点包括需要添加到ipython notebook命令的特定参数，例如启动查看器的命令。

将其存储为 HTML 文件在 Web 服务器上

笔记本的内置功能是将笔记本导出为不同格式，其中之一是 HTML。通过这种方式，你可以将笔记本导出为 HTML，并在做出更改时将文件复制到你的 Web 服务器上。

命令为jupyter nbconvert <notebook name>.ipynb --to html。

这将是一个非互动的只读版本的笔记本。

在网络服务器上安装 Jupyter

Jupyter 作为 Web 应用程序部署。如果你有直接访问 Web 服务器的权限，你可以在 Web 服务器上安装 Jupyter，创建笔记本，然后这些笔记本将可供其他人访问，且完全动态。

作为一个 Web 服务器，你也可以控制对 Web 服务器的访问，因此能够控制谁可以访问你的笔记本。

这是一个高级交互，需要与网站管理员合作，确定正确的处理方式。

如何确保笔记本的安全？

Jupyter 笔记本的安全性有两个方面：

确保只有特定用户能够访问你的笔记本
确保你的笔记本不会用于托管恶意代码

访问控制

尽管 Jupyter 的许多使用场景主要是用于教育他人，但也有一些情况下，所访问的信息是机密的，且应该保持机密。Jupyter 允许你通过多种方式为你的笔记本设置访问限制。

当我们识别用户时，我们实际上是在验证该用户的身份。通常，通过在允许进入之前提出登录挑战来完成此过程，用户需要输入用户名和密码。

如果 Jupyter 实例托管在 Web 服务器上，你可以使用 Web 服务器的访问控制来限制对笔记本的访问。此外，大多数支持笔记本托管的供应商提供了一种机制，允许限制特定用户的访问。

恶意内容

安全性的另一个方面是确保你的笔记本内容不包含恶意代码。你应该确保你的笔记本是安全的，如下所示：

确保 HTML 被清理（检查恶意 HTML 代码并消除它）
不要允许笔记本执行外部 JavaScript
检查可能包含恶意内容的单元格是否在服务器环境中受到挑战
清理单元格输出，以避免对用户机器产生不必要的影响

总结

在本章中，我们深入探讨了 Jupyter 用户界面的细节：它与哪些对象进行交互，Jupyter 可以执行哪些操作，显示内容告诉我们数据的哪些信息，以及有哪些工具可供使用？接下来，我们查看了来自多个行业的 R 和 Python 编码的实际案例。然后，我们了解了将笔记本与其他用户共享的几种方式，以及如何通过不同的安全机制保护我们的笔记本。

在下一章中，我们将看到如何在 Jupyter 笔记本中使用 Python。

第二章：在 Jupyter 上处理分析数据

Jupyter 并不做分析数据的繁重工作：所有工作都是由用所选编程语言编写的程序完成的。Jupyter 提供了一个框架，可以运行多种编程语言模块。因此，我们可以选择如何在 Jupyter 中分析数据。

Python 是数据分析编程的流行选择。Jupyter 完全支持 Python 编程。我们将查看一些可能对这种支持系统造成压力的编程解决方案，并看看 Jupyter 如何应对。

使用 Python 笔记本进行数据抓取

数据分析的常用工具是从公共源（如网站）收集数据。Python 擅长从网站抓取数据。在这里，我们看一个示例，加载来自 Google 财经的数据的股票价格信息。

特别地，给定一个股票符号，我们想要获取该符号的过去一年的价格范围。

Google 财经网站上的一页将提供一只证券公司的过去一年的价格数据。例如，如果我们对 超威半导体（AMD）的价格点感兴趣，我们将输入以下 URL：

www.google.com/finance/historical?q=NASDAQ:AMD

在这里，NASDAQ 是承载 AMD 股票的证券交易所。在生成的 Google 页面上，有一个包含感兴趣数据点的表格，如下图所示的部分截图所示。

就像你将尝试访问的许多网站一样，页面上还有很多其他信息，比如头部、尾部和广告，正如你在以下截图中看到的那样。网页是为人类读者构建的。幸运的是，Google 和其他公司意识到你在抓取他们的数据，并且会保持数据的相同格式，这样你就不必修改脚本。

提前警告你，如果你过于频繁地访问某个页面或整个网站，可能会被阻止访问。访问频率是你访问的特定网站需要讨论的问题。再次强调，这些网站知道你在抓取数据，并且只要不干扰正常的人工网络流量，他们是可以接受的。

在该网页上有一个清晰的表格。如果我们查看生成网页所用的底层 HTML，我们会发现有很多头部、尾部和侧边栏的信息，但更重要的是，我们会发现一个 id 为 price_data 的 HTML div 标签。在这个 div 标签中，我们看到一个 HTML 表格，每一行包含该数据的 date、opening price、high、low、close 和 volume 的值，正如屏幕上所看到的那样。

我们可以使用一个标准的 Python 库包，lxml，将网页文本加载并解析为 HTML Python 组件，以便我们进行操作。

然后，对于每一天的数据，我们提取出列信息，并将其添加到我们的数据列表中。

通常，你可能每天运行此脚本一次，并将最新一天的信息存储到本地数据库中以供进一步分析。在我们的案例中，我们只是将最后一天的值打印在屏幕上。

使用的 Python 脚本如下：

from lxml import html import requests from time import sleep # setup the URL for the symbol we are interested in exchange = "NASDAQ" ticker = "AMD" url = "https://www.google.com/finance/historical?q=%s:%s"%(exchange,ticker) # retrieve the web page response = requests.get(url) print ("Retrieving prices for %s from %s"%(ticker,url)) # give it a few seconds in case there is some delay sleep(3) # convert the text into an HTML Document parser = html.fromstring(response.text) # find the HTML DIV tag that has id 'prices' price_store = parser.get_element_by_id("prices") # we will store our price information in the price_data list price_data = [] # find the HTML TABLE element within the prices DIV for table in price_store:
 #every row (skip first row headings) of table has #  date, open, high, low, close, volume for row in table[1:]: #store tuples for a day together day = {"date":row[0].text.strip('\n'), \ "open":row[1].text.strip('\n'), \ "high":row[2].text.strip('\n'), \ "low":row[3].text.strip('\n'), \ "close":row[4].text.strip('\n'), \ "volume":row[5].text.strip('\n')}        #add day's information to our set
 price_data.append(day) print ("The last day of pricing information we have is:") print (price_data[0])

在 Jupyter 控制台运行此脚本，我们看到如下部分截图中的结果：

在 Jupyter 中使用重型数据处理函数

Python 有多个处理函数组，可能会占用计算机系统的计算资源。让我们在 Jupyter 中使用其中的一些，来检查其功能是否按预期执行。

在 Jupyter 中使用 NumPy 函数

NumPy 是 Python 中的一个包，提供多维数组和数组处理的常用函数。我们通过import * from numpy语句引入 NumPy 包。特别是，NumPy 包定义了array关键字，表示一个具有广泛功能的 NumPy 对象。

NumPy 数组处理函数从基础的min()和max()函数（它们提供给定数组维度上的最小值和最大值）到更有趣的实用函数，用于生成直方图和计算数据框元素之间的相关性。

使用 NumPy，你可以以多种方式操作数组。例如，我们将通过以下脚本介绍其中一些函数，使用 NumPy 来：

创建一个数组
计算数组中的最大值
计算数组中的最小值
确定第二个轴的总和

# numpy arrays
import numpy as np

# create an array 'a' with 3 3-tuples
a = np.array([[1, 1, 2], [3, 5, 8], [13, 21, 34]])
print("Array contents", a)

# determine the minimum value in array
print("max value = ", a.max())

# max value in array
print("min value = ", a.min())

# sum across the 2nd axis  
print("sum across 2nd axis", a.sum(axis = 1))

如果我们将此脚本转移到 Python 笔记本中，当执行单元格时，我们会看到如下显示：

我们可以使用以下脚本，通过更有趣的histogram和correlate函数来处理数组：

import numpy as np
import random

# build up 2 sets of random numbers

# setup empty array 2 columns, 1000 rows
numbers = np.empty([2,1000], int)

# set seed so we can repeat results
random.seed(137)

# populate the array
for num in range(0, 1000):
 numbers[0,num] = random.randint(0, 1000)
 numbers[1,num] = random.randint(0, 1000)

# produce a histogram of the data
(hist, bins) = np.histogram(numbers, bins = 10, range = (0,1000))
print ("Histogram is ",hist)

# calculate correlation between the 2 columns

corrs = np.correlate(numbers[:,1], numbers[:,2], mode='valid')
print ("Correlation of the two rows is ", corrs)

在这个脚本中，我们正在：

用随机数填充一个两列的数组
在 100 个点的范围内生成两个列的直方图
最后，确定两个列之间的相关性（这应该是非常高的相关性）

将此脚本输入到 Jupyter Notebook 并执行单元格后，我们得到如下输出。可以理解的是，这些桶的大小非常接近：

在 Jupyter 中使用 pandas

pandas 是一个开源的高性能数据分析工具库，在 Python 中可用。特别感兴趣的是以下功能：

读取文本文件
读取 Excel 文件
从 SQL 数据库读取
操作数据框

使用 pandas 在 Jupyter 中读取文本文件

最常见的包含分析数据的文本文件类型是 CSV 文件。互联网上有许多以此格式提供的数据集。我们将查看在vincentarelbundock.github.io/Rdatasets/csv/datasets/Titanic.csv中找到的泰坦尼克号生还者数据。

与大多数 pandas 函数一样，函数调用非常简单：

import pandas as pd
df = pd.read_csv ('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/Titanic.csv')
print (df.head)

然而，和许多 pandas 函数一样，read_csv函数有一组广泛的可选参数，可以传入，默认值是最常用的功能，这样我们就可以像之前一样写出简短的代码来完成工作。我们可以使用的一些附加参数允许我们：

跳过行
跳过/定义列标题
并且可以更改索引字段（Python 总是希望在数据框中保持一个主要索引字段，以加快访问速度）

在 Jupyter 中执行的脚本结果显示如下（注意，我只使用head函数打印表格的前后 30 行）：

在 Jupyter 中使用 pandas 读取 Excel 文件

同样，我们也可以轻松加载 Microsoft Excel 文件。例如，关于同一泰坦尼克号数据集的 Excel 文件可以在vandebilt.edu找到（完整链接在以下脚本中）。我们有以下脚本：

import pandas as pd
df = pd.read_excel('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.xls')
print (df.head)

读取 Excel 文件时也有一组广泛的可选参数，例如：

选择 Excel 文件中的工作表进行读取
跳过行
指定 NA 值的处理方式

在 Jupyter 中显示的结果流程如下。数据集与之前读取的 CSV 文件非常相似。

使用 pandas 操作数据框

一旦我们有了数据框，就有几个 pandas 方法可以进一步处理数据。我们将使用 pandas 来：

groupby函数
操作列
计算异常值

在数据框中使用 groupby 函数

groupby函数可以用来根据您的标准对数据框中的记录进行分组（并计数）。

继续使用我们的泰坦尼克号数据集，我们可以使用groupby按年龄计数乘客数量。

我们可以使用以下脚本：

# read in the titanic data set
import pandas as pd
df = pd.read_excel('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.xls')
# extract just the age column to its own dataset, 
# group by the age, and
# add a count for each age
ages = df[['age']].groupby('age')['age'].count()
print (ages)

在 Jupyter 中显示的结果如下。我没意识到船上有这么多婴儿。

操作数据框中的列

一个有趣的列操作是排序。我们可以使用sort_values函数对先前的年龄计数数据进行排序，以确定船上旅客的最常见年龄。

脚本如下：

import pandas as pd
df = pd.read_excel('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.xls')
# the [[]] syntax extracts the column(s) into a new dataframe
# we groupby the age column, and 
# apply a count to the age column
ages = df[['age']].groupby('age')['age'].count()
print("The most common ages")
print (ages.sort_values(ascending=False))

在 Jupyter 中的结果显示如下。从数据来看，船上有很多年轻的旅客。考虑到这一点，也就能理解为什么船上有这么多婴儿。

在数据框中计算异常值

我们可以使用标准计算来计算异常值，即判断与均值的差的绝对值是否大于 1.96 倍的标准差。（这假设数据服从正态高斯分布）。

例如，使用之前加载的泰坦尼克号数据集，我们可以根据年龄确定哪些乘客是异常值。

Python 脚本如下：

import pandas as pd

df = pd.read_excel('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.xls')

# compute mean age
df['x-Mean'] = abs(df['age'] - df['age'].mean())

# 1.96 times standard deviation for age
df['1.96*std'] = 1.96*df['age'].std()

# this age is an outlier if abs difference > 1.96 times std dev
df['Outlier'] = abs(df['age'] - df['age'].mean()) > 1.96*df['age'].std()

# print (results)
print ("Dataset dimensions", df.count)
print ("Number of age outliers", df.Outlier.value_counts()[True])

在 Jupyter 下，结果显示如下：

Number of age outliers 65

因为大约有 1300 名乘客，所以我们有大约 5% 的异常值，这意味着年龄可能符合正态分布。

在 Jupyter 中使用 SciPy

SciPy 是一个开源库，适用于数学、科学和工程领域。由于其广泛的应用范围，我们可以使用 SciPy 探索许多领域：

积分
优化
插值法
傅里叶变换
线性代数
还有一些其他强大的功能集，比如信号处理

在 Jupyter 中使用 SciPy 积分

一个标准的数学过程是对方程进行积分。SciPy 通过回调函数来逐步计算出函数的积分。例如，假设我们想要计算以下方程的积分：

我们将使用如下脚本。我们使用标准 math 包中的 pi 定义。

from scipy.integrate import quad
import math

def integrand(x, a, b):
 return a*math.pi + b

a = 2
b = 1
quad(integrand, 0, 1, args=(a,b))

再次强调，这段代码非常简洁，虽然几乎无法在许多语言中实现。在 Jupyter 中运行此脚本时，我们很快就能看到结果：

我很好奇在执行过程中 integrand 函数是如何使用的。我正在利用它来演练回调函数。为了观察其效果，我在脚本中添加了一些调试信息，统计每次调用时发生了多少次迭代，并显示每次调用的值：

from scipy.integrate import quad
import math

counter = 0
def integrand(x, a, b):
 global counter
 counter = counter + 1
 print ('called with x=',x,'a = ',a,'b = ', b)
 return a*math.pi + b

a = 2
b = 1
print(quad(integrand, 0, 1, args=(a,b)))
print(counter)

我们在全局范围内使用了一个计数器，因此在 integrand 函数内引用时需要使用 global 关键字。否则，Python 会假设它是该函数的局部变量。

结果如下：

该函数被调用了 21 次以缩小解决方案范围。

在 Jupyter 中使用 SciPy 优化

通过优化法，我们试图确定在多个变量下函数的最大值或最小值。比如，让我们使用一个包含有趣曲线的方程：

如果我们绘制该曲线并查看是否存在明显的最小值，可以使用以下脚本生成结果图表。（%mathplotlib inline 使得图表在 Jupyter 会话内显示，而不是在新窗口中创建图表。）

%matplotlib inline
from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np

def f(x):
 return x**4 - x**3 + x**2 + 1

x = np.linspace(-100, 50, 100)
plt.plot(x, f(x));

在 Jupyter 中运行这个脚本时，我们看到在 x = 0 处有一个自然最小值。

在 Jupyter 中使用 SciPy 插值法

通过插值法，我们根据一组离散点对函数的值进行猜测。例如，假设你的测试结果显示了如下内容：

%matplotlib inline
import matplotlib.pyplot as plt

x = [1, 3, 5, 7]
y = [0.5, 0.4, 0.35, 0.29]
plt.plot(x,y)

在这种情况下，我们可以使用如下脚本，插值求解当 x 为 4 时的函数值：

from scipy.interpolate import interp1d 
g = interp1d(x, y) 
print (g(4))

这给出了结果 0.375，看起来是正确的。

在 Jupyter 中使用 SciPy 傅里叶变换

SciPy 中有一组 FFT（傅里叶变换）函数。考虑到需要进行的处理量，它们非常易于使用。

我们可以通过以下代码执行快速傅里叶变换（FFT）：

from scipy.fftpack import fft import numpy as np x = np.array([2.0, 1.0, 2.0, 1.0, 2.0]) fft(x)

这里我们有一个小数据集进行分析。这些数据点代表我们需要评估的小信号集。当在 Jupyter 中处理时，我们得到以下显示：

请注意，即使是这个小数据集，变换操作也持续了几秒钟。

我们还可以使用生成的数据集，如下面的代码所示：

from scipy.fftpack import fft import numpy as np
# how many points n = 100 spacing = 1.0 / 250.0 x = np.linspace(0.0, n*spacing, n) y = np.sin(30.0 * np.pi * x) + 0.5 * np.sin(7.0 * np.pi * x) yf = fft(y) xf = np.linspace(0.0, 1.0/(2.0*spacing), n//2)
#plot the data to get a visual import matplotlib.pyplot as plt plt.plot(xf, 2.0/n * np.abs(yf[0:n//2])) plt.grid() plt.show()

在 Jupyter 中运行这个脚本会生成一个数据点的图形，并显示在新的屏幕上：

这与我们预期的结果相符，显示中有一个大波动和一个小波动。

在 Jupyter 中使用 SciPy 线性代数

有一整套线性代数函数可用。例如，我们可以通过以下步骤求解线性系统：

import numpy as np from scipy import linalg
A = np.array([[1, 1], [2, 3]]) print ("A array") print (A)
b = np.array([[1], [2]]) print ("b array") print (b)
solution = np.linalg.solve(A, b) print ("solution ") print (solution)
# validate results print ("validation of solution (should be a 0 matrix)") print (A.dot(solution) – b)

在 Jupyter 中，输出如下所示：

我们通过最终的零矩阵来验证结果。

扩展 Jupyter 中的 pandas 数据框

用于处理数据框的内置函数比我们迄今为止使用的更多。如果我们从本章之前的一个示例中取出 Titanic 数据集（来自 Excel 文件），我们可以使用更多的函数来帮助呈现和处理数据集。

重复一遍，我们使用脚本加载数据集：

import pandas as pd
df = pd.read_excel('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.xls')

然后，我们可以使用 info 函数检查数据框，该函数显示数据框的特征：

df.info()

一些有趣的点如下：

1309 条记录
14 列
body 列中有效数据的字段不多——大多数数据丢失了
能够很好地概述涉及的数据类型

我们还可以使用 describe 函数，它为数据框中的每一列数字型数据提供了统计 breakdown。

df.describe()

这会生成以下表格展示：

对于每一列数字数据，我们有：

计数
均值
标准差
25%、50% 和 75% 百分位点
项目的最小值和最大值

我们可以使用语法 df[12:13] 切片感兴趣的行，其中第一个数字（默认是数据框中的第一行）是要切割的第一行，第二个数字（默认是数据框中的最后一行）是要切割的最后一行。

运行这个切片操作，我们得到了预期的结果：

由于我们在选择数据框的列时实际上是在创建一个新的数据框，因此我们也可以在此基础上使用 head 函数：

在 Jupyter/IPython 中排序和筛选数据框

数据框自动允许你轻松地对数据集进行排序和筛选，使用数据框本身已有的功能。

过滤数据框

我们可以根据条件选择/过滤特定的行（使用相同的 Titanic 数据框）：

print(df[df.age < 5])

这意味着我们查看数据框，并选择年龄小于五岁的人所在的行。（再次强调，这会创建一个新的数据框，可以根据需要进行操作。）

如果你考虑一下这个问题，你几乎可以对数据框应用任何过滤器。然后，你可以做一些事情，比如选择一个数据框的一部分，并将其与另一个数据框的部分进行合并/连接。很快，你就会得出类似 SQL 的操作，可以在数据库表上执行。换个角度看，你将能进行比基础数据框所能提供的更广泛的数据操作。

排序数据框

在大多数语言中，排序意味着重新组织你正在处理的数据集。在数据框中，可以通过选择另一个索引来访问数据框，从而实现排序。所有数据框在一开始都会有一个由 NumPy 内置的基本增量行索引。你可以改变用于访问数据框的索引，并有效地按照你希望的方式对数据框进行排序。

如果我们查看（Titanic）数据框的显示内容，我们会注意到没有命名的第一列序号值：

如果我们为数据框分配另一个索引来使用，我们将按该索引对数据框进行排序。例如：

df.set_index('name').head()

记住，由于我们尚未分配这个新的数据框（命名为 index），我们仍然保留着原始数据框。

与前一节相关，实际上可以对数据框执行排序操作，使用 sort_values 方法。例如，如果我们使用以下脚本：

print(df.sort_values(by='home.dest', ascending=True).head())

这个脚本将数据框按 home.dest 列升序排序，并打印前五条记录（按此顺序）

我们将看到如下结果：

总结

在本章中，我们研究了在 Jupyter 中可能执行的一些计算密集型任务。我们使用 Python 抓取网站以收集数据进行分析。我们使用了 Python 的 NumPy、pandas 和 SciPy 函数来进行深入计算。我们进一步探索了 pandas，并研究了如何操作数据框。最后，我们看到了一些排序和过滤数据框的例子。

在下一章中，我们将做出一些预测，并使用可视化来验证我们的预测。

第三章：数据可视化与预测

做预测通常是危险的。然而，已经有一些方法可以让你对结果充满信心。在 Jupyter 下，我们可以使用 Python 和/或 R 来进行预测，并且这些功能是现成的。

使用 scikit-learn 做预测

scikit-learn 是一个基于 Python 构建的机器学习工具包。该包的一部分是监督学习，其中样本数据点具有属性，允许你将数据点分配到不同的类别中。我们使用一种估算器，它将数据点分配到一个类别，并对其他具有相似属性的数据点做出预测。在 scikit-learn 中，估算器提供了两个函数，fit() 和 predict()，分别提供了分类数据点和预测其他数据点类别的机制。

作为例子，我们将使用来自 uci.edu/ 的住房数据（我认为这些数据来自波士顿地区）。有许多因素，包括价格因素。

我们将采取以下步骤：

我们将把数据集分成训练集和测试集
从训练集中，我们将生成一个模型
我们将使用该模型在测试集上进行评估，并验证我们的模型如何拟合实际数据以预测住房价格

数据集中的属性（在我们的数据框中按相应顺序排列）为：

CRIM	每个城镇的人均犯罪率
`ZN`	住宅用地比例，规划为超过 25,000 平方英尺的地块
`INDUS`	每个城镇非零售商业用地的比例
`CHAS`	查尔斯河虚拟变量（`= 1` 如果地区边界为河流；否则为 `0`）
`NOX`	氮氧化物浓度（每千万分之一）
`RM`	每个住宅的平均房间数
`AGE`	1940 年之前建造的自有住宅单元的比例
`DIS`	到波士顿五个就业中心的加权距离
`RAD`	辐射高速公路的可达性指数
`TAX`	每 $10,000 的全额财产税率
`PTRATIO`	每个城镇的师生比例
`B`	1000(Bk - 0.63)² 其中 Bk 是每个城镇的黑人居民比例
`LSTAT`	低阶层人口的百分比
`MEDV`	自有住宅的中位数价值（单位：$1,000）

以下的代码之后是对所使用的算法和结果的讨论：

#define all the imports we are using import matplotlib.pyplot as plt import numpy as np import pandas as pd import random from sklearn import datasets, linear_model from sklearn.cross_validation import train_test_split # load the data set df = pd.read_table('http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data', sep='\s+') # add column names df.columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', \
 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PRATIO',\ 'B', 'LSTAT', 'MDEV'] #produce basic statistics to make sure things are lined up df.head()

head() 函数的结果是数据框的前几行：

df.describe()

同样，前面的 describe 语句为我们提供了一些数据框的快速统计信息：

当将数据集分割为训练集和测试集时，我们使用两者之间的随机分配。这为我们提供了一组无偏的数据进行处理。然而，为了让你能复现这里展示的结果，你需要使用相同的随机种子/初始值。这就是为什么要调用random.seed()的原因。在实践中，你可以省略这个方法调用：

#we are going to be splitting up the data set 'randomly', #however we need to reproduce results so set the seed random.seed(3277) #split the data into training and testing (25% for testing) training, testing = train_test_split(df, test_size = 0.25) #need this step to create an instance of the lreg model regr = linear_model.LinearRegression() # Train the model using the training set (MDEV=target) training_data = training.drop('MDEV', axis=1) training_test = training.iloc[:,-1] #training.loc[:,['MDEV']] #look at coefficients in the model to validate regr.fit(training_data,training_test) print('Coefficients: \n', regr.coef_) 'Coefficients: \n', array([
 -1.18763385e-01,   4.19752612e-02,  -1.18584543e-02, 5.53125252e-01,  -1.19774970e+01,   3.80050180e+00, -8.55663104e-03,  -1.46613256e+00,   3.86772585e-01, -1.53024705e-02,  -9.55933426e-01,   1.31347272e-02, -5.28183554e-01]))

这些大多数是小数字，除了与第 6 项（房间数）有 3.8 的正相关，以及与第 8 项（从商业中心的距离）有 -1.5 的负相关。人们如此重视靠近工作地点的价值真是有趣：

#split up our test set testing_data = testing.loc[:,['CRIM', 'ZN', 'INDUS', 'CHAS',\ 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PRATIO', 'B',\ 'LSTAT']] testing_test = testing[['MDEV']].as_matrix() #make our prediction prediction_of_test = regr.predict(testing_data) # compute MSE # would usually use the built-in mse function, # but the test_test and prediction have diff # cols sum = 0 rows = len(testing_test) for i in range(rows):
 test = testing_test[i] prediction = prediction_of_test[i] diff = (test - prediction) ** 2 sum = sum + diff mse = sum / rows print("MSE ", mse) ('MSE ', array([ 23.1571225]))

MSE 为 23，相较于处理的数字大小来说，这个值似乎非常低。现在，让我们绘制结果图形，以便直观了解发生了什么：

**#this preceding line is needed to display inline on Jupyter** **#plot the tests and predictions** **plt.scatter(testing_test, prediction_of_test, color='black')
****#draw a line through the middle showing the fit** **x0 = min(testing_test)** **x1 = max(testing_test)** **y0 = min(prediction_of_test)** **y1 = max(prediction_of_test)** **plt.plot([x0,x1],[y0,y1], color="red")
****#add labels** **plt.xlabel("Actual Price")** **plt.ylabel("Predicted Price")** **plt.title("Actual Price vs Predicted Price")
****plt.show()**

**从视觉上看，我们似乎有一个不错的拟合。大部分数据点都与通过的轴对齐。正如往常一样，也有一些明显的异常值，比如 20，50：

** **# 使用 R 进行预测

我们可以使用 R 在笔记本中执行相同的分析。不同语言的函数不同，但功能非常接近。

我们使用相同的算法：

加载数据集
将数据集划分为训练和测试分区
基于训练分区开发一个模型
使用模型从测试分区进行预测
比较预测值与实际测试值

编码如下：

#load in the data set from uci.edu (slightly different from other housing model)
housing <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data")

#assign column names
colnames(housing) <- c("CRIM", "ZN", "INDUS", "CHAS", "NOX",
 "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO",
 "B", "LSTAT", "MDEV")
#make sure we have the right data being loaded
summary(housing)
 CRIM                ZN             INDUS            CHAS 
 Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000 
 1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000 
 Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000 
 Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917 
 3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000 
 Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000 
 NOX               RM             AGE              DIS 
 Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130 
 1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100 
 Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207 
 Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795 
 3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188 
 Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127 
...

确保数据集按正确顺序排列以便进行建模。

housing <- housing[order(housing$MDEV),]

#check if there are any relationships between the data items
plot(housing)

数据显示如下，展示了每个变量与数据集中的其他变量的关系。我在查看是否有任何漂亮的 45 度“线”显示出变量之间的高度对称性，目的是看是否应该移除一个变量，因为另一个变量就足够作为一个贡献因子。有趣的项目是：

CHAS：查尔斯河访问权限，但这是一个二值值。
LSTAT（低收入群体）和 MDEV（价格）有反向关系——但价格不会成为一个因素。
NOX（烟雾）和 DIST（到工作地点的距离）有反向关系。我认为我们需要这个。
否则，数据项之间似乎没有任何关系：

我们像之前一样强制设置随机种子，以便能够复现结果。然后，我们使用 createDataPartitions 函数将数据划分为训练和测试分区。接下来我们可以训练我们的模型并测试结果模型进行验证：

#force the random seed so we can reproduce results
set.seed(133)

#caret package has function to partition data set
library(caret)
trainingIndices <- createDataPartition(housing$MDEV, p=0.75, list=FALSE)
#break out the training vs testing data sets
housingTraining <- housing[trainingIndices,]
housingTesting <- housing[-trainingIndices,]
#note their sizes
nrow(housingTraining)
nrow(housingTesting)
#note there may be warning messages to update packages
381
125

#build a linear model
linearModel <- lm(MDEV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE +
 DIS + RAD + TAX + PRATIO + B + LSTAT, data=housingTraining)
summary(linearModel)
Call:
lm(formula = MDEV ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + 
 DIS + RAD + TAX + PRATIO + B + LSTAT, data = housingTraining)

Residuals:
 Min       1Q   Median       3Q      Max 
-15.8448  -2.7961  -0.5602   2.0667  25.2312 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)  36.636334   5.929753   6.178 1.72e-09 ***
CRIM         -0.134361   0.039634  -3.390 0.000775 ***
ZN            0.041861   0.016379   2.556 0.010997 * 
INDUS         0.029561   0.068790   0.430 0.667640 
CHAS          3.046626   1.008721   3.020 0.002702 ** 
NOX         -17.620245   4.610893  -3.821 0.000156 ***
RM            3.777475   0.484884   7.790 6.92e-14 ***
AGE           0.003492   0.016413   0.213 0.831648 
DIS          -1.390157   0.235793  -5.896 8.47e-09 ***
RAD           0.309546   0.078496   3.943 9.62e-05 ***
TAX          -0.012216   0.004323  -2.826 0.004969 ** 
PRATIO       -0.998417   0.155341  -6.427 4.04e-10 ***
B             0.009745   0.003300   2.953 0.003350 ** 
LSTAT        -0.518531   0.060614  -8.555 3.26e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.867 on 367 degrees of freedom
Multiple R-squared:  0.7327,  Adjusted R-squared:  0.7233 
F-statistic:  77.4 on 13 and 367 DF,  p-value: < 2.2e-16

有趣的是，这个模型也发现查尔斯河景观对价格有高额溢价的影响。此外，这个模型还提供了 p-value（对模型有良好的信心）：

# now that we have a model, make a prediction
predicted <- predict(linearModel,newdata=housingTesting)
summary(predicted)

#visually compare prediction to actual
plot(predicted, housingTesting$MDEV)

它看起来像是一个相当好的相关性，非常接近 45 度的映射。唯一的例外是预测值略高于实际值：

** **# 互动式可视化

有一个 Python 包，Bokeh，可以用来在你的笔记本中生成一个图形，用户可以互动并改变图形。

在这个例子中，我使用了本章稍后直方图示例中的相同数据（也包含在本章的文件集中），以显示一个互动式的 Bokeh 直方图。

编码如下：**

**output_notebook()** **# load the counts from other histogram example** **from_counts = np.load("from_counts.npy")** **# convert array to a dataframe for Histogram** **df = pd.DataFrame({'Votes':from_counts})** **# make sure dataframe is working correctly** **print(df.head())
** **Votes** **0     23** **1     29** **2     23** **3    302** **4     24** **# display the Bokeh histogram** **hist = Histogram(from_counts, \** **title="How Many Votes Made By Users", \** **bins=12)** **show(hist)**

**我们可以看到如下所示的直方图。对于图表的清理，自动执行的操作较少，比如移动计数器或不重要的坐标轴标签。我猜测Histogram函数有选项可以进行进一步的更改：

注意图像顶部的控件小部件：

左侧是一个 Bokeh 图标
右侧是以下图标：
- 将图像移动到屏幕的另一部分
- 放大
- 调整大小
- 滚轮缩放- 滚动滚轮进行放大/缩小
- 将图像保存到磁盘
- 刷新图像
- Bokeh 功能的交互式帮助** **# 使用 Plotly 绘图

Plotly 是一个有趣的混合体。它是一个订阅网站，提供强大的数据分析图形功能。你可以使用该软件的免费版本，但仍需要使用凭证登录。图形功能支持多种语言，从 Node.js 到 Python 等。

此外，生成的图形可以在 Plotly 和本地笔记本中查看。如果你将图形标记为公开，那么就可以像其他互联网上的图形一样从笔记本中访问它。同样，作为网页图形，你可以从显示中选择并根据需要保存到本地。

在这个例子中，我们再次使用投票直方图，但这次使用的是 Plotly 的功能。

脚本如下：

import plotly
import plotly.graph_objs as go
import plotly.plotly as py
import pandas as pd
import numpy as np

#once you set credentials they are stored in local space and referenced automatically
#you need to subscribe to the site to get the credentials
#the api key would need to be replaced with your key
#plotly.tools.set_credentials_file(username='DemoAccount', api_key='lr17zw81')

#we are generating a graphic that anyone can see it
plotly.tools.set_config_file(world_readable=True, sharing='public')

# load voting summary from other project
from_counts = np.load("from_counts.npy")
print(from_counts.shape)
(6110,)

#plotly expects a list in a data block
from_count_list = []
for from_count in from_counts:
 from_count_list.append(from_count)

data = [go.Histogram(x=from_count_list)]

# plot on plot.ly site
py.iplot(data, filename='basic histogram')

我认为这是我见过的使用开箱即用的选项/设置绘制的直方图中最漂亮的一种。我们看到了之前看到的相同直方图，只是采用了更具视觉吸引力的属性：

** **# 创建人类密度图

我原本计划制作一个全球人类密度图，但现有的图形不支持设置每个国家的颜色。因此，我构建了一个美国的密度图。

算法是：

获取每个州的图形形状。
获取每个州的密度。
决定颜色范围，并将最低密度应用于范围的一端，最高密度应用于另一端。
对于每个状态：
- 确定它的密度
- 查找该密度值在范围中的位置并选择颜色
- 绘制状态

这段代码如下（代码进行时嵌入了注释）：**

**import matplotlib.pyplot as plt** **from mpl_toolkits.basemap import Basemap** **from matplotlib.patches import Polygon** **import pandas as pd** **import numpy as np** **import matplotlib** **# create the map** **map = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
** **projection='lcc',lat_1=33,lat_2=45,lon_0=-95)** **# load the shapefile, use the name 'states'** **# download from https://github.com/matplotlib/basemap/tree/master/examples/st99_d00.dbf,shx,shp** **map.readshapefile('st99_d00', name='states', drawbounds=True)** **# collect the state names from the shapefile attributes so we can** **# look up the shape obect for a state by it's name** **state_names = []** **for shape_dict in map.states_info:
** **state_names.append(shape_dict['NAME'])** **ax = plt.gca() # get current axes instance** **# load density data drawn from** **# https://en.wikipedia.org/wiki/List_of_U.S._states_by_population_density** **df = pd.read_csv('states.csv')** **print(df.head())** **State        rank density/mi2  density/km2  pop_rank   2015_pop** **New Jersey      1       1,218          470        11  8,958,013** **Rhode Island    2       1,021          394        43  1,056,298** **Massachusetts   3         871          336        15  6,794,422** **Connecticut     4         741          286        29  3,590,886** **Maryland        5         618          238        19  6,006,401  
****land_rank area_mi2   area_km2** **0         46    7,354  19,046.80** **1         50    1,034   2,678.00** **2         45    7,800  20,201.90** **3         48    4,842  12,540.70** **4         42    9,707  25,141.00** **# determine the range of density values** **max_density = -1.0** **min_density = -1.0** **for index, row in df.iterrows():
** **d = row['density/mi2']** **density = float(d.replace(',' , ''))** **if (max_density==-1.0) or (max_density<density):** **max_density = density** **if (min_density==-1.0) or (min_density>density):** **min_density = density** **print('max',max_density)** **print('min',min_density)** **range_density = max_density - min_density** **print(range_density)** **('max', 1218.0)** **('min', 1.0)** **1217.0** **# we pick a color for the state density out of color map** **cmap = matplotlib.cm.get_cmap('Spectral')** **# for each state get the color for it's density** **for index, row in df.iterrows():
** **state_name = row['State']** **d = row['density/mi2']** **density = float(d.replace(',' , ''))** **color = cmap((density - min_density)/range_density)** **seg = map.states[state_names.index(state_name)]** **poly = Polygon(seg, facecolor=color, edgecolor=color)** **ax.add_patch(poly)** **plt.show()**

**我们在下图中看到了一个颜色编码的密度图。我不确定为什么明尼苏达州和威斯康星州与数据不匹配（它们在地图上显示没有颜色）。数据文件看起来是正确的，并且确实映射到图像点上。

这个示例中使用的包需要安装，因为它们不是标准包的一部分：

** **# 绘制社会数据的直方图

有各种各样的社交网站会生成数据集。在这个示例中，我们将收集其中一个数据集并从数据中生成一个直方图。具体的数据集是来自snap.stanford.edu/data/wiki-Vote.html的 WIKI 投票行为数据。每一项数据展示了用户编号N为用户编号X投票。因此，我们通过生成直方图来分析投票行为，方法如下：

收集所有发生的投票
每一票：
- 增加一个计数器，记录谁投了票
- 增加一个计数器，记录谁被投票
- 对数据进行处理，以便我们可以将其显示在两个直方图中

代码如下：**

**# import all packages being used** **import matplotlib.pyplot as plt** **import pandas as pd** **import numpy as np** **import matplotlib
****# load voting data drawn from https://snap.stanford.edu/data/wiki-Vote.html** **df = pd.read_table('wiki-Vote.txt', sep=r"\s+", index_col=0)
****# produce standard summary info to validate** **print(df.head())** **print(df.describe())**

**Python 会自动将第一列指定为表格的索引，不管索引是否被重用（如本例所示）。你可以在describe()的结果中看到，只有ToNodeId列被提到：

 ToNodeId FromNodeId 30              1412 30              3352 30              5254 30              5543 30              7478
 ToNodeId count  103689.000000 mean     3580.347018 std      2204.045658 min         3.000000 25%      1746.000000 50%      3260.000000 75%      5301.000000 max      8297.000000

接下来，我们根据每个人的投票数量和每个人收到的投票数量来生成分组总计。我认为应该有一个内建函数可以更好地完成这个任务，但我没有找到：

from_counter = {} to_counter = {} for index, row in df.iterrows():
 ton = row['ToNodeId'] fromn = index    #add the from entry
 if from_counter.has_key(fromn): # bump entry from_counter[fromn] = from_counter.get(fromn) + 1 else: # create entry from_counter[fromn] = 1    #add the to entry
 if to_counter.has_key(ton): # bump entry to_counter[ton] = to_counter.get(ton) + 1 else: # create entry to_counter[ton] = 1print(from_counter) print(to_counter) {3: 23, 4: 29, 5: 23, 6: 302, 7: 24, 8: 182, 9: 81, 10: 86, 11: 743,…

我们已经可以看到其中一些较大的数字，比如743：

#extract the count values from_counts = from_counter.values() to_counts = to_counter.values()
print("Most votes by a user",max(from_counts)) print("Most voted for",max(to_counts)) ('Most votes by a user', 893) ('Most voted for', 457)
#make histogram of number of references made by a user plt.hist(from_counts) plt.title("How Many Votes Made By Users") plt.xlabel("Value") plt.ylabel("Frequency") plt.show()

我们看到以下图表，展示了用户投票的熟悉布局。我认为这是我见过的最简洁的布局之一：

现在我们通过以下代码生成一个用户的引用直方图：

#make histogram of number of references made for a user plt.hist(to_counts) plt.title("How Many Votes Made for User") plt.xlabel("Value") plt.ylabel("Frequency") plt.show()

我们看到如下的用户投票图表。我没想到结果会如此失衡：只有少数人投了很多票，只有少数人收到了大量的投票：

** **# 绘制三维数据

许多数据分析包（如 R、Python 等）具有强大的数据可视化功能。一个有趣的功能是将数据展示为三维图形。通常，当使用三维时，会出现意想不到的可视化效果。

在这个示例中，我们使用来自uci.edu/的汽车数据集。它是一个使用广泛的数据集，包含多个车辆属性，例如mpg、weight和acceleration。如果我们将这三项数据属性一起绘制出来，看看能否识别出任何明显的规律呢？

第四章：数据挖掘与 SQL 查询

PySpark 将 Spark 编程模型暴露给 Python。Spark 是一个快速的、通用的大规模数据处理引擎。我们可以在 Jupyter 中使用 Python。因此，我们可以在 Jupyter 中使用 Spark。

安装 Spark 需要在你的机器上安装以下组件：

Java JDK。
从 www.scala-lang.org/download/ 下载 Scala。
推荐从 continuum.io 下载带 Python 的 Anaconda。
从 spark.apache.org/downloads.html 下载 Spark。
winutils：这是一个命令行工具，它将 Linux 命令暴露给 Windows。32 位和 64 位版本可以在以下链接找到：
- 32 位 winutils.exe 可在 code.google.com/p/rrd-hadoop-win32/source/checkout 下载
- 64 位 winutils.exe 可在 github.com/steveloughran/winutils/tree/master/hadoop-2.6.0/bin 下载

然后设置环境变量，指明前述组件的位置：

JAVA_HOME：你安装 JDK 的 bin 目录
PYTHONPATH：Python 安装目录
HADOOP_HOME：winutils 所在的目录
SPARK_HOME：Spark 安装的位置

这些组件可以通过互联网轻松获得，适用于各种操作系统。我已经在 Windows 和 Mac 环境下成功安装了这些组件。

一旦你安装了这些组件，你应该能够在命令行窗口中运行命令 pyspark，并可以使用 Python（并且能够访问 Spark）的 Jupyter Notebook。在我的安装中，我使用了以下命令：

pyspark

由于我将 Spark 安装在根目录的 \spark 目录下。是的，pyspark 是 Spark 的内置工具。

Windows 安装的特别说明

Spark（其实是 Hadoop）需要一个临时存储位置来存放其工作数据集。在 Windows 上，默认位置是 \tmp\hive。如果在 Spark/Hadoop 启动时该目录不存在，它会自动创建。不幸的是，在 Windows 上，安装程序并未内置正确的工具来设置该目录的访问权限。

你应该能够在 winutils 中运行 chmod 来设置 hive 目录的访问权限。然而，我发现 chmod 功能并未正确运行。

一个更好的做法是以管理员模式自行创建 tmp\hive 目录。然后，再次以管理员模式将该目录的完全权限授予所有用户。

如果没有这个更改，Hadoop 会立刻失败。当你启动 pyspark 时，输出（包括任何错误）会显示在命令行窗口中。其中一个错误会是该目录的访问权限不足。

使用 Spark 分析数据

访问 Spark 的第一步是创建一个 SparkContext。SparkContext 初始化所有 Spark 并设置对 Hadoop 的任何访问（如果你也在使用 Hadoop）。

最初使用的对象是 SQLContext，但最近已被弃用，推荐使用 SparkContext，它更为开放。

我们可以使用一个简单的示例来读取文本文件，如下所示：

from pyspark import SparkContext
sc = SparkContext.getOrCreate()

lines = sc.textFile("B05238_04 Spark Total Line Lengths.ipynb")
lineLengths = lines.map(lambda s: len(s))
totalLength = lineLengths.reduce(lambda a, b: a + b)
print(totalLength)

在此示例中：

我们获取一个 SparkContext
使用上下文，读取一个文件（本示例中的 Jupyter 文件）
我们使用 Hadoop map 函数将文本文件拆分为不同的行并收集行的长度。
我们使用 Hadoop reduce 函数计算所有行的长度。
我们展示我们的结果

在 Jupyter 中，效果如下所示：

另一个 MapReduce 示例

我们可以在另一个示例中使用 MapReduce，从文件中获取单词计数。这是一个标准问题，但我们使用 MapReduce 来完成大部分繁重的工作。我们可以使用此示例的源代码。我们可以使用类似这样的脚本来计算文件中单词的出现次数：

import pyspark
if not 'sc' in globals():
 sc = pyspark.SparkContext()

text_file = sc.textFile("Spark File Words.ipynb")
counts = text_file.flatMap(lambda line: line.split(" ")) \
 .map(lambda word: (word, 1)) \
 .reduceByKey(lambda a, b: a + b)
for x in counts.collect():
 print x

我们有相同的代码前言。

然后我们将文本文件加载到内存中。

text_file 是一个 Spark RDD（弹性分布式数据集），而不是数据框。

假设文件是巨大的，且内容分布在多个处理程序中。

一旦文件加载完成，我们将每行拆分成单词，然后使用 lambda 函数标记每个单词的出现。代码实际上为每个单词出现创建了一个新记录，例如 at appears 1，at appears 1。例如，如果单词 at 出现了两次，每次出现都会添加一个记录，像 at appears 1。这个想法是暂时不聚合结果，只记录我们看到的出现情况。这个过程可以被多个处理器分担，每个处理器生成这些低级信息。我们并不关心优化这个过程。

一旦我们拥有所有这些记录，我们就根据提到的单词出现次数来减少/汇总记录集。

counts 对象在 Spark 中也是 RDD。最后的 for 循环会对 RDD 执行 collect()。如前所述，这个 RDD 可能分布在多个节点上。collect() 函数会将 RDD 的所有副本拉到一个位置。然后我们遍历每条记录。

当我们在 Jupyter 中运行时，我们会看到类似这样的显示：

该列表被简化，因为单词列表会继续一段时间。

前面的示例在 Python 3 中效果不佳。直接在 Python 中编写代码时有解决方法，但在 Jupyter 中没有。

使用 SparkSession 和 SQL

Spark 提供了许多类似 SQL 的操作，可以对数据框执行。例如，我们可以加载一个包含产品销售信息的 CSV 文件数据框：

from pyspark.sql import SparkSession 
spark = SparkSession(sc) 

df = spark.read.format("csv") \
 .option("header", "true") \
 .load("productsales.csv");
df.show()

示例：

启动一个 SparkSession（大多数数据访问需要）
使用会话读取一个包含头部记录的 CSV 格式文件
显示初始行

在销售数据中，我们有一些有趣的列：

按部门划分的产品实际销售额
按部门划分的产品预测销售额

如果这是一个更大的文件，我们可以使用 SQL 来确定产品列表的范围。接下来是用 Spark SQL 确定产品列表的代码：

df.groupBy("PRODUCT").count().show()

数据框的groupBy函数与 SQL 中的Group By子句非常相似。Group By根据指定列中的值将数据集中的项目进行分组。在本例中是PRODUCT列。Group By的结果是生成一个包含结果的数据集。作为数据集，我们可以使用count()函数查询每个组中的行数。

所以，groupBy的结果是对应分组元素的项目数量。例如，我们按CHAIR分组，发现有 288 个这样的项目：

所以，显然我们没有真实的产品数据。任何公司不太可能在每个产品线中有完全相同数量的产品。

我们可以查看数据集，以使用filter()命令查看不同部门在实际销售与预测销售方面的表现，示例如下：

df.filter(df['ACTUAL'] > df['PREDICT']).show()

我们将一个逻辑测试传递给filter命令，该命令将在数据集的每一行上执行。如果该行数据通过测试，则返回该行；否则，该行会被从结果中删除。

我们的测试只关心实际销售额超过预测值的情况。

在 Jupyter 下，显示如下：

因此，我们得到了一个缩小的结果集。再次强调，这个结果是由filter函数生成的一个数据框，并可以像其他数据框一样调用show。注意，之前显示的第三条记录不再出现，因为它的实际销售额低于预测值。通常，做一个快速调查是个好主意，以确保得到正确的结果。

如果我们想进一步追踪，找出公司内表现最佳的区域该怎么办？

如果这是一个数据库表格，我们可以创建另一列，存储实际销售与预测销售之间的差异，然后根据这一列对结果进行排序。在 Spark 中我们可以执行类似的操作。

使用数据框，我们可以用类似这样的代码：

#register dataframe as temporary SQL table
df.createOrReplaceTempView("sales")
# pull the values by the difference calculated
sqlDF = spark.sql("SELECT *, ACTUAL-PREDICT as DIFF FROM sales ORDER BY DIFF desc")
sqlDF.show()

第一条语句是在当前上下文中创建一个视图/数据框，以便进一步操作。此视图是惰性计算的，除非采取特定步骤，否则不会持久化，最重要的是，可以作为 Hive 视图访问。该视图可以直接通过SparkContext访问。

然后，我们使用创建的新销售视图创建一个新的数据框，并计算新的列。在 Jupyter 下，显示如下：

再次，我认为我们没有现实的值，因为差异与预测值相差甚远。

创建的数据框是不可变的，不像数据库表。

合并数据集

所以，我们已经看到将数据框移动到 Spark 进行分析。这看起来非常接近 SQL 表。在 SQL 中，标准做法是不在不同的表中重复项。例如，产品表可能会有价格，订单表只会通过产品标识符引用产品表，以避免重复数据。因此，另一个 SQL 做法是将表连接或合并，以便获得完整的所需信息。继续使用订单的类比，我们将所有涉及的表合并，因为每个表都有订单完成所需的数据。

创建一组表并使用 Spark 连接它们会有多难？我们将使用 Product、Order 和 ProductOrder 示例表：

表格	列
产品	产品 ID, 描述, 价格
订单	订单 ID, 订单日期
产品订单	订单 ID, 产品 ID, 数量

所以，一个 Order 关联有一组 Product/Quantity 值。

我们可以填充数据框并将它们移入 Spark：

from pyspark import SparkContext
from pyspark.sql import SparkSession

sc = SparkContext.getOrCreate()
spark = SparkSession(sc) 

# load product set
productDF = spark.read.format("csv") \
 .option("header", "true") \
 .load("product.csv");
productDF.show()
productDF.createOrReplaceTempView("product")

# load order set
orderDF = spark.read.format("csv") \
 .option("header", "true") \
 .load("order.csv");
orderDF.show()
orderDF.createOrReplaceTempView("order")

# load order/product set
orderproductDF = spark.read.format("csv") \
 .option("header", "true") \
 .load("orderproduct.csv");
orderproductDF.show()
orderproductDF.createOrReplaceTempView("orderproduct")

现在，我们可以尝试在它们之间执行类似 SQL 的 JOIN 操作：

# join the tables
joinedDF = spark.sql("SELECT * " \
 "FROM orderproduct " \
 "JOIN order ON order.orderid = orderproduct.orderid " \
 "ORDER BY order.orderid")
joinedDF.show()

在 Jupyter 中执行这些操作的结果如下所示：

我们的标准导入获取了一个 SparkContext 并初始化了一个 SparkSession。注意 SparkContext 的 getOrCreate 方法。如果你在 Jupyter 外运行这段代码，则没有上下文，会创建一个新的上下文。在 Jupyter 中，Spark 启动时会为所有脚本初始化一个上下文。我们可以在任何 Spark 脚本中使用该上下文，而不必自己创建一个。

加载我们的 product 表：

加载 order 表：

加载 orderproduct 表。请注意，至少有一个订单包含多个产品：

我们在结果集中有来自 order 和 orderproduct 的 orderid 列。我们可以在查询中进行更精确的筛选，指定我们想要返回的确切列：

我曾尝试使用 Spark 的 join() 命令，但未能成功。

我在网上找到的文档和示例都是过时、稀疏且不正确的。

使用该命令也显示了一个持续的错误，即任务没有及时返回结果。从底层的 Hadoop 看，我预计处理任务通常会被拆分为独立的任务。我猜 Spark 也会类似地将函数拆分到不同的线程中去完成。为什么这些小任务没有完成并不清楚，因为我并没有要求它执行任何超常的操作。

将 JSON 加载到 Spark 中

Spark 还可以访问 JSON 数据进行操作。这里我们有一个示例：

将 JSON 文件加载到 Spark 数据框中
检查数据框的内容并显示其明显的模式
与前面提到的其他数据框一样，将数据框移入上下文以便 Spark 会话可以直接访问
显示一个在 Spark 上下文中访问数据框的示例

列表如下：

我们的标准包括 Spark 的内容：

from pyspark import SparkContext
from pyspark.sql import SparkSession 
sc = SparkContext.getOrCreate()
spark = SparkSession(sc)

读取 JSON 并显示我们找到的内容：

#using some data from file from https://gist.github.com/marktyers/678711152b8dd33f6346
df = spark.read.json("people.json")
df.show()

我在将标准 JSON 加载到 Spark 中时遇到了一些困难。Spark 似乎期望 JSON 文件中的每个列表记录都是一条数据记录，而我所见的许多 JSON 文件格式将记录布局进行缩进等处理。

注意使用了 null 值，表示某个实例未指定属性。

显示数据的解释模式：

df.printSchema()

所有列的默认值为nullable。你可以更改列的属性，但不能更改列的值，因为数据值是不可变的。

将数据框移入上下文并从中访问：

df.registerTempTable("people")
spark.sql("select name from people").show()

此时，people表像 Spark 中的任何其他临时 SQL 表一样工作。

使用 Spark 透视表

pivot()函数允许你将行转换为列，同时对某些列执行聚合操作。如果你仔细想想，这就像是围绕一个支点调整表格的坐标轴。

我想到了一个简单的例子来展示这一切是如何工作的。我认为这是一项功能，一旦你看到它的实际应用，就会意识到可以在很多地方应用它。

在我们的示例中，我们有一些原始的股票价格数据，我们希望通过透视表将其转换，按股票每年生成平均价格。

我们示例中的代码是：

from pyspark import SparkContext
from pyspark.sql import SparkSession
from pyspark.sql import functions as func

sc = SparkContext.getOrCreate()
spark = SparkSession(sc)

# load product set
pivotDF = spark.read.format("csv") \
 .option("header", "true") \
 .load("pivot.csv");
pivotDF.show()
pivotDF.createOrReplaceTempView("pivot")

# pivot data per the year to get average prices per stock per year
pivotDF \
 .groupBy("stock") \
 .pivot("year",[2012,2013]) \
 .agg(func.avg("price")) \
 .show()

在 Jupyter 中的展示如下：

所有标准包括我们需要的内容，以便 Spark 初始化SparkContext和SparkSession：

我们从 CSV 文件加载股票价格信息。重要的是，至少有一只股票在同一年份有多个价格：

我们按股票符号对信息进行分组。透视表基于数据集中包含两个值的年份——2012 年和 2013 年。我们正在计算每年每只股票的平均价格。

总结

在本章中，我们熟悉了如何获取SparkContext。我们看到了使用 Hadoop MapReduce 的示例。我们用 Spark 数据进行了 SQL 操作。我们结合了数据框，并对结果集进行了操作。我们导入了 JSON 数据并用 Spark 进行了处理。最后，我们查看了如何使用透视表来汇总数据框的信息。

在下一章，我们将学习如何在 Jupyter 中使用 R 编程。

第五章：R 与 Jupyter

在本章中，我们将使用 R 语言代码在 Jupyter 中进行编程。我认为 R 是 Jupyter 中预期使用的主要语言之一。Jupyter 用户可以使用 R 的完整功能。

如何为 Jupyter 设置 R

过去，必须安装 Jupyter、Python 等的独立组件才能构建一个工作系统。但使用 Continuum Analytics 后，安装 Jupyter 并将 R 引擎添加到 Jupyter 的解决方案集中的过程变得非常简单，并且支持 Windows 和 Mac 系统。

假设你已经安装了 conda，我们只需要一个命令就能为 Jupyter 添加对 R 编程的支持：

conda install -c r r-essentials

此时，当你启动 Jupyter 时，列出的内核中将会有 R。

2016 年美国选举人口统计数据分析

为了了解 R 开发者可用的资源，我们可以查看 2016 年的选举数据。在这个案例中，我引用了维基百科的资料(en.wikipedia.org/wiki/United_States_presidential_election,_2016)，特别是名为“2016 年各人口子群体的总统选举投票”表格。以下是我们编写的代码：

定义一个辅助函数，以便我们可以轻松地打印出值。新的printf函数接受传入的任何参数(...)并将其传递给sprintf：

printf <- function(...)print(sprintf(...))

我已将不同的人口统计统计数据存储在不同的TSV（制表符分隔值）文件中，可以使用以下代码读取。对于每个表格，我们使用read.csv函数，并将字段分隔符指定为制表符，而不是默认的逗号。然后，我们使用head函数显示有关已加载数据框的信息：

age <- read.csv("Documents/B05238_05_age.tsv", sep="\t")
head(age)

education <- read.csv("Documents/B05238_05_education.tsv", sep="\t")
head(education)

gender <- read.csv("Documents/B05238_05_gender.tsv", sep="\t")
head(gender)

ideology <- read.csv("Documents/B05238_05_ideology.tsv", sep="\t")
head(ideology)

income <- read.csv("Documents/B05238_05_income.tsv", sep="\t")
head(income)

orientation <- read.csv("Documents/B05238_05_orientation.tsv", sep="\t")
head(orientation)

party <- read.csv("Documents/B05238_05_party.tsv", sep="\t")
head(party)

race <- read.csv("Documents/B05238_05_race.tsv", sep="\t")
head(race)

region <- read.csv("Documents/B05238_05_region.tsv", sep="\t")
head(region)

religion <- read.csv("Documents/B05238_05_religion.tsv", sep="\t")
head(religion)

顶部显示了各列的名称。每一行显示该行的各列值。每次read操作的结果如下所示：

现在，我们可以找出每个投票率的主导特征：

printf("Most Clinton voters from %s",age[which.max(age$Clinton),'age'])
printf("Most Clinton voters from %s",education[which.max(education$Clinton),'education'])
printf("Most Clinton voters from %s",gender[which.max(gender$Clinton),'gender'])
printf("Most Clinton voters from %s",ideology[which.max(ideology$Clinton),'ideology'])
printf("Most Clinton voters from %s",income[which.max(income$Clinton),'income'])
printf("Most Clinton voters from %s",orientation[which.max(orientation$Clinton),'orientation'])
printf("Most Clinton voters from %s",party[which.max(party$Clinton),'party'])
printf("Most Clinton voters from %s",race[which.max(race$Clinton),'race'])
printf("Most Clinton voters from %s",region[which.max(region$Clinton),'region'])
printf("Most Clinton voters from %s",religion[which.max(religion$Clinton),'religion'])

printf("Most Trump voters from %s",age[which.max(age$Trump),'age'])
printf("Most Trump voters from %s",education[which.max(education$Trump),'education'])
printf("Most Trump voters from %s",gender[which.max(gender$Trump),'gender'])
printf("Most Trump voters from %s",ideology[which.max(ideology$Trump),'ideology'])
printf("Most Trump voters from %s",income[which.max(income$Trump),'income'])
printf("Most Trump voters from %s",orientation[which.max(orientation$Trump),'orientation'])
printf("Most Trump voters from %s",party[which.max(party$Trump),'party'])
printf("Most Trump voters from %s",race[which.max(race$Trump),'race'])
printf("Most Trump voters from %s",region[which.max(region$Trump),'region'])
printf("Most Trump voters from %s",religion[which.max(religion$Trump),'religion'])

克林顿的结果如下：

特朗普的结果如下：

有趣的是，支持两位候选人的主要群体之间没有重叠。我认为各党派故意针对不同的群体，以避免重叠。一定有一套花钱的指南，通过广告向已知的、感兴趣的群体宣传，而不是争夺那些模糊不清的人口群体。

分析 2016 年选民登记与投票情况

同样，我们可以查看选民登记与实际投票（使用来自www.census.gov/data/tables/time-series/demo/voting-and-registration/p20-580.html的普查数据）。

首先，我们加载数据集并显示头部信息，以便直观地检查加载是否准确：

df <- read.csv("Documents/B05238_05_registration.csv")
summary(df)

因此，我们有一些按州分类的注册和投票信息。使用 R 自动使用plot命令以x和y格式绘制所有数据：

plot(df)

我们特别关注注册选民与实际投票之间的关系。我们可以从以下图表中看到，大多数数据高度相关（如大多数关系的 45 度角所示）：

我们可以使用 Python 产生类似的结果，但图形显示却相差甚远。

导入我们用于示例的所有包：

from numpy import corrcoef, sum, log, arange from numpy.random import rand from pylab import pcolor, show, colorbar, xticks, yticks import pandas as pd import matplotlib from matplotlib import pyplot as plt

在 Python 中读取 CSV 文件非常类似。我们调用 pandas 来读取文件：

df

如果数据框中存在字符串数据，pandas 将抛出错误，因此只需删除该列（带有州名）：

del df['state'] #pandas do not deal well with strings

一个近似的 Python 函数是corr()函数，它打印出数据框中所有交叉相关项的数值。你需要扫描数据，寻找接近1.0的相关值：

#print cross-correlation matrix print(df.corr())

类似地，我们有corrcoef()函数，它为数据框中类似相关的项提供颜色强度。我没有找到标记相关项的方法：

#graph same fig = plt.gcf() fig.set_size_inches(10, 10)
# plotting the correlation matrix R = corrcoef(df) pcolor(R) colorbar() yticks(arange(0.5,10.5),range(0,10)) xticks(arange(0.5,10.5),range(0,10)) show()

我们希望看到注册和投票之间的实际数值相关性。我们可以通过调用cor函数传入两个感兴趣的数据点来做到这一点，如下所示：

cor(df$voted,df$registered)
0.998946393424037

实际的相关值可能因机器类型而异。这会对后续数值产生影响。

由于相关性达到 99%，我们几乎完美。

我们可以使用数据点来得出回归线，使用lm函数，其中我们说明lm(y ~ (predicted by) x)：

fit <- lm(df$voted ~ df$registered)
fit
Call:
lm(formula = df$voted ~ df$registered)

Coefficients:
 (Intercept)  df$registered 
 -4.1690         0.8741

根据注册数，我们将其乘以 87%并减去 4 以获取实际选民人数。再次强调，数据是相关的。

我们可以通过调用plot函数并传入fit对象来显示回归线的特征（par函数用于布局输出，本例中为 2x2 矩阵样式的四个图形）：

par(mfrow=c(2,2))
plot(fit)

然后，在 2x2 显示的第二部分中，我们有这两个图形：

从前述图表中我们可以看到以下内容：

残差值在非常大的数值之前接近于零。
大部分范围的理论量化与标准化残差非常匹配。
拟合值与标准化残差相比我预期的不是那么接近。
标准化残差与杠杆的关系显示出大多数数值之间的紧密配置。

总体而言，我们对数据仍然有一个非常好的模型。我们可以通过调用residuals函数查看回归的残差如何：

... (for all 50 states)

残差值比我预期的要大，对于如此高度相关的数据，我本以为会看到非常小的数字，而不是数百的值。

一如既往，我们展示我们得到的模型摘要：

summary(fit)
Call:
lm(formula = df$voted ~ df$registered)

Residuals:
 Min      1Q  Median      3Q     Max 
-617.33  -29.69    0.83   30.70  351.27 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)   -4.169018  25.201730  -0.165    0.869 
df$registered  0.874062   0.005736 152.370   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 127.9 on 49 degrees of freedom
Multiple R-squared:  0.9979,  Adjusted R-squared:  0.9979 
F-statistic: 2.322e+04 on 1 and 49 DF,  p-value: < 2.2e-16

你可能会看到不同的结果，具体取决于你运行脚本的机器类型。

在此展示的模型中有许多数据点：

再次，残差值显示出相当大的范围。
系数与我们之前看到的相同，但截距的标准误差较高。
预期 R 平方接近 1。
p值最小是预期的。

我们也可以使用 Python 来得出线性回归模型：

import numpy as np
import statsmodels.formula.api as sm

model = sm.ols(formula='voted ~ registered', data=df)
fitted = model.fit()
print (fitted.summary())

我们看到回归结果呈标准格式如下：

警告信息提示数据中存在一些问题：

我们直接使用协方差矩阵，因此不清楚如果要以其他方式指定该矩阵该怎么做。
我想数据中可能存在较强的多重共线性，因为数据只有这两个项。

我们还可以使用 Python 脚本绘制实际值与拟合值的对比图：

plt.plot(df['voted'], df['registered'], 'ro')
plt.plot(df['voted'], fitted.fittedvalues, 'b')
plt.legend(['Data', 'Fitted model'])
plt.xlabel('Voted')
plt.ylabel('Registered')
plt.title('Voted vs Registered')
plt.show()

我觉得我更喜欢 Python 版本的图表，而不是 R 的显示。左下角添加的强度圆圈让我感到困惑。

分析大学招生变化

我们可以查看近几年大学录取率的趋势。为此分析，我使用了www.ivywise.com/ivywise-knowledgebase/admission-statistics上的数据。

首先，我们读取数据集并显示摘要信息，从头部数据验证：

df <- read.csv("Documents/acceptance-rates.csv")
summary(df)
head(df)

我们看到学校录取率的摘要数据如下：

有趣的是，录取率差异如此之大，从 2017 年的 5%低点到 41%的高点。

让我们再次查看数据图，以验证数据点是否正确：

plot(df)

从相关性图表来看，似乎无法使用 2007 年的数据点。图表显示 2007 年与其他年份之间存在较大的差异，而其他三年的相关性较好。

所以，我们有来自 25 所美国主要大学的连续 3 年数据。我们可以通过几个步骤将数据转换为时间序列。

首先，我们创建一个包含 2015 至 2017 年期间这些高校平均接受率的向量。我们使用均值函数来确定数据框中所有高校的平均值。数据框中有一些NA值，因此我们需要告诉均值函数忽略这些值（na.rm=TRUE）：

myvector <- c(mean(df[["X2015"]],na.rm=TRUE),
 mean(df[["X2016"]],na.rm=TRUE),
 mean(df[["X2017"]],na.rm=TRUE))

接下来，我们将该向量转换为时间序列。时间序列通过传入向量来使用起始和结束点，以及数据点的频率。在我们的例子中，频率是每年一次，所以frequency = 1：

ts <- ts(myvector, start=c(2015), end=c(2017), frequency=1)

然后绘制时间序列，以获得一个好的可视化效果：

plot(ts)

因此，明显的趋势是全面降低接受率，正如我们看到的，初始接受率为.15，在 2017 年稳步下降至.14。

数据看起来非常好，非常合适，因为数据点排列整齐。我们可以使用这个时间序列来预测未来几年。Holt-Winters 算法的版本可以基于级别数据、级别数据加上趋势分量以及级别数据加上趋势分量和季节性分量进行预测。我们有一个趋势，但没有季节性：

# double exponential - models level and trend
fit <- HoltWinters(ts, gamma=FALSE)
fit
Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
HoltWinters(x = ts, gamma = FALSE)

Smoothing parameters:
 alpha: 0.3
 beta : 0.1
 gamma: FALSE

Coefficients:
 [,1]
a  0.14495402
b -0.00415977

我们的指数平滑系数为七分之一，接近零意味着我们并没有大幅降低录取率，但它们正在下降。

现在我们有了现有数据的良好时间序列模型，我们可以预测未来三年的数据并绘制出来：

install.packages("forecast", repos="http://cran.us.r-project.org")
library(forecast)
forecast(fit, 3)
plot(forecast(fit, 3))

趋势显然是负面的，但如前所述，这不是一个急剧的下降，大约每年下降半个百分点。我们还可以查看 Python 的类似编码，如下所示。导入我们将使用的所有 Python 包：

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

将大学录取数据读入数据框架：

data = pd.read_csv('Documents/acceptance-rates.csv')
print (data.head())
 School  2017  2016  2015  2007
0      Amherst College   NaN  0.14  0.14  0.18
1       Boston College  0.32  0.32  0.28  0.27
2     Brown University  0.08  0.09  0.09  0.15
3  Columbia University  0.06  0.06  0.06  0.12
4   Cornell University  0.13  0.14  0.15  0.21

删除School列，因为 Python 无法从字符串中计算：

del data['School']
print (data.head())
 2017  2016  2015  2007
0   NaN  0.14  0.14  0.18
1  0.32  0.32  0.28  0.27
2  0.08  0.09  0.09  0.15
3  0.06  0.06  0.06  0.12
4  0.13  0.14  0.15  0.21

按年份将数据转换为集合：

data = data.transpose()
print (data.head())

我们看到数据集转置为我们期望的形状如下：

查看数据的样子：

plt.plot(data);

我们看到相同的录取率略微下降趋势。

在 Python 中使用 Holt-Winters 预测存在问题，因为它需要进一步转换数据。总体而言，在前一节中 R 中的相同处理变得更加复杂。

预测飞机到达时间

R 具有内置功能，可以在训练集和测试集之间拆分数据框架，基于训练集构建模型，使用模型和测试集预测结果，然后可视化模型的工作效果。

对于这个例子，我使用了 2008 年从stat-computing.org/dataexpo/2009/the-data.html获取的航空公司到达和出发时间与预定到达和出发时间。数据集分发为.bz2文件，解压后为 CSV 文件。我喜欢这个数据集，因为初始行数超过 700 万，并且在 Jupyter 中运行得非常顺利。

我们首先读取飞机数据并显示摘要。数据集中还有其他我们没有使用的列：

df <- read.csv("Documents/2008-airplane.csv")
summary(df)
...
CRSElapsedTime      AirTime          ArrDelay          DepDelay 
 Min.   :-141.0   Min.   :   0     Min.   :-519.00   Min.   :-534.00 
 1st Qu.:  80.0   1st Qu.:  55     1st Qu.: -10.00   1st Qu.:  -4.00 
 Median : 110.0   Median :  86     Median :  -2.00   Median :  -1.00 
 Mean   : 128.9   Mean   : 104     Mean   :   8.17   Mean   :   9.97 
 3rd Qu.: 159.0   3rd Qu.: 132     3rd Qu.:  12.00   3rd Qu.:   8.00 
 Max.   :1435.0   Max.   :1350     Max.   :2461.00   Max.   :2467.00 
 NA's   :844      NA's   :154699   NA's   :154699    NA's   :136246 
 Origin             Dest            Distance          TaxiIn 
 ATL    : 414513   ATL    : 414521   Min.   :  11.0   Min.   :  0.00 
 ORD    : 350380   ORD    : 350452   1st Qu.: 325.0   1st Qu.:  4.00 
 DFW    : 281281   DFW    : 281401   Median : 581.0   Median :  6.00 
 DEN    : 241443   DEN    : 241470   Mean   : 726.4   Mean   :  6.86 
 LAX    : 215608   LAX    : 215685   3rd Qu.: 954.0   3rd Qu.:  8.00 
 PHX    : 199408   PHX    : 199416   Max.   :4962.0   Max.   :308.00 
 (Other):5307095   (Other):5306783                    NA's   :151649 
 TaxiOut         Cancelled       CancellationCode    Diverted 
 Min.   :  0.00   Min.   :0.00000    :6872294        Min.   :0.000000 
 1st Qu.: 10.00   1st Qu.:0.00000   A:  54330        1st Qu.:0.000000 
 Median : 14.00   Median :0.00000   B:  54904        Median :0.000000 
 Mean   : 16.45   Mean   :0.01961   C:  28188        Mean   :0.002463 
 3rd Qu.: 19.00   3rd Qu.:0.00000   D:     12        3rd Qu.:0.000000 
 Max.   :429.00   Max.   :1.00000                    Max.   :1.000000 
 NA's   :137058 
 CarrierDelay      WeatherDelay        NASDelay       SecurityDelay 
 Min.   :   0      Min.   :   0      Min.   :   0      Min.   :  0 
 1st Qu.:   0      1st Qu.:   0      1st Qu.:   0      1st Qu.:  0 
 Median :   0      Median :   0      Median :   6      Median :  0 
 Mean   :  16      Mean   :   3      Mean   :  17      Mean   :  0 
 3rd Qu.:  16      3rd Qu.:   0      3rd Qu.:  21      3rd Qu.:  0 
 Max.   :2436      Max.   :1352      Max.   :1357      Max.   :392 
 NA's   :5484993   NA's   :5484993   NA's   :5484993   NA's   :5484993 
 LateAircraftDelay
 Min.   :   0 
 1st Qu.:   0 
 Median :   0 
 Mean   :  21 
 3rd Qu.:  26 
 Max.   :1316 
 NA's   :5484993

许多数据点具有NA值。我们需要删除这些值以建立准确的模型：

# eliminate rows with NA values
df <- na.omit(df)

让我们创建我们的分区：

# for partitioning to work data has to be ordered
times <- df[order(df$ArrTime),]
nrow(times)
1524735

# partition data - 75% training
library(caret)
set.seed(1337)
trainingIndices <- createDataPartition(df$ArrTime,p=0.75,list=FALSE)
trainingSet <- df[trainingIndices,]
testingSet <- df[-trainingIndices,]
nrow(trainingSet)
nrow(testingSet)
1143553
381182

让我们根据以下字段建立到达时间(ArrTime)的模型：

CRSArrTime: 计划到达时间
ArrDelay: 到达延误
DepDelay: 出发延误
Diverted: 飞机是否使用了转移路线
CarrierDelay: 航空公司系统的延误
WeatherDelay：由于天气原因的延误
NASDelay：由于 NAS 的延误
SecurityDelay：由于安全原因的延误
LateAircraftDelay：由于其他延误导致飞机晚点

其中有两个数据项只是标志（0/1），很不幸。最大的预测因素似乎是计划的到达时间。其他各种延误因素的影响较小。我认为这就像是安检或类似的事情多花了额外的 20 分钟；对于旅行者来说，这是件大事。

现在我们有了模型，让我们使用测试集来进行预测：

predicted <- predict(model, newdata=testingSet)
summary(predicted)
summary(testingSet$ArrTime)
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -941    1360    1629    1590    1843    2217 
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 1    1249    1711    1590    2034    2400

绘制预测数据与实际数据的对比图，以评估模型的准确性：

plot(predicted,testingSet$ArrTime)

从视觉上看，预测数据与实际数据非常吻合，几乎呈 45 度线对齐。然而，图表右下角的整组预测点是令人困扰的。似乎有很多预测数据远低于实际值。我认为必定还有其他因素在影响预测结果，因为我本来预期所有数据会聚集在一个区域，而不是分布在两个区域。

总结

在本章中，我们首先将 R 设置为可用于笔记本的引擎之一。然后我们使用一些基础的 R 语言分析了总统选举中的选民人口统计数据。我们分析了选民登记与实际投票情况。接下来，我们分析了大学招生趋势。最后，我们研究了使用预测模型来判断航班是否会延误。

在下一章，我们将探讨在 Jupyter 中以不同方式处理数据。

第六章：数据清理

在本章中，我们查看了不同形式的数据，并提取了有用的统计信息。用于访问数据的工具已经得到很好的开发，并允许某些记录中缺失标题或数据点。

读取 CSV 文件

文件格式的标准之一是 CSV。在本节中，我们将通过读取 CSV 文件并调整数据集的过程，得出一些关于数据的结论。我使用的数据来自加利福尼亚房屋的供暖系统选择数据集，网址为 vincentarelbundock.github.io/Rdatasets/datasets.html：

#read in the CSV file as available on the site
heating <- read.csv(file="Documents/heating.csv", header=TRUE, sep=",")
# make sure the data is laid out the way we expect
head(heating)

数据似乎符合预期；然而，许多列有缩写名称，并且有些列内容重复。让我们修改我们感兴趣的名称，使其更易读，并移除我们不打算使用的多余内容：

# change the column names to be more readable
colnames(heating)[colnames(heating)=="depvar"] <- "system"
colnames(heating)[colnames(heating)=="ic.gc"] <- "install_cost"
colnames(heating)[colnames(heating)=="oc.gc"] <- "annual_cost"
colnames(heating)[colnames(heating)=="pb.gc"] <- "ratio_annual_install"

# remove columns which are not used
heating$idcase <- NULL
heating$ic.gr <- NULL
heating$ic.ec <- NULL
heating$ic.hp <- NULL
heating$ic.er <- NULL
heating$oc.gr <- NULL
heating$oc.ec <- NULL
heating$oc.hp <- NULL
heating$oc.er <- NULL
heating$pb.gr <- NULL
heating$pb.ec <- NULL
heating$pb.er <- NULL
heating$pb.hp <- NULL

# check the data layout again now that we have made changes
head(heating)

现在我们有了一个更精简的数据集，让我们开始查看数据：

# get rough statistics on the data
summary(heating)

总结中有几点突出：

有五种不同类型的供暖系统，其中以燃气制冷最为普遍
成本变化远比预期的要大
数据涵盖了加利福尼亚的四个大区域
年度成本与初始成本的比率变化远比预期的要大

数据关系并不显而易见，但我们可以使用 R 中的 plot() 函数提供一个快速的快照，显示任何重要的内容：

plot(heating)

再次，几个有趣的事实浮现出来：

初始成本在不同类型的系统中变化很大
年度成本在不同类型的系统中也有所不同
成本在客户收入、年龄、房屋房间数和地区范围内变化很大

变量之间唯一的直接关系似乎是系统的初始成本和年度成本。通过协方差，我们可以衡量两个变量之间如何相互变化。如果我们对安装成本和年度成本进行协方差分析，我们得到：

cov(heating$install_cost, heating$annual_cost) 
2131

我不确定我是否见过更高的协方差结果。

读取另一个 CSV 文件

我们可以查看同一数据集中的另一个 CSV 文件，看看我们会遇到什么问题。使用我们之前从同一网站下载的所有大联盟棒球运动员的年度击球记录，我们可以使用类似以下的代码开始分析数据：

players <- read.csv(file="Documents/baseball.csv", header=TRUE, sep=",")head(players)

这会生成以下的表头显示：

这个数据集中包含了很多棒球运动员的统计数据。也有很多 NA 值。R 对 NA 值的处理非常好。让我们首先使用以下方法查看数据的统计信息：

summary(players)

这会生成所有相关字段的统计数据（还有一些未在此显示的字段）：

在之前的显示中有一些有趣的点，值得注意：

每个球员大约有 30 个数据点
有趣的是，球员数据追溯到 1871 年
每支队伍大约有 1,000 个数据点
美国联赛和国家联赛显然更受欢迎
一些数据点的范围令人惊讶：
- 打击数范围从 0 到 700
- 跑垒数（r）范围从 0 到 177
- 击球数（h）范围从 0 到 257

如果我们只绘制一些主导性数据点，我们可以看到如下内容：

plot(players$h, type="l", col="blue")
# doubles lines(players$X2b, type="l", col="yellow")
# home runs lines(players$hr, type="l", col="green")
# triples lines(players$X3b, type="l", col="red")
# Create a title with a red, bold/italic font title(main="Hits", font.main=4)#, xlab="Players", ylab="Hits and Home Runs")

统计图按此顺序显示，以避免较小的值被较大的值覆盖（例如，击球数（最大数字）首先显示，而三垒打（最小数字）最后显示在先前数字的上方）。

我们显示了球员随时间变化的击球类型，如下所示。

我觉得很有趣的是，三垒打如此少见。而且，由于数据的组织方式是按时间顺序的，三垒打的数量已经有一段时间在减少。也许现在更注重击球员争取本垒打，而不是击三垒打。

使用 `dplyr` 操作数据

R 的 dplyr 包被描述为一个提供数据操作语法的包。它具备你期望的所有数据框操作入口，整合在一个包中。我们将使用 dplyr 包处理本章前面使用的棒球选手统计数据。

我们读取球员数据并显示前几行：

players <- read.csv(file="Documents/baseball.csv", header=TRUE, sep=",") head(players)

我们将使用 dplyr 包，因此我们需要将该包导入我们的笔记本：

library(dplyr)

将数据框转换为 `dplyr` 表

dplyr 包有函数可以将数据对象转换为 dplyr 表。dplyr 表以紧凑的格式存储数据，使用的内存远小于普通数据框。大多数其他 dplyr 函数也可以直接在该表上操作。

我们可以通过以下方式将数据框转换为表：

playerst <- tbl_df(players) playerst

这将导致一个非常相似的显示模式：

快速了解数据值范围

dplyr 还有一个可用的函数是 glimpse()。它会显示每一列的值范围。我们可以按以下方式使用该函数：

glimpse(playerst)

这有如下显示效果：

我喜欢这个方法，除了 summary/head 的显示方式，它能让你感知涉及的变量，而不是涉及的行，实际上是倒转了数据集。

抽样数据集

dplyr 包有一个函数可以从数据集中抽取样本，sample()。你传入数据集和你想要抽取的样本数量，sample_n()，以及样本比例，sample_frac()，如下例所示：

data <- sample_n(players, 30) glimpse(data)

我们可以看到如下截图中的结果：

注意，结果集中有 30 个观察值，如请求的那样。

在数据框中筛选行

我们可以使用的另一个函数是filter函数。filter函数以数据框作为参数并使用过滤语句。该函数在数据框的每一行上运行，并返回满足过滤语句的行：

#filter only players with over 200 hits in a season over200 <- filter(players, h > 200) head(over200) nrow(over200)

看起来很多球员每赛季都能达到 200 次击球。如果我们看看那些在赛季中还能超过 40 次本垒打的球员呢？

over200and40hr <- filter(players, h > 200 & hr > 40) head(over200and40hr) nrow(over200and40hr)

这是一个非常小的清单。我知道球员的名字有些混乱，但你可以认出一些，比如贝贝·卢斯（Babe Ruth）。

我想知道有没有球员在一个赛季中击中超过 300 次。

filter(players, h > 300)

有趣的是，没有记录符合我们的filter条件，但结果处理程序需要一些列，并抛出错误，因为在这种情况下没有列。通常，R 中的错误是由于编程错误造成的。R 很少因我认为正常的无结果数据而生成错误。

向数据框添加列

mutate函数可以使用您在其他地方看到过的熟悉的 R 编程语法向数据框添加列。在本例中，我们正在向数据框添加一个列，该列显示球员上场击中的百分比：

pct <- mutate(players, hitpct = h / ab) head(pct)

你可以在上面的显示中看到右侧的新列。如果我们运行summary，我们将得到所有字段（包括新的hitpct）的汇总：

summary(pct)

最大值为1.0，这意味着一些球员每次上场都能击中。同样地，那些值为0.0的球员从未击中。看起来在 20%范围内有一个狭窄的区间。假设1.0表示一次上场一次击中，因为所有其他值都是用多个小数点来衡量的。

在计算字段上获得汇总

我们可以使用summarize函数更直接地得到列的汇总值。该函数接受一个数据框和一个计算出的结果。我们可以通过以下脚本看到相同的结果：

summarize(pct, mean(hitpct, na.rm = TRUE))

在函数之间进行数据传输

我们可以通过piping数据在函数之间传递来获得相同的结果。在 R 编程中，使用符号%>%来表示管道。它是从magrittr包中获取的。管道符号通常被认为是英语中then的同义词。例如，R 语句data %>% function()意味着获取数据对象，然后将其传递给function()，就像我们输入了语句function(data)一样。

要使用管道对计算字段进行相同的汇总，我们将编写以下内容（意思是获取pct数据集，然后将其管道传递到summarize函数，并获取hitpct字段并将其管道传递给均值函数）：

library(magrittr) pct %>% summarize(hitpct %>% mean(na.rm = TRUE))

获得 99%分位数

我们可以使用quantile()函数查看 99%分界点的值。使用相同的样本数据，我们可以使用：

quantile(pct$hitpct, probs = 0.99, na.rm = TRUE)

这将产生相应的输出：

99%: 0.470588235294118

因此，47%的击球率是数据中 99%水平的分界点。考虑到三分之四的百分位数是 28%（如前面的hitpct图示），最后四分之一的数据点的表现差异很大——也就是说，确实有一些优秀的棒球选手。

我们可以使用以下方法列出那些在击球率前 25%的球员：

top_players <- filter(pct, hitpct > 0.47) top_players <- top_players[order(top_players$hitpct) , ] head(top_players) nrow(top_players) 198

如果球员按击球率降序排列，则会显示出击球率完美的球员，但他们的打击次数都不到 10 次。

我们可以如下查看数据点：

所以，我们在数据集中排名前 25%的有 200 个（198）球员，这意味着 1%的球员位于击球表现的前 25%。我没想到数据会如此失衡。

获取分组数据的摘要

好的，前面的步骤告诉我们一些关于单个球员的信息。总有一种说法，X 队总是比其他队更强。那么如果我们能按团队获得击球率并进行比较呢？

在这个例子中，我们按团队对球员进行分组，然后计算整个团队的平均击球率：

teamhitpct <- summarize(group_by(pct, team), mean(hitpct, na.rm = TRUE)) names(teamhitpct) <- c("team", "hitpct") summary(teamhitpct)

那么，历史上最强的队伍是谁呢？我们可以通过队伍的击球率对数据进行排序（-teamhitpct子句意味着结果应该按降序排列），方法如下：

teamhitpct <- teamhitpct[order(-teamhitpct$hitpct) , ] head(teamhitpct)

我不确定 CNU 队是什么，我猜它是一个早期的芝加哥队。其他队伍可以识别为费城、多伦多和波士顿。

我们可以使用之前为球员击球表现使用过的quantile函数，找出排名前几的队伍：

quantile(teamhitpct$hitpct, probs = 0.99)

这将给我们以下结果：

99%: 0.340577141193618

与前面的表格相比，我们可以看到只有两个队伍（在 130 个队伍中）位于顶尖表现组（大约是 1%的水平）。

使用 tidyr 整理数据

tidyr包可用于清理/整理你的数据集。tidyr的作用是重新排列你的数据，使其：

每一列是一个变量
每一行是一个观察值

当你的数据以这种方式排列时，分析起来会容易得多。有很多数据集发布时混合了列和行以及数值。如果你在实际应用中使用这些数据，你就必须根据需要进行调整。

tidyr提供了三个函数来清理你的数据：

gather
separate
spread

gather()函数将你的数据整理成键值对，类似于 Hadoop 数据库模型。我们使用股票价格在某个日期的标准例子，数据如下：

library(tidyr)
stocks <- data_frame(
 time = as.Date('2017-08-05') + 0:9, X = rnorm(10, 20, 1), #how many numbers, mean, std dev Y = rnorm(10, 20, 2), Z = rnorm(10, 20, 4) )

这将生成如下的数据：

每一行都有一个时间戳和当时三个股票的价格。

我们首先使用gather()将股票的键值对拆分出来。gather()函数与其将处理的数据框、输出列名以及需要忽略的列（-time）一起调用。因此，我们可以使用以下方法获得一个包含不同时间、股票和价格的行：

stocksg <- gather(stocks, stock, price, -time) head(stocksg)

这将生成如下的head()显示：

separate()函数用于拆分同一条目点中的值。

我们将使用来自 UCI 的道琼斯工业价格历史数据（archive.ics.uci.edu/ml/datasets/Dow+Jones+Index）：

dji <- read.csv("Documents/dow_jones_index.data") dj <- dji[,c("stock","date","close")] summary(dj) head(dj)

我们只关注stock、date和close列。

因此，我们已经收集了日期数据以开始（如果我们一开始有未整理的数据，我们会使用gather函数将其整理到此为止）。

spread()函数将从gather()函数得到的键值对拆分成多个列。我们使用包含源数据框的日期、作为列的值以及每个日期/列的数据点来调用spread()。继续我们的示例，我们可以使用以下方法按日期展开所有证券：

prices <- dj %>% spread(stock, close) summary(prices) head(prices)

这会导致如下的summary显示（缩短为仅显示前几个证券）：

它还会产生如下的head显示，按日期显示所有证券的价格：

我们还可以通过列出每行股票的所有价格来以类似的方式重新组织数据，使用：

times <- dj %>% spread(date, close) summary(times) head(times)

在这里，我们的行驱动项是股票，列标题是日期值，而每个股票/日期的数据点是该日期证券的收盘价。

从summary（简化版）中，我们可以看到以下内容：

以下是显示我们所需格式的head示例：

摘要

在本章中，我们读取了 CSV 文件并进行了快速数据分析，包括可视化以帮助理解数据。接下来，我们考虑了dplyr包中的一些函数，包括绘制数据项范围的概览、对数据集进行采样、筛选数据、使用 mutate 添加列以及生成摘要。在此过程中，我们还开始使用管道操作，便于将一个操作的结果传递到另一个操作中。最后，我们研究了tidyr包，通过使用相应的 gather、separate 和 spread 函数，将数据整理成独立的列和观测值。

在下一章中，我们将查看如何在 Jupyter 下生成仪表板。

第七章：Jupyter 仪表盘

仪表盘是一种将多个显示组合在一起用于演示的机制。例如，您可以将来自不同数据的 3 个生产线图形显示在一个屏幕上，放在同一个框架或仪表盘中。通过 Jupyter，我们可以使用许多机制来检索和可视化数据，然后将它们组合成一个演示。

可视化已准备好的字形数据

字形是一个符号。在这一部分，我们希望在图表的不同位置显示字形，而不是标准的点，因为字形应该为观众提供更多的视觉信息。通常，数据点的某些属性可以用来将其转化为有用的字形，正如我们将在以下示例中看到的那样。

ggplot2包在以多种方式可视化数据时非常有用。ggplot被描述为 R 的一个绘图系统。我们将看一个示例，展示全球范围内的火山数据点。我使用了来自www.ngdc.noaa.gov/nndc的国家环境信息中心的数据。我选择了 1964 年后的火山信息。

这生成了一组数据，我将其复制到本地 CSV 文件中：

#read in the CSV file as available as described previously
volcanoes = read.csv("volcanoes.csv")
head(volcanoes)

如果我们仅在世界地图上绘制点，就能看到火山的位置。我们使用mapdata包来生成地图：

最新的 R 语法要求在install.packages()函数调用中指定仓库（或镜像）的地址，以找到包的位置。

#install.packages("mapdata",repos = "http://cran.us.r-project.org")
library(maps)
map("world")
points(volcanoes$Longitude, volcanoes$Latitude, pch=16, col="red", cex=1)

有一些零散的火山，但主要的模式出现在西美洲和东南亚，正好与地壳板块线相对应。

我们可以将一些数据点用作字形，并使用ggplot2中的qplot（快速绘图）函数以类似的图形展示这些信息：

    library(ggplot2)
    qplot(Longitude, Latitude, data=volcanoes, color=Type, size=Elev, alpha=0.5)

这些参数是：

x，y，经度和纬度
我们参考的数据集
我们使用颜色来表示火山类型
数据点的大小与火山的海拔高度相关
设置 alpha 为中间值将绘制半阴影颜色点，以便我们更好地看到重叠部分

可能有一种方法可以将图形数据叠加在前面的地理地图上。

大多数火山位于海拔超过 2000 英尺的地方。东南亚和南太平洋的火山往往位于水下。

同样，我们可以展示标准的iris数据集中的字形数据：

library(ggplot2) 
qplot(Sepal.Length, Petal.Length, data=iris, color=Species, size=Petal.Length)

在这个图形中，我们使用了花瓣长度和物种。图形清楚地表明，物种可以通过植物部件的尺寸来确定。

或者，使用内置的钻石数据集，我们可以从字形数据中获取一些信息：

library(ggplot2) 
dsmall <- diamonds[sample(nrow(diamonds), 100), ]

(d <- ggplot(dsmall, aes(carat, price)) + geom_point(aes(shape = cut))) +
scale_shape(solid = FALSE)

我没想到尺寸（克拉）会是价格的主要决定因素。我知道每个人都听说过切工、净度和颜色会有所不同，但实际数据的图表似乎并没有显示这些因素有任何区别。

我们也可以用类似的方式查看汽车每加仑英里数的数据集：

library(ggplot2) 
ggplot(data=mpg) + 
 geom_point(mapping = aes(x=displ, y=hwy, color=class, size=class))

在这个例子中，我们将美学数据点映射，其中x、y位置基于发动机排量和高速公路每加仑英里数，但我们还根据车辆类别对数据点进行颜色和大小的标注。

这些符号使其更清晰：

无论发动机排量大小如何，SUV 的油耗较低
小型和紧凑型车辆的油耗更好
有趣的是，发动机尺寸在任何车辆类别中对油耗几乎没有影响

我知道这些是直观的数据点，但使用符号确认信息后，变得更加清晰。

发布笔记本

你可以使用 markdown 发布笔记本/仪表盘。Markdown 涉及向笔记本中的单元格添加注释，Jupyter 会将其解释并转换为你在其他已发布材料中看到的更标准的 HTML 表示形式。

在所有情况下，我们创建带有 markdown 类型的单元格。然后，我们在单元格中输入 markdown 的语法。一旦运行 markdown 单元格，单元格的显示会变成有效的 markdown 表示形式。你还应该注意，markdown 单元格没有行号标注，因为 markdown 单元格中没有代码执行。

字体 markdown

你可以使用 斜体 和 粗体 HTML 标记调整字体样式。例如，如果我们有如下格式的单元格代码，你可以使用带有斜体（<i>）和粗体（<b>）标签的 markdown：

当我们运行单元格时，我们看到有效的 markdown 如下：

列表 markdown

我们可以使用以下列表，其中我们开始一个无序列表（我们本可以使用 nl 来创建有序列表），并用列表项（<li>）标签将两个列表项括起来：

当我们运行这个单元格时，结果的 markdown 会显示如下：

标题 markdown

有一组层次化的标题可以使用，就像这个 markdown 示例一样，Jupyter 特别关注使用井号符号 # 的标题标记：

奇怪的是，Jupyter 决定特别处理这个 markdown，而不是其他的 markdown。

结果显示如下：

表格 markdown

我们可以像在 HTML 中一样生成表格。HTML 表格以 table 标签开始。每一行的表格以 tr 标签开始。每行中的列以 td 标签开始。或者，你可以使用列标题标签 th 来代替 td 标签。在这个示例中，我们有一个简单的表格 markdown：

显示结果：

表格默认在 markdown 中居中。我们可以使用更多的 markdown 将表格改为左对齐。我们还可以通过向 table 标签添加 'borders=1' 来为表格单元格加上边框。注意，列标题默认是加粗的。

代码 markdown

有时在展示中展示实际的代码语句是很有用的。我们通过在代码序列前添加三个反引号和语言名称来解释代码的含义。通过在末尾加上三个反引号来关闭代码，如以下示例所示：

显示的结果对应于语言特定的关键字高亮：

创建 Shiny 仪表盘

Shiny 是一个用于 R 的 Web 应用框架。它不要求用户编写 HTML 代码。通常有两个代码集：服务器端和 用户界面（UI）。这两个代码集都运行在 Shiny 服务器上。Shiny 服务器可以托管在你的机器上，或者通过多个托管公司托管在云端。

Shiny 服务器代码集处理数据访问、计算结果、获取用户输入和与其他服务器代码集交互以更改结果。UI 代码集处理展示的布局。

例如，如果你有一个应用程序生成数据的直方图，服务器集会获取数据并生成结果，展示结果，并与用户交互以更改结果——例如，可能会更改显示的数据桶的数量或范围。UI 代码集将严格处理布局。

Shiny 代码不能在 Jupyter 中运行。你可以使用 RStudio 开发代码。RStudio 是一个 集成开发环境（IDE），用于开发 R 脚本。

使用 RStudio，你可以开发 server.R 和 ui.R 组件的代码集，然后运行它们。运行一个 Shiny 应用程序将：

打开一个新的浏览器窗口
使用代码生成相应的 HTML 代码
在新的浏览器窗口中显示生成的 HTML

一旦你的应用程序工作正常，你也可以进一步将应用程序部署到互联网上的 R 服务器上，在那里你可以与其他用户分享与你的应用程序相关的 URL。

R 应用编码

当你启动 RStudio 时，你会看到（希望）熟悉的显示界面。显示界面中有四个主要窗口（从左上角开始顺时针方向）：

源代码窗口。在此窗口中，你可以输入特定 R 脚本文件的代码。在源代码窗口中保存文件会更新存储在磁盘上的文件（并且该文件可以在文件窗口中查看）。
环境/历史窗口。在这个窗口中，任何已加载的数据集或已建立的变量都会显示出来。
文件窗口。在此窗口中，当前的文件集会显示出来（并且可以通过点击文件将其拉入源代码窗口）。
执行窗口。在此窗口中，会向 R 服务器发送实际的命令。

对于简单的例子，我们将使用 server.R 和 io.R 文件，这些文件与 Shiny 网页上的文件非常相似。

app.R 文件包含以下命令：

## app.R ##
# load the shiny dashboard library
library(shinydashboard)
# create a dashboard with title, sidebar and body
ui <- dashboardPage(
 dashboardHeader(title = "Shiny dashboard"),
 dashboardSidebar(),
 dashboardBody(
 fluidRow(
 box(plotOutput("plot1", height = 250)), 
 box(
 title = "Controls",
 sliderInput("slider", "Number of observations:", 1, 100, 50)
 )
 )
 )
)
# when the server is called upon by system, will generate normalized data and plot it
# further, every interaction with the user will rerun this set of code generating new data
server <- function(input, output) {
 set.seed(122)
 histdata <- rnorm(500)
 output$plot1 <- renderPlot({
 data <- histdata[seq_len(input$slider)]
 hist(data)
 })
}
# start the application - this calls upon the browser to get output for display
shinyApp(ui, server)

UI 文件仅设置页面的布局：

library(shinydashboard)
# only calls upon shiny to create a dashboard page with header, sidebar and body
# each of these components is defined in the preceding app.R code.
dashboardPage(
 dashboardHeader(),
 dashboardSidebar(),
 dashboardBody()
)

如果我们在 RStudio 中运行 app.R 程序，我们将得到一个如下所示的页面。页面有一个头部区域、侧边栏和主体部分。主体包含显示区域和滑块控制。我们可以通过在源代码窗口中选择 app.R 文件，然后点击源代码窗口顶部的“运行应用”按钮来运行应用程序。点击“运行应用”按钮的效果是向控制台窗口输入命令 runApp('shiny')（所以你本来也可以自己输入该命令）。一旦应用程序运行，你会在控制台窗口看到类似于“监听 http://127.0.0.1:6142”的消息，这意味着 R 决定将你的应用程序部署到你的机器上。

在新打开的浏览器窗口中，我们可以看到展示的 Shiny 应用程序：

请注意，Shiny 显示是响应式的。响应式网页应用程序会根据窗口物理尺寸的变化调整组件的布局。正如之前的展示一样，你会看到组件按出现的顺序水平排列。当我调整显示的尺寸以适应打印页面时，布局重新组织了。

页面是完全动态的。当你移动滑块时，图形会重新生成。

发布你的仪表盘

一旦你对你的仪表盘感到满意，就可以将其发布到互联网上。源代码窗口右上角附近的“部署”按钮将把你的应用程序集推送到你所连接的服务。

在我的例子中，我使用了一个免费的账户，注册于www.shinyapps.io/。还有其他服务可以托管你的 Shiny 应用程序。

一旦你选择发布应用程序，RStudio 将提示你输入发布所需的服务和凭据（凭据会在你注册时提供给你）。

之后，文件会被上传到你的主机，应用程序也开始运行。通过访问我的演示账户 dantoomeysoftware.shinyapps.io/shiny，你会看到与前面展示完全相同的界面，只是我们是在托管的服务器上运行。

所以，我们有一个易于使用的系统，可以使用各种显示机制创建我们的仪表盘。总体而言，比我预期的要容易得多——考虑到我在其他系统（如微软工具）中开发此类演示的经验。

构建独立的仪表盘

使用 Node.js，开发人员提出了一种在 jupyter-dashboard-server 上托管你的仪表盘/笔记本而不使用 Jupyter 的方法。

安装需要安装 Node.js（因为服务器是用 Node.js 编写的）。这是一个较大的安装包。

一旦安装了 Node.js，其中一个安装的工具是 npm-node 产品管理器。你可以使用 npm 通过以下命令安装仪表盘服务器：

npm install -g jupyter-dashboards-server

安装完成后，你可以使用以下命令运行服务器：

C:\Users\Dan>jupyter-dashboards-server --KERNEL_GATEWAY_URL=http://my.gateway.com

mygateway.com 是一个占位符。你需要使用你的网关服务器（如果需要）。此时，服务器已在你提到的环境中运行，并将输出几行信息：

Using generated SESSION_SECRET_TOKEN
Jupyter dashboard server listening on 127.0.0.1:3000

你可以在浏览器中打开 URL (http://127.0.0.1:3000/dashboards)，查看服务器控制台的样子：

至于开发你可以托管在服务器上的仪表盘，我们还需要安装更多内容：

conda install jupyter_dashboards -c conda-forge

然后启用扩展（这是一个笔记本扩展）：

jupyter nbextension enable jupyter_dashboards --py --sys-prefix

然后需要使用以下命令安装布局扩展：

pip install jupyter_dashboards_bundlers
jupyter bundlerextension enable --sys-prefix --py dashboards_bundlers

此时，你可以使用以下命令上传你的仪表盘笔记本文件：

POST /_api/notebooks/[PATH/]NAME

URL 前缀是你使用的托管站点，PATH 是可选的（默认情况下指向根位置），NAME 是你确定的名称。

小结

在本章中，我们使用符号图形化地展示数据，强调数据的关键方面。我们使用 Markdown 注释了笔记本页面。我们使用 Shiny 生成了一个交互式应用程序。并且我们看到了在 Jupyter 之外托管笔记本的方法。

在下一章中，我们将介绍 Jupyter 下的统计建模。

第八章：统计建模

在本章中，我们将使用原始数据并通过构建统计模型来尝试解读信息。一旦我们构建了模型，通常会更容易看到共同点，我们也可以确定趋势。

将 JSON 转换为 CSV

对于本章，我们将使用来自 www.yelp.com/dataset/challenge 的 Yelp 数据集。此部分使用的是挑战的 第 9 轮 数据集。背景是，Yelp 是一个用于评估不同产品和服务的站点，Yelp 会将评估结果发布给用户。

数据集文件是一个非常大的（几 GB）评级数据。下载的内容包括多个评级信息集——包括商业评级、评论、建议（比如这会是一个不错的地方）和用户集。我们对评论数据感兴趣。

在处理如此大的文件时，可能需要找到并使用一个大型文件编辑器，以便能够深入查看数据文件。在 Windows 上，大多数标准编辑器限制为几个兆字节。我使用了 Large Text File Viewer 程序来打开这些 JSON 文件。

所有文件都是 JSON 格式的。JSON 是一种人类可读的格式，包含结构化元素——例如，一个城市对象包含街道对象。虽然 JSON 格式易于阅读，但在处理大量元素时，这种格式显得笨拙。在 reviews 文件中有几百万行。因此，我们首先使用此脚本将 JSON 转换为平坦的 CSV 格式，以便更容易进行处理：

import time import datetime import json, csv print( datetime.datetime.now().time()) headers = True #with open('c:/Users/Dan/reviews.json') as jsonf, open('c:/Users/Dan/reviews.csv', "wb") as csvf: filein = 'c:/Users/Dan/yelp_academic_dataset_review.json' fileout = 'c:/Users/Dan/yelp_academic_dataset_review.csv' with open(filein) as jsonf, open(fileout, "wb") as csvf:
 for line in jsonf: data = json.loads(line)        #remove the review text
 data.pop('text') if headers: w = csv.DictWriter(csvf, data.keys()) w.writeheader() headers = False w.writerow(data)print( datetime.datetime.now().time())

我正在打印出开始和结束时间，以便了解此过程的耗时。对于我的机器，转换文件花费了 1.5 分钟。我在使前述代码以令人满意的速度工作之前，尝试了几个版本的代码。在开发这个脚本时，我取了原始数据文件的一个小子集（2000 行），并一直在这个文件上工作，直到进展足够为止。

如你所见，我正在读取来自 Yelp 的原始 JSON 文件，并写出一个相应的 CSV 文件。

该脚本读取每一行 JSON（每一行包含一个完整的对象），并输出相应的 CSV。我去除了评论文本，因为我并未对评论文本进行评估，而且评论文本占用了大量空间。使用此编码后，评论文件的大小从 3GB 降低到了 300MB。除此之外，我们确保将头信息作为第一条记录写入 CSV。我然后使用一个单独的脚本/笔记本条目读取 CSV 并进行处理。

评估 Yelp 评论

我们使用这个脚本读取处理后的 Yelp 评论，并打印出一些数据统计信息：

reviews <- read.csv("c:/Users/Dan/yelp_academic_dataset_review.csv")

我通常在加载数据后查看一些数据，以便直观检查是否一切按预期工作。我们可以通过调用 head() 函数来实现：

head(reviews)

数据总结

所有列似乎都正确加载了。现在，我们可以查看数据的总结统计信息：

summary(reviews)

总结中有几点值得注意：

我原本以为一些数据点会是TRUE/FALSE，0/1，但实际上它们是区间值；例如，funny的最大值超过 600；useful的最大值为 1100，cool的最大值为 500。
所有的 ID（用户、公司）都被篡改了。我们可以使用用户文件和公司文件来找到确切的引用。
星级评分是1到5，这符合预期。然而，平均值和中位数大约为4，我认为这表明很多人只愿意花时间写好评。

评论分布

我们可以通过简单的直方图来了解评论的分布，使用如下代码：

hist(reviews$stars)

这将生成一个内联直方图，如下所示：

再次，我们看到人们倾向于只写好评（但差评的数量还是相当大的）。

寻找评分最高的公司

我们可以使用以下脚本，找到评分最多的公司（5 星）。这个脚本使用类似 SQL 的数据框访问方式，SQL 有内置的机制来按需搜索、选择和排序数据组件。

在这个脚本中，我们正在构建一个计算数据框，包含两列：business_id和评论count。该数据框按评分从高到低排列。创建完成后，我们展示数据框的头部，以便查看数据集中的评分最高的公司：

#businesses with most 5 star ratings
#install.packages("sqldf", repos='http://cran.us.r-project.org')
library(sqldf)
five_stars = sqldf("select business_id, count(*) from reviews where stars = 5 group by business_id order by 2 desc")
head(five_stars)

值得注意的是，排名前五的公司在评分数量上存在如此明显的偏差（排名第1的公司几乎是排名第6公司评分数量的两倍）。你会怀疑，是否在评分过程中存在某种串通现象，导致这种差异。再次提醒，名字到目前为止都被篡改了。

寻找评价最高的公司

所以，这些公司得到了最多的好评。那哪些公司得到了最多的评论呢？我们可以使用类似的脚本来确定最被评价的公司：

#which places have most ratings
library(sqldf)
most_ratings = sqldf("select business_id, count(*) from reviews group by business_id order by 2 desc")
head(most_ratings)

在此脚本中，我们不对评分进行限定，以确定它的归属，最终结果是：

同样地，我们可以看到少数几家公司有远超平均值的评分数。而且，虽然名字被篡改，但我们确实看到四家公司既是评分最高的公司，同时也出现在最被评价的公司名单中。

查找评分最高的公司所有的评分

如果我们查看其中一家评分最高的公司，看看评分分布情况如何呢？我们可以使用以下脚本：

# range of ratings for business with most ratings
library(sqldf)
most_rated = sqldf("select * from reviews where business_id = '4JNXUYY8wbaaDmk3BPzlWw' ")
hist(most_rated$stars)

该脚本获取一个评分最高的 ID，访问其所有评分，并显示一个相关的直方图：

所以，对于某些评分最高的公司，几乎没有低评分。

确定评分与评论数量之间的相关性

我们可以使用以下脚本查看星级评分与评论数量之间的相关性：

# correlation number of reviews and number of stars
library(sqldf)
reviews_stars = sqldf("select stars,count(*) as reviews from reviews group by stars")
reviews_stars
cor(reviews_stars)

所以，我们看到 5 星评价是 1 星评价的三倍多。我们还看到了评论数量和星级数量之间非常高的相关性（0.8632361）。人们只愿意为好的公司评分。这使得仅使用 Yelp 来确定公司是否被评论过变得非常有趣。如果公司没有评分（或者评分很少），那么那些未书面的评论就会很差。

我们可以使用以下脚本可视化公司评分和评论数量之间的关系：

#correlation business and rating
library(sqldf)
business_rating = sqldf("select business_id, avg(stars) as rating from reviews group by business_id order by 2 desc")
head(business_rating)
hist(business_rating$rating)

其中，business_rating 数据框是一个商业和平均星级评分的列表。生成的直方图如下：

这看起来像是一个泊松分布。有趣的是，公司的评级分布呈现出如此自然的离散性。

构建评论模型

我们可以从数据集中构建一个模型，估算一个评分可能涉及多少颗星。然而，评论中的数据点只有：

funny
useful
cool

这些似乎不是良好的评分指示器。我们可以使用一个模型，例如：

model <- lm(stars ~ funny + useful + cool, data=reviews)
summary(model)

这生成了模型的统计数据：

如预期的那样，我们没有足够的信息来进行处理：

超过四百万的自由度，几乎每个评论一个自由度
P 值非常小——我们估算正确的概率几乎为零
3.7 截距（接近范围的中点）
这么低的影响率（每个因素不到一次），意味着我们距离截距没有太大的变化。

使用 Python 比较评分

在之前的示例中，我们使用 R 来处理从转换后的 JSON 到 CSV 文件构建的数据框。如果我们使用 Yelp 的企业评分文件，可以直接使用 Python，因为它更小，且产生类似的结果。

在这个示例中，我们从 Yelp 文件中收集包含餐厅的商业类别的菜品。我们为所有菜品累计评分，然后为每个菜品生成平均分。

我们将 JSON 文件读取为单独的行，并将每行转换为 Python 对象：

我们将每行转换为 Unicode，使用 errors=ignore 选项。这是因为数据文件中存在许多错误字符。

import json
#filein = 'c:/Users/Dan/business.json'
filein = 'c:/Users/Dan/yelp_academic_dataset_business.json'
lines = list(open(filein))

我们为菜品的评分使用字典。字典的键是菜品的名称，字典的值是该菜品的评分列表：

ratings = {}
for line in lines:
 line = unicode(line, errors='ignore')
 obj = json.loads(line)
 if obj['categories'] == None:
 continue
 if 'Restaurants' in obj['categories']:
 rating = obj['stars']
 for category in obj['categories']:
 if category not in ratings:
 ratings[category] = []
 clist = ratings.get(category)
 clist.append(rating)

现在我们收集了所有评分后，可以生成一个新的菜品字典，包含平均评分。我们还会累计总数以生成整体平均值，并跟踪评分最高的菜品：

cuisines = {}
total = 0
cmax = ''
maxc = 0
for cuisine in ratings:
 clist = ratings[cuisine]
 if len(clist) < 10:
 continue
 avg = float(sum(clist))/len(clist)
 cuisines[cuisine] = avg
 total = total + avg
 if avg > maxc:
 maxc = avg
 cmax = cuisine

print ("Highest rated cuisine is ",cmax," at ",maxc)
print ("Average cuisine rating is ",total/len(ratings))

print (cuisines)

有趣的是，Personal Chefs（私人厨师）得分最高。我只听说过名人会有私人厨师，但数据显示这可能是值得的。所有菜系的平均评分为 1.6，极为低下。早些时候我们看到数据似乎没有高低评分的平衡。然而，仔细查看结果输出时，会发现很多条目并非菜系，尽管Restaurants（餐馆）键存在。我曾尝试通过仅计算评分数达到 10 条以上的菜系来剔除坏数据，虽然剔除了一部分坏数据，但仍然有许多错误记录存在。

按菜系可视化平均评分

现在我们已经计算出各菜系的平均值，可以通过直方图展示它们的分布情况。我们首先将字典转换为数据框。然后将数据框中的Rating列绘制为直方图：

我们使用五个箱子来对应五个可能的评分。

import pandas as pd
import numpy as np
df = pd.DataFrame(columns=['Cuisine', 'Rating'])
for cuisine in cuisines:
 df.loc[len(df)]=[cuisine, cuisines[cuisine]]
hist, bin_edges = np.histogram(df['Rating'], bins=range(5))

import matplotlib.pyplot as plt
plt.bar(bin_edges[:-1], hist, width = 1)
plt.xlim(min(bin_edges), max(bin_edges))
plt.show()

再次，我们看到明显偏向高平均值的趋势。我曾尝试改善数据展示的梯度，但未能成功。

随机评分搜索

由于数据已经以易于加载的格式存在，我们可以根据任意条件进行搜索，比如Personal Chefs（私人厨师）是否允许狗——也许他们会为你的狗定制烹饪。

我们可以使用如下脚本：for line in lines:

 line = unicode(line, errors='ignore')
 obj = json.loads(line)
 if obj['categories'] == None:
 continue
 if 'Personal Chefs' in obj['categories']:
 if obj['attributes'] == None:
 continue
 for attr in obj['attributes']:
 print (attr)

在这里我们对筛选出的项目做一些有用的处理。此脚本仅会显示Personal Chefs（私人厨师）的属性。可以在以下显示中看到：

我们也可以轻松地进行某些计算或其他操作，快速缩小范围，专注于数据的特定部分。

确定评分数量与评分之间的关系

根据前述结果，似乎人们大多数时间只是进行积极投票。我们可以查看公司收到的投票数量与其评分之间是否存在关系。

首先，我们使用以下脚本来累积数据集，提取每家公司的投票数和评分：

#determine relationship between number of reviews and star rating
import pandas as pd
from pandas import DataFrame as df 
import numpy as np 

dfr2 = pd.DataFrame(columns=['reviews', 'rating'])
mynparray = dfr2.values

for line in lines:
 line = unicode(line, errors='ignore')
 obj = json.loads(line)
 reviews = int(obj['review_count'])
 rating = float(obj['stars'])
 arow = [reviews,rating]
 mynparray = np.vstack((mynparray,arow)) 

dfr2 = df(mynparray)
print (len(dfr2))

这段代码仅构建了我们两个变量的数据框。我们使用 NumPy，因为它更容易将一行添加到数据框中。当我们处理完所有记录后，我们会将 NumPy 数据框转换回 pandas 数据框。

列名在转换过程中丢失了，因此我们将它们重新添加并绘制一些总结统计数据：

dfr2.columns = ['reviews', 'rating']
dfr2.describe()

在以下输出中，我们可以看到收集的评论和评分数据的布局。Yelp 没有对该数据集的数据输入做出任何限制。评分应有 5 个独特的值：

接下来，我们绘制数据以获得关系的视觉线索，使用以下内容：

#import matplotlib.pyplot as plt
dfr2.plot(kind='scatter', x='rating', y='reviews')
plt.show()

因此，最终数据似乎呈现出一个清晰的泊松分布，与之前的business_rating直方图相比。

接下来，我们计算回归参数：

#compute regression
import statsmodels.formula.api as smf

# create a fitted model in one line
lm = smf.ols(formula='rating ~ reviews', data=dfr2).fit()

# print the coefficients
lm.params

我们计算了所有评分值的截距。我原本预期会得到一个单一的值。

现在，我们使用以下方法来确定观察数据的范围：

#min, max observed values
X_new = pd.DataFrame({'reviews': [dfr2.reviews.min(), dfr2.reviews.max()]})
X_new.head()

所以，正如我们之前猜测的那样，一些企业有着大量的评论。

Now, we can make predictions based on the extent data points:
#make corresponding predictions
preds = lm.predict(X_new)
preds

我们看到的预测值范围远比预期的大。接下来，绘制出观察数据和预测数据：

# first, plot the observed data
dfr2.plot(kind='scatter', x='reviews', y='rating')

# then, plot the least squares line
plt.plot(X_new, preds, c='red', linewidth=2)
plt.show()

在以下展示的图中，似乎没有发现评论数量和企业评分之间的关系。这看起来像是一场数字游戏——如果你让人们评论你的企业，平均而言，他们会给你一个高分。

似乎没有发现评论数量和企业评分之间的关系。

总结

在这一章中，我们使用了数据并将 JSON 文件转换为 CSV 文件。我们评估了 Yelp 的餐饮评论数据集，确定了评分最高和评论最多的企业。我们查看了评分的分布。我们使用 Python 对 Yelp 企业评分进行了类似的评估，发现数据的分布非常相似。

在下一章中，我们将学习如何在 Jupyter 下使用机器学习。

第九章：使用 Jupyter 进行机器学习

在本章中，我们将使用几种机器学习算法在 Jupyter 下进行演示。我们在 R 和 Python 中编写代码，以展示 Jupyter 开发者可用的各种选择。

朴素贝叶斯

朴素贝叶斯是一种使用概率通过贝叶斯定理对数据进行分类的算法，前提是特征之间有较强的独立性。贝叶斯定理基于先验条件估计事件的概率。因此，整体上，我们使用一组特征值来估算一个值，假设在这些特征具有相似值时，条件保持一致。

使用 R 的朴素贝叶斯

我们的第一个朴素贝叶斯实现使用了 R 编程语言。该算法在 R 中的实现被编码在e1071库中。e1071似乎是该包开发时所在学校的部门标识符。

我们首先安装包，并加载库：

#install.packages("e1071", repos="http://cran.r-project.org") 
library(e1071) 
library(caret) 
set.seed(7317) 
data(iris)

关于这些步骤的一些备注：

install.packages调用被注释掉了，因为我们不想每次运行脚本时都运行它。
e1071是朴素贝叶斯算法包。
caret包包含一个方法，用于随机分割数据集。
我们设置了seed，以便能够重现结果。
我们在这个示例中使用的是iris数据集。具体来说，使用其他iris特征来预测物种。

包的调用方式如下：

model <- naiveBayes(response ~ ., data=training) 
prediction <- predict(model, test, type="class")

naiveBayes的参数如下：

公式形式为 y ~ x1 + x2 .... ——尝试根据x1, x2, ...预测y。
数据框
可选的拉普拉斯平滑
基于布尔过滤器的可选数据子集
可选的处理na值的函数（na.action）——默认是通过

一旦我们建立了模型，就可以使用predict()函数并传入参数来尝试预测：

模型（来自前面的调用）
数据框
输入数据是类别数据还是原始数据（条件）

因此，我们继续使用iris示例：

trainingIndices <- createDataPartition(iris$Species, p=0.75, list=FALSE) 
training <- iris[trainingIndices,] 
testing <- iris[-trainingIndices,] 
nrow(training) 
nrow(testing) 
114 
36

在这里，我们将数据拆分为 75%的训练集和 25%的测试集，从每个数据框中的行数可以看到这一点。

接下来，我们构建模型——我们正在尝试根据数据框的其他特征/列预测Species：

model <- naiveBayes(Species ~ ., data=training) 
model

有趣的是，Apriori 假设假设可能性之间是均等的。花萼长度、宽度和花瓣长度对物种有很强的影响。

我们基于模型对测试数据进行预测：

prediction <- predict(model, testing, type="class")

现在，我们需要衡量模型的准确性。通常我们可以使用散点图，x来自实际值，y来自预测值，但由于我们有分类数据，我们可以建立一个实际值与预测值的向量，并在一个新的结果数据框中比较这两个向量：

results <- data.frame(testing$Species, prediction) 
results["accurate"] <- results['testing.Species'] == results['prediction'] 
nrow(results) 
nrow(results[results$accurate == TRUE,]) 
36 
35

最终我们得到了一个提供 97%准确率的模型（35/36）。这是一个非常好的性能水平，几乎处于优秀的统计边界内（+/- 2%）。

使用 Python 的朴素贝叶斯

算法的 Python 实现位于 sklearn 库中。整个过程要简单得多。首先，加载 iris 数据集：

from sklearn import datasets 
irisb = datasets.load_iris() 
iris = irisb['data'] 
iris.shape

调用内建的高斯朴素贝叶斯估算器来一步完成模型和预测：

from sklearn.naive_bayes import GaussianNB 
gnb = GaussianNB() 
y_pred = gnb.fit(irisb.data, irisb.target).predict(irisb.data)

确定模型的准确度：

print("Number of errors out of a total %d points : %d"  
      % (irisb.data.shape[0],(irisb.target != y_pred).sum())) 
Number of errors out of a total 150 points : 6

我们得到了非常相似的估算准确性结果。

使用 R 实现最近邻

在这个例子中，我们使用来自 ics.edu 的房价数据。首先，我们加载数据并指定列名：

housing <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data") 
colnames(housing) <- c("CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO", "B", "LSTAT", "MDEV") 
summary(housing)

我们重新排序数据，使得关键值（房价 MDEV）按升序排列：

housing <- housing[order(housing$MDEV),]

现在，我们可以将数据拆分为训练集和测试集：

#install.packages("caret") 
library(caret) 
set.seed(5557) 
indices <- createDataPartition(housing$MDEV, p=0.75, list=FALSE) 
training <- housing[indices,] 
testing <- housing[-indices,] 
nrow(training) 
nrow(testing) 
381 
125

我们使用这两个数据集来构建最近邻模型：

library(class) 
knnModel <- knn(train=training, test=testing, cl=training$MDEV) 
knnModel 
10.5 9.7 7 6.3 13.1 16.3 16.1 13.3 13.3...

让我们看看结果：

plot(knnModel)

在数据的较高点附近存在轻微的泊松分布。我认为这符合自然数据的特点。起始和结束的尾部在页面上显得异常偏离。

那么，这个模型的准确度如何呢？我没有找到将 knnModel 中预测的因素转化为数值的清晰方法，所以我将它们提取到一个平面文件中，然后单独加载：

predicted <- read.table("housing-knn-predicted.csv") 
colnames(predicted) <- c("predicted") 
predicted

预测值
10.5
9.7
7.0

然后我们可以构建一个 results 数据框：

results <- data.frame(testing$MDEV, predicted)

计算我们的准确率：

results["accuracy"] <- results['testing.MDEV'] / results['predicted'] 
head(results) 
mean(results$accuracy) 
1.01794816307793

testing.MDEV	预测值	准确率
`5.6`	`10.5`	`0.5333333`
`7.2`	`9.7`	`0.7422680`
`8.1`	`7.0`	`1.1571429`
`8.5`	`6.3`	`1.3492063`
`10.5`	`13.1`	`0.8015267`
`10.8`	`16.3`	`0.6625767`

因此，我们在测试数据的 2% 内（1.01）进行估算。

使用 Python 实现最近邻

在 Python 中，我们有非常相似的步骤来生成最近邻估算。

首先，我们导入要使用的包：

from sklearn.neighbors import NearestNeighbors 
import numpy as np 
import pandas as pd

Numpy 和 pandas 是标准库。最近邻是 sklearn 的一项特性。

现在，我们加载我们的房价数据：

housing = pd.read_csv("http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data",
                       header=None, sep='\s+') 
housing.columns = ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", \ 
"DIS", "RAD", "TAX", "PRATIO", \ 
"B", "LSTAT", "MDEV"] 
housing.head(5)

这是我们之前在 R 中看到的相同数据。

让我们看看它的大小：

len(housing) 
506

然后将数据划分为 训练 和 测试 集：

mask = np.random.rand(len(housing)) < 0.8 
training = housing[mask] 
testing = housing[~mask] 
len(training) 
417 
len(testing) 
89

查找最近邻：

nbrs = NearestNeighbors().fit(housing)

显示它们的索引和距离。索引的变化较大。距离似乎呈现出带状分布：

distances, indices = nbrs.kneighbors(housing) 
indices 
array([[  0, 241,  62,  81,   6], 
       [  1,  47,  49,  87,   2], 
       [  2,  85,  87,  84,   5], 
       ...,  
       [503, 504, 219,  88, 217], 
       [504, 503, 219,  88, 217], 
       [505, 502, 504, 503,  91]], dtype=int32) 
distances 
array([[  0\.        ,  16.5628085 , 17.09498324,18.40127391, 
         19.10555821], 
       [  0\.        ,  16.18433277, 20.59837827, 22.95753545, 
         23.05885288] 
       [  0\.        ,  11.44014392, 15.34074743, 19.2322435 , 
         21.73264817], 
       ...,  
       [  0\.        ,   4.38093898,  9.44318468, 10.79865973, 
         11.95458848], 
       [  0\.        ,   4.38093898,  8.88725757, 10.88003717, 
         11.15236419], 
       [  0\.        ,   9.69512304, 13.73766871, 15.93946676, 
         15.94577477]])

从 训练 集构建一个最近邻模型：

from sklearn.neighbors import KNeighborsRegressor 
knn = KNeighborsRegressor(n_neighbors=5) 
x_columns = ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO", "B", "LSTAT"] 
y_column = ["MDEV"] 
knn.fit(training[x_columns], training[y_column]) 
KNeighborsRegressor(algorithm='auto', leaf_size=30, metric='minkowski', 
          metric_params=None, n_jobs=1, n_neighbors=5, p=2, 
          weights='uniform')

有趣的是，使用 Python 时我们不需要将模型单独存储。方法是有状态的。

做出我们的预测：

predictions = knn.predict(testing[x_columns]) predictions array([[ 20.62],
 [ 21.18], [ 23.96], [ 17.14], [ 17.24], [ 18.68], [ 28.88],

确定我们预测房价的准确性：

columns = ["testing","prediction","diff"] 
index = range(len(testing)) 
results = pd.DataFrame(index=index, columns=columns) 

results['prediction'] = predictions 

results = results.reset_index(drop=True) 
testing = testing.reset_index(drop=True) 
results['testing'] = testing["MDEV"] 

results['diff'] = results['testing'] - results['prediction'] 
results['pct'] = results['diff'] / results['testing'] 
results.mean() 
testing       22.159551 
prediction    22.931011 
diff          -0.771461 
pct           -0.099104

我们的平均差异为¾，而平均值为22。这应该意味着大约 3%的平均百分比差异，但计算出来的平均百分比差异接近 10%。因此，我们在 Python 中并没有做出很好的估计。

决策树

在本节中，我们将使用决策树来预测值。决策树具有逻辑流程，用户根据属性做出决策，沿着树形结构向下直到根节点，在根节点提供分类结果。

在这个示例中，我们使用汽车特征，例如车辆重量，来判断该车是否能提供良好的油耗。这些信息来自页面 alliance.seas.upenn.edu/~cis520/wiki/index.php?n=Lectures.DecisionTrees。我将数据复制到 Excel 中，并将其作为 CSV 文件用于本示例。

R 中的决策树

我们加载用于rpart和caret的库。rpart包含决策树建模包，caret包含数据划分功能：

library(rpart) 
library(caret) 
set.seed(3277)

我们加载mpg数据集，并将其分为训练集和测试集：

carmpg <- read.csv("car-mpg.csv") 
indices <- createDataPartition(carmpg$mpg, p=0.75, list=FALSE) 
training <- carmpg[indices,] 
testing <- carmpg[-indices,] 
nrow(training) 
nrow(testing) 
33 
9

我们开发一个模型，根据其他因素预测mpg的可接受性：

fit <- rpart(mpg ~ cylinders + displacement + horsepower + weight + acceleration +
             modelyear + maker, method="anova", data=training) 
fit 
n= 33  

node), split, n, deviance, yval 
      * denotes terminal node 

1) root 33 26.727270 1.909091   
2) weight>=3121.5 10  0.000000 1.000000 * 
3) weight< 3121.5 23 14.869570 2.304348   
6) modelyear>=78.5 9  4.888889 1.888889 * 
7) modelyear< 78.5 14  7.428571 2.571429 *

显示的是决策树的文本表示形式。你可以通过以下方式图形化查看决策树：

plot(fit) 
text(fit, use.n=TRUE, all=TRUE, cex=.5)

这似乎是一个非常简单的模型。1980 年可能发生了油耗的变化，因为这是决策树的主要驱动因素。

最后，我们预测值并将其与我们的testing集进行比较：

predicted <- predict(fit, newdata=testing) 
predicted 
testing

看起来这个包已经将Bad、OK和Good转换为数字等效值，其中1代表Bad，其他则为OK或Good。总体来说，我们不确定这个模型是否有效。显然可用的数据量很小，增加测试集的大小可能会澄清模型。

Python 中的决策树

我们可以在 Python 中执行相同的分析。加载一些需要用到的库：

import pandas as pd 
import numpy as np 
from os import system 
import graphviz #pip install graphviz  
from sklearn.cross_validation import train_test_split 
from sklearn.tree import DecisionTreeClassifier 
from sklearn.metrics import accuracy_score 
from sklearn import tree

读取mpg数据文件：

carmpg = pd.read_csv("car-mpg.csv") 
carmpg.head(5)

将数据分解为因素和结果：

columns = carmpg.columns 
mask = np.ones(columns.shape, dtype=bool) 
i = 0 #The specified column that you don't want to show 
mask[i] = 0 
mask[7] = 0 #maker is a string 
X = carmpg[columns[mask]] 
Y = carmpg["mpg"]

将数据分为训练集和测试集：

X_train, X_test, y_train, y_test  
= train_test_split( X, Y, test_size = 0.3,  
random_state = 100)

创建一个决策树模型：

clf_gini = tree.DecisionTreeClassifier(criterion = "gini",  
random_state = 100, max_depth=3, min_samples_leaf=5)

计算模型的拟合度：

clf_gini.fit(X_train, y_train) 
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=3, 
            max_features=None, max_leaf_nodes=None, 
            min_impurity_split=1e-07, min_samples_leaf=5, 
            min_samples_split=2, min_weight_fraction_leaf=0.0, 
            presort=False, random_state=100, splitter='best')

绘制出决策树：

#I could not get this to work on a Windows machine 
#dot_data = tree.export_graphviz(clf_gini, out_file=None,  
#                         filled=True, rounded=True,   
#                         special_characters=True)   
#graph = graphviz.Source(dot_data)   
#graph

神经网络

我们可以将房屋数据建模为神经网络，其中不同的数据元素作为系统的输入，而网络的输出是房价。通过神经网络，我们最终得到一个图形化的模型，展示了每个输入应应用的因素，以得出房价。

R 中的神经网络

R 中有一个神经网络包。我们加载它：

#install.packages('neuralnet', repos="http://cran.r-project.org") 
library("neuralnet")

加载房屋数据：

filename = "http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data" 
housing <- read.table(filename) 
colnames(housing) <- c("CRIM", "ZN", "INDUS", "CHAS", "NOX",  
                       "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO", 
                       "B", "LSTAT", "MDEV")

将房屋数据分为训练集和测试集（我们在之前的示例中已经看到过这种编码）：

housing <- housing[order(housing$MDEV),] 
#install.packages("caret") 
library(caret) 
set.seed(5557) 
indices <- createDataPartition(housing$MDEV, p=0.75, list=FALSE) 
training <- housing[indices,] 
testing <- housing[-indices,] 
nrow(training) 
nrow(testing) 
testing$MDEV

计算我们的neuralnet模型：

nnet <- neuralnet(MDEV ~ CRIM + ZN + INDUS + CHAS + NOX  
                  + RM + AGE + DIS + RAD + TAX + PRATIO  
                  + B + LSTAT, 
                  training, hidden=10, threshold=0.01) 
nnet

neuralnet模型的显示信息相当广泛。以下是第一组显示内容。尚不清楚这些点是否有用：

$call 
neuralnet(formula = MDEV ~ CRIM + ZN + INDUS + CHAS + NOX + RM +  
    AGE + DIS + RAD + TAX + PRATIO + B + LSTAT, data = training,  
    hidden = 10, threshold = 0.01) 

$response 
    MDEV 
399  5.0 
406  5.0 
... 
$covariate 
           [,1]  [,2]  [,3] [,4]   [,5]  [,6]  [,7]    [,8] [,9] [,10] [,11] 
  [1,] 38.35180   0.0 18.10    0 0.6930 5.453 100.0  1.4896   24   666  20.2 
  [2,] 67.92080   0.0 18.10    0 0.6930 5.683 100.0  1.4254   24   666  20.2 
  [3,]  9.91655   0.0 18.10    0 0.6930 5.852  77.8  1.5004   24   666  20.2 
....

显示模型：

plot(nnet, rep="best")

这只是图表的上半部分。如您所见，每个因素都被调整到模型中以得出我们的房价。这并不有用——每个因素不可能都这么重要。

确定我们在该模型中的准确性：

results <- compute(nnet, testing[,-14]) 
diff <- results$net.result - testing$MDEV 
sum( (diff - mean(diff) )² ) #sum of squares 
9275.74672

鉴于该模型似乎非常不准确，我不确定在 Python 中按照相同的步骤操作是否会有益。

随机森林

随机森林算法尝试多棵随机决策树，并提供在所用参数下最有效的树。

R 中的随机森林

在 R 中，我们包括了将要使用的包：

install.packages("randomForest", repos="http://cran.r-project.org") 
library(randomForest)

加载数据：

filename = "http://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data" 
housing <- read.table(filename) 
colnames(housing) <- c("CRIM", "ZN", "INDUS", "CHAS", "NOX",  
                       "RM", "AGE", "DIS", "RAD", "TAX", "PRATIO", 
                       "B", "LSTAT", "MDEV")

将数据分割：

housing <- housing[order(housing$MDEV),] 
#install.packages("caret") 
library(caret) 
set.seed(5557) 
indices <- createDataPartition(housing$MDEV, p=0.75, list=FALSE) 
training <- housing[indices,] 
testing <- housing[-indices,] 
nrow(training) 
nrow(testing)

计算我们的模型：

forestFit <- randomForest(MDEV ~ CRIM + ZN + INDUS + CHAS + NOX  
                  + RM + AGE + DIS + RAD + TAX + PRATIO  
                  + B + LSTAT, data=training) 
forestFit 
Call: 
 randomForest(formula = MDEV ~ CRIM + ZN + INDUS + CHAS + NOX +      RM + AGE + DIS + RAD + TAX + PRATIO + B + LSTAT, data = training)  
               Type of random forest: regression 
                     Number of trees: 500 
No. of variables tried at each split: 4 

          Mean of squared residuals: 11.16163 
                    % Var explained: 87.28

这是关于模型的更多信息性显示——我们看到模型解释了 87%的变量。

做出我们的预测：

forestPredict <- predict(forestFit, newdata=testing) 
See how well the model worked: 
diff <- forestPredict - testing$MDEV 
sum( (diff - mean(diff) )² ) #sum of squares 
1391.95553131418

这是我们在本章中产生的模型中，最小的平方和之一。

摘要

在本章中，我们使用了几种机器学习算法，其中一些在 R 和 Python 中进行对比。我们使用朴素贝叶斯算法来确定数据可能的使用方式。我们以不同的方式应用最近邻算法来查看我们的结果。我们使用决策树来构建预测算法。我们尝试用神经网络来解释房价。最后，我们使用随机森林算法做了相同的尝试——并取得了最佳结果！

在下一章中，我们将研究如何优化 Jupyter 笔记本。

第十章：优化 Jupyter Notebook

在开发 Jupyter Notebook 之前，你应该处理在公众开始访问之前需要进行的优化。这些优化包括从特定语言问题（使用最佳实践的 R 编码风格）到将笔记本部署到高可用环境中的一系列选项。

部署笔记本

Jupyter Notebook 是一个网站。你可以在用来显示此文档的计算机上托管网站。你所在部门可能有一台机器正在用作 Web 服务器。

如果你在本地机器上进行部署，你将拥有一个单用户网站，其他用户将被阻止访问，或者会发生冲突。发布笔记本的第一步是使用提供多用户访问的托管服务。

部署到 JupyterHub

当前主流的 Jupyter 托管产品是 JupyterHub。为了明确，JupyterHub 是安装在你控制的机器上。它为你的笔记本提供多用户访问。这意味着你可以在你的环境中的一台机器上安装 JupyterHub，并且只有内部用户（多个内部用户）可以访问它。

当 JupyterHub 启动时，它会开始一个 hub 或控制代理。该 hub 将启动一个监听器或代理实例来处理 Jupyter 请求。当代理收到 Jupyter 请求时，它将把请求交给 hub。如果 hub 判断这是一个新用户，它将生成一个新的 Jupyter 服务器实例，并将该用户与 Jupyter 之间的所有后续交互绑定到该用户的版本服务器上。此功能提供了多用户访问。

安装 JupyterHub

JupyterHub 需要 Python 3.3 或更高版本，我们将使用 Python 工具 pip3 来安装 JupyterHub。安装命令是：

npm install -g configurable-http-proxy
pip3 install jupyterhub

访问 JupyterHub 安装

现在我们可以直接从命令行启动 JupyterHub：

jupyterhub

替代你一直在使用的典型 jupyter notebook。

如前所述，JupyterHub 是用户访问笔记本的入口点。它在背后为用户实例化 Jupyter 实例。

用户界面也有一些小差异，如下图所示：

注意屏幕右上方的额外按钮：

控制面板：用于控制 JupyterHub
注销：注销将在 JupyterHub 的 SSL 功能中用于控制对笔记本的访问。

控制面板提供了两个选项：

停止我的服务器：停止 JupyterHub 的进一步响应。如果需要升级系统的某些部分，这个选项很有用。
我的服务器：返回到 JupyterHub 的主页（上一张截图）。

如果你在桌面上启动两个浏览器实例并访问 JupyterHub 上的笔记本，你将看到为每个用户呈现不同的信息，因为用户信息是分别分配给每个用户的。

Jupyter 托管

如果你想将笔记本分享给更广泛的观众，一个更好的机制是使用公共托管服务来展示你的笔记本。这可以减轻你自己做这件事时需要承担的一些大量安全性和可用性要求。

提供 Jupyter Notebook 托管服务的（托管）供应商越来越多。一些主要的供应商包括：

Rackspace：Rackspace 特别面向教育笔记本，提供为教育领域客户运行时的特殊处理
Azure：由微软提供的服务器托管服务—提供完整的笔记本托管
GitHub：主要的编程工件源仓库之一，例如你的笔记本

这些托管服务的一个伟大特点是，你不再需要安装从笔记本脚本中访问的任何库。托管公司已经下载（并维护和更新）了你通常会使用的所有库。GitHub 负责这些维护工作。其他托管公司则是纯粹的网络托管公司，期望你安装任何必要的软件。

优化你的脚本

你可以对笔记本脚本进行优化，使其运行更高效。这些优化是依赖于脚本语言的。我们已经在笔记本中讨论了 Python 和 R 脚本，并将讨论可以为这两种语言进行的优化。

Jupyter 确实支持其他语言，如 Scala 和 Spark。其他语言有自己的优化工具和策略。

优化你的 Python 脚本

性能调优你的 Python 脚本可以使用多种工具：

timeit
Python 正则表达式
字符串处理
循环优化
hotshot 性能分析

确定脚本执行所需的时间

Python 中的 timeit 函数接收一行代码并确定执行所需的时间。你还可以重复执行相同的脚本，查看是否有需要解决的启动问题。

timeit 以这种方式使用：

import timeit
t = timeit.Timer("myfunction('Hello World')", "import myfunction") 
t.timeit() 
3.32132323232
t.repeat(2, 2000000) 
[#, #, #]

repeat() 函数告诉 timeit 执行定时指令两次，每次执行 2,000,000 次。repeat() 调用后的数字是三次时间的展示，作为重复调用的代表。

你通常会用它来测试完整的函数，而不是与用户进行任何交互的部分。

使用 Python 正则表达式

有很多情况下你需要让脚本评估一个数据字符串，可能是由用户输入的。是的，你可以手动完成这个过程，但 Python 正则表达式处理非常高效，并能捕捉到你可能没有注意到的所有边界情况。

使用 Python 字符串处理

字符串在内存分配方面表现较差，每当对字符串进行微小的修改时，就会重复进行内存分配。此外，像将字母大写这样的小操作也会导致重新分配。最好一次性处理整个字符串（比如大写化），这样就能将重新分配的次数减少到一次。

最小化循环操作

在很多情况下，你可能会开发一个脚本，并将一些步骤放在循环内，因为在开发时这样做较为方便，而且对较小的测试数据运行时，性能差异可能不大。例如，在本书中的任何脚本里，我通常会提取一小部分数据行来处理。这样可以确保我正确地访问数据。然而，我可能会调整循环操作，使标记操作发生在其中。在大约 20 行的小数据集上，这种影响不大。然而，当我开始使用真正的数据集，可能包含数百万行时，每一行都需要设置标记，这将影响整个操作的性能。

然而，一旦处理较大的数据集时，那些每次循环执行时都会发生的小操作变得非常昂贵。通常有些操作可以从循环中提取出来，放到循环外执行一次。例如，如果我们想在循环中找到最大的数字，我们可以在循环外将结果初始化为一个不现实的值，在循环内部，如果看到这个不现实的值，就用第一个结果来初始化结果。这个测试会在每次循环操作时进行。然而，我们可以在进入循环之前就设置好第一个记录的结果。

分析脚本

热点分析器（hotshot profiler）在 Python 中可用，能够完整地展示脚本执行的过程。你需要先安装 Hotshot 才能执行。

它可以这样使用：

import hotshot
import myfunction
prof = hotshot.Profile('my_hotshot _stats')
prof.run('myfunction').close()

然后，要查看分析器的结果总结，可以使用以下命令：

import hotshot.stats
hotshot.stats.load('my_hotshot_stats').strip_dirs().sort_stats('time').print_stats()

优化你的 R 脚本

R 也提供了一些工具，帮助定位 R 代码的性能问题：

microbenchmark
修改常用函数
优化名称查找
优化数据框值提取
R 实现
更改算法

使用 microbenchmark 来分析 R 脚本

microbenchmark 是作为 R 库的一部分提供的。一旦将其包含在脚本中，你就可以将待分析的代码用 microbenchmark 标签包裹起来，执行后该工具会输出该脚本的性能分析信息。

例如，我们可以这样使用：

library(microbenchmark)
x <- runif(125) 
microbenchmark( mean(x) )

这样会对包含的代码进行 125 次（默认 100 次）执行，并输出诸如以下的分析信息：

Unit: nanoseconds
 expr  min     lq    mean median     uq   max neval
 sqrt(x)  825  860.5 1212.79  892.5  938.5 12905   100
 x⁰.5 3015 3059.5 3776.81 3101.5 3208.0 15215   100

我们关注的是均值，它能很好地指示每次迭代所需的时间。我们还应该注意到，若最小值和最大值相差很大，意味着存在明显的偏差——这正是我们在这里看到的情况。

修改提供的功能

R 允许你更改大多数对象的行为，包括提供的大多数知名函数。例如，一个极端的例子是重新设计mean()函数的工作方式。也许你对正在处理的数据有深入的理解，可以在特定情况下相应地提高性能。

与前面的工具一样，你可以先执行提供的功能，再进行你的实现，并比较 microbenchmark 提供的性能分析信息。

优化名称查找

R 允许动态命名对象，因此每次你访问一个变量时，R 必须遍历一系列作用域，找到对象当前所在的位置。

你也可以使用本地缓存机制直接访问确切的值，从而避免 R 在所有作用域中查找变量。

优化数据框值提取

R 提供多种编码风格，允许你以多种方式访问数据框中的特定值。如果你的代码频繁访问大型数据框以提取单个值，那么值得对不同的替代方法进行基准测试，看看哪种访问方式能提供最佳性能。

更改 R 实现

R 是由几家公司实现的，每个实现都有其自己的性能范围。在极端情况下，尝试不同的 R 实现可能会很有价值，以确定哪个最适合你。更改 Jupyter 使用的 R 引擎可能会比较困难。

更改算法

每当你编写笔记本代码时，你已采用某种方法来解决问题。这种方法被称为算法。你的算法可能包括遍历记录或直接查询数据库中的记录来获取感兴趣的记录。很多时候，你会在过程初期选择一个看似足够有效的算法。也许还有其他算法能够以更高效的方式解决问题。

对任何编程实现来说，最大的提升来自于改变整体处理方法。不幸的是，这是最难实现的变动，因为需要重新设计和重写代码。即便如此，像其他技巧一样，你仍然需要比较不同方法的基准，确保新方法更好。

监控 Jupyter

与本章之前关于优化的讨论一样，你也可以使用编程工具来监控笔记本的整体交互。Linux/Mac 环境下主要的工具是memory_profiler。如果你启动这个工具，然后再启动笔记本，分析器将会跟踪笔记本的内存使用情况。

根据此信息记录，如果您发现发生了大量内存使用，您可能能够调整程序的内存分配以减小内存占用。例如，性能分析器可能会突出显示您在一个循环中不断创建（和丢弃）一个大内存项。当您回到代码时，您会意识到这个内存访问可以从循环中提取出来，或者可以轻松地减少分配的大小。

缓存您的笔记本

缓存是一种常见的编程实践，用于加速性能。如果计算机不需要重新加载某一段代码、变量或文件，而是可以直接从缓存中访问，这将提高性能。

如果您将笔记本部署到 Docker 环境中，则有一种机制可以缓存您的笔记本。Docker 是一种在一台机器上虚拟化代码的机制，在 Java 编程领域，采用这种方法已经成为常规做法。幸运的是，Docker 非常灵活，已经确定了一种方法，可以将 Jupyter 用在 Docker 中。一旦进入 Docker，自动缓存页面就变得相对简单。底层工具是memcached，这也是一个常见的缓存工具，用于缓存任何内容，这里是 Jupyter 笔记本。

笔记本的安全

笔记本的安全性可以通过多种方法实现，例如：

管理授权
笔记本内容的安全

管理笔记本授权

笔记本可以通过用户名/密码授权进行保护。默认情况下，您的笔记本启用了授权。在 Jupyter 中，使用的是令牌/密码，而不是用户名/密码，因为令牌在解释时更灵活。请参阅 Jupyter 文档，了解如何实现授权，因为这一点随着时间的推移有所变化。

笔记本内容的安全

一个笔记本可能存在安全问题，Jupyter 会自动保护其中一些标准内容：

不信任的 HTML 会被清理
不信任的 JavaScript 不会被执行
Markdown 单元中的 HTML 和 JavaScript 不被信任
笔记本输出不被信任
笔记本中的其他 HTML 或 JavaScript 不被信任

信任问题归结为一个问题：是用户执行了这个操作，还是 Jupyter 脚本执行了？不信任意味着不会生成内容。

被清理过的代码会被包装，强制其值仅以文本显示—不会生成可执行代码。例如，如果您的笔记本单元格生成 HTML，如额外的H1标题标签，Jupyter 会清理输出，使其呈现原始 HTML，而不是页面上显示的 HTML 标题效果，原始 HTML 可能是类似<H1>Additional Heading</H1>的内容。

扩展 Jupyter 笔记本

扩展是指提供大量并发用户访问一个笔记本，而不降低性能。今天，唯一实现这一点的供应商是 Azure。他们每天都有成千上万的页面和用户在大规模使用。

最令人惊讶的是，这是一项免费的服务。

共享 Jupyter 笔记本

Jupyter 笔记本可以通过将笔记本放置在服务器上（有多种类型）或将笔记本转换为其他格式来共享（虽然它将不再具备交互性，但内容仍然可用）。

在笔记本服务器上共享 Jupyter 笔记本

笔记本配置内置了可以直接公开笔记本服务器的扩展功能。可以使用以下命令生成笔记本配置：

Jupyter Notebook -generate-config

在生成的 jupyter_notebook_config.py 文件中，有一些设置可以用于配置：

笔记本的 IP/端口地址
加密证书位置
密码

通过设置此项并启动 Jupyter，您应该能够在网络中的其他计算机上通过指定的 IP 地址访问笔记本。

在这样做之前，您应该与您的网络安全人员协作。

在笔记本服务器上共享加密的 Jupyter 笔记本

如果在前述配置文件中正确指定了证书信息，则笔记本只能通过 HTTPS 或安全的加密通道访问。

在 Web 服务器上共享笔记本

配置文件的另一部分是 tornado_settings。这一系列设置描述了将 web 流量引导到您的笔记本的 Web 服务器。再次强调，一旦这些设置完成，您可以通过 Web 服务器的 IP 地址访问您的笔记本。

这在将笔记本访问呈现为您现有网站的一部分时可能会很有用。

在 Docker 上共享笔记本

Docker 是一个框架，通过一个小的配置文件 Dockerfile 提供软件实例的托管服务。Docker 允许根据需要自动实例化多个软件实例。因此，我们将有多个笔记本实例可供用户使用。用户无法区分多个实例，因为他们只会引用同一个笔记本。Docker 根据用户最初连接到系统的情况，将流量引导到某个实例。

Dockerfile 包含环境设置，告诉 Docker 系统需要在实例中存在哪些系统组件，以便可以执行所引用的对象，在本例中是笔记本。

转换笔记本

您还可以通过将笔记本转换为可供接收者阅读的格式与他人共享笔记本。可以使用笔记本文件菜单中的“另存为”功能将笔记本转换为多种格式。

笔记本可以转换为以下格式：

<language> format：此选项取决于用于创建笔记本的语言。例如，R 笔记本可以选择将其下载为 R 脚本。
HTML：这种表示方式是 HTML 编码，用于显示笔记本中呈现页面的 HTML 构造。
Markdown：Markdown 是一种简单的显示标签格式，某些较旧的 Linux 系统使用该格式。
reST：另一种 Markdown 格式，它比 HTML 更简洁。
PDF。

版本控制笔记本

在编程领域，一个常见的做法是维护程序更改的历史记录。随着时间的推移，程序的不同版本会被保存在一个软件仓库中，程序员可以检索先前的版本，以便回到程序的旧有、可用状态。

在上一节中，我们提到过将笔记本放置在 GitHub 上。Git 是一个广泛使用的软件仓库，GitHub 是 Git 的一个基于互联网的实例。一旦你将任何软件放入 Git 中，它会自动进行版本管理。下次你更新 GitHub 中的笔记本时，Git 会将当前实例存储为历史版本，并将新实例放为当前版本——这样，任何访问你 GitHub 仓库的人都会默认看到最新版本。

总结

在本章中，我们将笔记本部署到了一组不同的环境中。我们探讨了可以对笔记本脚本进行的优化。我们学习了不同的分享笔记本的方式。最后，我们研究了如何将笔记本转换为没有 Jupyter 访问权限的用户使用的格式。

posted @ 2025-10-08 11:37 绝不原创的飞龙阅读(9) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈