IPython-交互式计算和可视化秘籍第二版-三-

IPython 交互式计算和可视化秘籍第二版（三）

原文：annas-archive.org/md5/62516af4e05f6a3b0523d8aa07fd5acb

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章：高级可视化

本章将涵盖以下主题：

• 使用 prettyplotlib 创建更漂亮的 matplotlib 图形

• 使用 seaborn 创建漂亮的统计图

• 使用 Bokeh 创建交互式 Web 可视化

• 在 IPython 笔记本中使用 D3.js 可视化 NetworkX 图形

• 使用 mpld3 将 matplotlib 图形转换为 D3.js 可视化

• 使用 Vispy 开始进行高性能交互式数据可视化

介绍

可视化是本书的一个核心主题。我们在大多数食谱中创建图形，因为这是传达定量信息的最有效方式。在大多数情况下，我们使用 matplotlib 来创建图表。在本章中，我们将看到 Python 中更高级的可视化功能。

首先，我们将看到几个包，它们可以让我们改善 matplotlib 图形的默认样式和类似 MATLAB 的 pyplot 接口。还有其他一些高层次的可视化编程接口，在某些情况下可能更方便。

此外，Web 平台正在越来越接近 Python。IPython 笔记本就是这种趋势的一个很好例子。在本章中，我们将看到一些技术和库，帮助我们在 Python 中创建交互式 Web 可视化。这些技术让我们能够将 Python 在数据分析方面的强大功能与 Web 在交互性方面的优势结合起来。

最后，我们将介绍 Vispy，这是一个新的高性能交互式可视化库，专为大数据而设计。

使用 prettyplotlib 创建更漂亮的 matplotlib 图形

matplotlib 有时因其图形的默认外观而受到批评。例如，默认的颜色映射既不具美学吸引力，也没有清晰的感知信息。

有很多方法试图规避这个问题。在本食谱中，我们将介绍prettyplotlib，这是由 Olga Botvinnik 创建的。这款轻量级 Python 库显著改善了许多 matplotlib 图形的默认样式。

准备工作

你可以在项目页面找到 prettyplotlib 的安装说明，网址是 github.com/olgabot/prettyplotlib。你基本上只需在终端中执行 pip install prettyplotlib 即可。

如何实现……

1. 首先，让我们导入 NumPy 和 matplotlib：

   In [1]: import numpy as np
           import matplotlib.pyplot as plt
           import matplotlib as mpl
           %matplotlib inline
   

2. 然后，我们用 matplotlib 绘制几条曲线：

   In [2]: np.random.seed(12)
           for i in range(8):
               x = np.arange(1000)
               y = np.random.randn(1000).cumsum()
               plt.plot(x, y, label=str(i))
           plt.legend()
   

   

   注意

   如果你正在阅读这本书的印刷版，你将看不到颜色。你可以在本书的网站上找到彩色图像。

3. 现在，我们用 prettyplotlib 创建完全相同的图表。我们只需将 matplotlib.pyplot 命名空间替换为 prettyplotlib：

   In [3]: import prettyplotlib as ppl
           np.random.seed(12)
           for i in range(8):
               x = np.arange(1000)
               y = np.random.randn(1000).cumsum()
               ppl.plot(x, y, label=str(i))
           ppl.legend()
   

   

4. 让我们通过一个图像展示另一个示例。我们首先使用 matplotlib 的 pcolormesh() 函数将 2D 数组显示为图像：

   In [4]: np.random.seed(12)
           plt.pcolormesh(np.random.rand(16, 16))
           plt.colorbar()
   

   

   默认的彩虹颜色映射被认为会导致可视化数据被误解。

5. 现在，我们使用 prettyplotlib 显示完全相同的图像：

   In [5]: np.random.seed(12)
           ppl.pcolormesh(np.random.rand(16, 16))
   

   

   这种可视化方式要更加清晰，因为高值或低值比彩虹色图更明显。

它是如何工作的…

prettyplotlib 只是对 matplotlib 的默认样式选项进行了微调。它的绘图接口基本上与 matplotlib 相同。要了解如何修改 matplotlib 的样式，值得查看 prettyplotlib 的代码。

还有更多内容…

改进 matplotlib 样式的其他方法有很多：

• Randal Olson 的一篇博客文章解释了如何使用 matplotlib 创建干净且美观的图表；你可以在 www.randalolson.com/2014/06/28/how-to-make-beautiful-data-visualizations-in-python-with-matplotlib/ 阅读到这篇文章。

• matplotlib 正在进行一些工作，添加样式表支持；更多信息可以在 github.com/matplotlib/matplotlib/blob/master/doc/users/style_sheets.rst 找到。

• 关于为什么彩虹色图是误导性的，可以查看 eagereyes.org/basics/rainbow-color-map

另见

• 使用 seaborn 创建美丽的统计图 方案

使用 seaborn 创建美丽的统计图

matplotlib 配有一个名为 pyplot 的高级绘图 API。受 MATLAB 启发（MATLAB 是一个广泛使用的数值计算商业软件），这个接口对于科学家来说可能有些过于底层，因为它可能会导致难以阅读和维护的样板代码。然而，它可能是科学 Python 社区中最广泛使用的绘图接口之一。

存在更高级、更便捷的绘图接口。在这个方案中，我们介绍了 seaborn，它由 Michael Waskom 创建。这个库提供了一个高层次的绘图 API，专门为统计图形量身定制，同时与 pandas 紧密集成。

准备工作

你可以在项目页面 github.com/mwaskom/seaborn 上找到 seaborn 的安装说明。你只需要在终端输入 pip install seaborn。

如何操作…

1. 让我们导入 NumPy、matplotlib 和 seaborn：

   In [1]: import numpy as np
           import matplotlib.pyplot as plt
           import seaborn as sns
           %matplotlib inline
   

2. 我们生成一个随机数据集（参考 seaborn 网站上的示例 nbviewer.ipython.org/github/mwaskom/seaborn/blob/master/examples/linear_models.ipynb）：

   In [2]: x1 = np.random.randn(80)
           x2 = np.random.randn(80)
           x3 = x1 * x2
           y1 = .5 + 2 * x1 - x2 + 2.5 * x3 + \
                3 * np.random.randn(80)
           y2 = .5 + 2 * x1 - x2 + 2.5 * np.random.randn(80)
           y3 = y2 + np.random.randn(80)
   

3. Seaborn 实现了许多易于使用的统计绘图函数。例如，下面是如何创建一个小提琴图。这种类型的图可以展示数据点的详细分布，而不仅仅是像箱形图那样显示四分位数：

   In [3]: sns.violinplot([x1,x2, x3])
   

   

4. Seaborn 也实现了全功能的统计可视化函数。例如，我们可以使用一个单独的函数（regplot()）来执行并且显示两个变量之间的线性回归：

   In [4]: sns.regplot(x2, y2)
   

   

5. Seaborn 内建对 pandas 数据结构的支持。在这里，我们显示了 DataFrame 对象中所有变量之间的成对相关性：

   In [5]: df = pd.DataFrame(dict(x1=x1, x2=x2, x3=x3, 
                                  y1=y1, y2=y2, y3=y3))
           sns.corrplot(df)
   

   

还有更多…

除了 seaborn，还有其他高级绘图接口：

• The Grammar of Graphics 是 Dr. Leland Wilkinson 撰写的一本书，影响了许多高级绘图接口，如 R 的 ggplot2、Python 中 yhat 的 ggplot 等。

• Vega，由 Trifacta 提供，是一种声明式可视化语法，可以转换为 D3.js（一种 JavaScript 可视化库）。此外，Vincent 是一个 Python 库，让我们使用 Vega 创建可视化。

• Tableau 的 VizQL 是一种面向商业数据库的可视化语言。

这里还有更多参考资料：

• Vega 可在 trifacta.github.io/vega/ 获取

• Vincent 可在 vincent.readthedocs.org/en/latest/ 获取

• ggplot2 可在 ggplot2.org/ 获取

• Python 的 ggplot 可在 blog.yhathq.com/posts/ggplot-for-python.html 获取

• VizQL 可在 www.tableausoftware.com/fr-fr/products/technology 获取

参见

• 使用 prettyplotlib 美化 matplotlib 图形 的配方

使用 Bokeh 创建交互式网页可视化

Bokeh 是一个用于在浏览器中创建丰富交互式可视化的库。图形在 Python 中设计，并完全在浏览器中渲染。在本配方中，我们将学习如何在 IPython notebook 中创建并渲染交互式 Bokeh 图形。

准备工作

按照网站上的说明在 bokeh.pydata.org 安装 Bokeh。原则上，你可以在终端中输入 pip install bokeh。在 Windows 上，你也可以从 Chris Gohlke 的网站下载二进制安装程序，网址为 www.lfd.uci.edu/~gohlke/pythonlibs/#bokeh。

如何实现…

1. 让我们导入 NumPy 和 Bokeh。我们需要调用 output_notebook() 函数来告诉 Bokeh 在 IPython notebook 中渲染图形：

   In [1]: import numpy as np
           import bokeh.plotting as bkh
           bkh.output_notebook()
   

2. 我们创建一些随机数据：

   In [2]: x = np.linspace(0., 1., 100)
           y = np.cumsum(np.random.randn(100))
   

3. 让我们画一条曲线：

   In [3]: bkh.line(x, y, line_width=5)
           bkh.show()
   

   一个交互式图形已在 notebook 中渲染。我们可以通过点击图形上方的按钮来平移和缩放：

   

   使用 Bokeh 创建的交互式图形

4. 让我们继续另一个例子。我们首先加载一个示例数据集（鸢尾花）。我们还根据花的种类生成一些颜色：

   In [4]: from bokeh.sampledata.iris import flowers
           colormap = {'setosa': 'red',
                       'versicolor': 'green',
                       'virginica': 'blue'}
           flowers['color'] = flowers['species'].map(
                                      lambda x: colormap[x])
   

5. 现在，我们渲染一个交互式散点图：

   In [5]: bkh.scatter(flowers["petal_length"], 
                       flowers["petal_width"],
                       color=flowers["color"], 
                       fill_alpha=0.25, size=10,)
           bkh.show()
   

   

   使用 Bokeh 的交互式散点图

还有更多…

即使没有 Python 服务器，Bokeh 图表在笔记本中也是互动的。例如，我们的图表可以在 nbviewer 中进行互动。Bokeh 还可以从我们的图表生成独立的 HTML/JavaScript 文档。更多示例可以在图库中找到，网址为 bokeh.pydata.org/docs/gallery.html。

Bokeh 提供了一个 IPython 扩展，简化了在笔记本中集成互动图表的过程。这个扩展可以在 github.com/ContinuumIO/bokeh/tree/master/extensions 找到。

同样地，我们还要提到 plot.ly，这是一个在线商业服务，提供用于 Web 基于的互动可视化的 Python 接口，网址为 plot.ly。

另见

• 使用 mpld3 将 matplotlib 图表转换为 D3.js 可视化 这个食谱

使用 D3.js 在 IPython 笔记本中可视化 NetworkX 图表

D3.js (d3js.org) 是一个流行的 Web 互动可视化框架。它是用 JavaScript 编写的，允许我们基于 Web 技术（如 HTML、SVG 和 CSS）创建数据驱动的可视化。虽然还有许多其他 JavaScript 可视化和图表库，但在本食谱中我们将重点介绍 D3.js。

作为纯 JavaScript 库，D3.js 原则上与 Python 无关。然而，基于 HTML 的 IPython 笔记本可以无缝集成 D3.js 可视化。

在这个食谱中，我们将使用 Python 的 NetworkX 创建一个图，并在 IPython 笔记本中使用 D3.js 进行可视化。

准备工作

本食谱需要你了解 HTML、JavaScript 和 D3.js 的基础知识。

如何操作…

1. 让我们导入所需的包：

   In [1]: import json
           import numpy as np
           import networkx as nx
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们加载了一个著名的社交图，该图于 1977 年发布，名为 Zachary's Karate Club graph。这个图展示了空手道俱乐部成员之间的友谊。俱乐部的主席和教练发生了争执，导致该小组分裂。在这里，我们只用 matplotlib（使用 networkx.draw() 函数）显示图表：

   In [2]: g = nx.karate_club_graph()
           nx.draw(g)
   

   

3. 现在，我们将使用 D3.js 在笔记本中显示这个图。第一步是将图表导入到 JavaScript 中。我们选择将图表导出为 JSON 格式。D3.js 通常期望每个边都是一个具有源（source）和目标（target）的对象。此外，我们还指定了每个成员所代表的“俱乐部”属性（club）。NetworkX 提供了一个内置的导出功能，我们可以在这里使用：

   In [3]: from networkx.readwrite import json_graph
           data = json_graph.node_link_data(g)
           with open('graph.json', 'w') as f:
               json.dump(data, f, indent=4)
   

4. 下一步是创建一个 HTML 对象，容纳可视化图表。我们在笔记本中创建了一个 <div> 元素，并为节点和链接（也叫边）指定了一些 CSS 样式：

   In [4]: %%html
           <div id="d3-example"></div>
           <style>
           .node {stroke: #fff; stroke-width: 1.5px;}
           .link {stroke: #999; stroke-opacity: .6;}
           </style>
   

5. 最后一步比较复杂。我们编写 JavaScript 代码，从 JSON 文件加载图表并使用 D3.js 显示它。这里需要了解 D3.js 的基础知识（请参阅 D3.js 的文档）。代码较长，您可以在本书网站上找到完整代码。在这里，我们突出了最重要的步骤：

   In [5]: %%javascript
       // We load the d3.js library. 
       require(["d3"], function(d3) {
           // The code in this block is executed when the 
           // d3.js library has been loaded.
           [...]
           // We create a force-directed dynamic graph 
           // layout.
           var force = d3.layout.force().charge(-120).
                  linkDistance(30).size([width, height]);
           [...]
           // In the <div> element, we create a <svg> graphic
           // that will contain our interactive 
           // visualization.
           var svg = d3.select("#d3-example").select("svg");
           [...]
           // We load the JSON file.
           d3.json("graph.json", function(error, graph) {
               // We create the graph here.
               force.nodes(graph.nodes).links(graph.links)
                    .start();
   
               // We create a <line> SVG element for each
               // link in the graph.
               var link = svg.selectAll(".link")
                       .data(graph.links)
                       .enter().append("line")
                       .attr("class", "link");
   
               // We create <circle> SVG elements for the 
               // nodes.
               var node = svg.selectAll(".node")
                   .data(graph.nodes)
                   .enter().append("circle")
                   [...]
                   .style("fill", function(d) {
                       return color(d.club); 
                    })
                   .call(force.drag);
                   [...]
               });
           });
   

   当我们执行此单元格时，前一个单元格中创建的 HTML 对象会被更新。图表是动画的并且是互动的；我们可以点击节点，查看其标签，并在画布中移动它们：

   

   使用 D3.js 在笔记本中的互动图

还有更多…

D3.js 的画廊包含更多 Web 上美丽的互动可视化示例。它们可以在 github.com/mbostock/d3/wiki/Gallery 找到。

在这个示例中，我们通过一个静态数据集创建了一个 HTML/JavaScript 互动可视化。使用 IPython 2.0 及以上版本，我们还可以创建涉及浏览器与 Python 内核之间双向通信的动态、实时可视化。Brian Granger 提供了一个实验性实现，可以在 nbviewer.ipython.org/github/ellisonbg/talk-2014-strata-sc/blob/master/Graph%20Widget.ipynb 访问。

另外，我们还要提到 Vincent，一个 Python 到 Vega 的翻译器。Vega 是一种基于 JSON 的可视化语法，可以被翻译为 D3.js。Vincent 使得在 Python 中设计互动可视化并在浏览器中渲染成为可能。更多信息可以在 vincent.readthedocs.org/en/latest/ 找到。

另见

• 使用 Bokeh 创建互动 Web 可视化 示例

• 将 matplotlib 图形转换为 D3.js 可视化图表，使用 mpld3 示例

使用 mpld3 将 matplotlib 图形转换为 D3.js 可视化

mpld3 库会自动将 matplotlib 图形转换为互动 D3.js 可视化图表。在本示例中，我们将展示如何在笔记本中使用该库。

准备工作

要安装 mpld3 库，只需在终端中输入 pip install mpld3。更多信息请参见主网站 mpld3.github.io。

如何做…

1. 首先，我们像往常一样加载 NumPy 和 matplotlib：

   In [1]: import numpy as np
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 然后，我们通过一个函数调用在笔记本中启用 mpld3 图形：

   In [2]: from mpld3 import enable_notebook
           enable_notebook()
   

3. 现在，让我们使用 matplotlib 创建一个散点图：

   In [3]: X = np.random.normal(0, 1, (100, 3))
           color = np.random.random(100)
           size = 500 * np.random.random(100)
           plt.scatter(X[:,0], X[:,1], c=color,
                       s=size, alpha=0.5, linewidths=2)
           plt.grid(color='lightgray', alpha=0.7)
   

   matplotlib 图形通过 D3.js 渲染，而不是使用标准的 matplotlib 后端。特别是，图形是互动的（我们可以平移和缩放图形）：

   

   使用 mpld3 的互动 matplotlib 图

4. 现在，我们创建一个更复杂的示例，其中包含多个子图，表示一个 3D 数据集的不同 2D 投影。我们使用 matplotlib 的 subplots() 函数中的 sharex 和 sharey 关键字来自动绑定不同图形的 x 和 y 轴。在任何子图上进行平移和缩放都会自动更新所有其他子图：

   In [4]: fig, ax = plt.subplots(3, 3, figsize=(6, 6),
                                  sharex=True, sharey=True)
           fig.subplots_adjust(hspace=0.3)
           X[::2,2] += 3
           for i in range(3):
               for j in range(3):
                   ax[i,j].scatter(X[:,i], X[:,j], c=color,
                       s=.1*size, alpha=0.5, linewidths=2)
                   ax[i,j].grid(color='lightgray', alpha=0.7)
   

   这个用例完全可以通过 mpld3 处理；D3.js 子图是动态联动的：

   

   在 mpld3 中的互动联动子图

它是如何工作的…

mpld3 的工作原理是首先爬取并将 matplotlib 图形导出为 JSON（在 mplexporter 框架的上下文中）。然后，库从该 JSON 表示生成 D3.js 代码。这种架构可以支持除了 D3.js 之外的其他 matplotlib 后端。

还有更多…

这里有一些参考资料：

• mplexporter 可用：github.com/mpld3/mplexporter

• mpld3 在 GitHub 上可用：github.com/jakevdp/mpld3

另见

• 使用 Bokeh 创建交互式 Web 可视化 食谱

• 在 IPython 笔记本中使用 D3.js 可视化 NetworkX 图 食谱

入门 Vispy 进行高性能互动数据可视化

大多数现有的 Python 绘图或可视化库可以显示小型或中型数据集（包含不超过几万个点的数据集）。在 大数据 时代，有时需要显示更大的数据集。

Vispy (vispy.org) 是一个年轻的 2D/3D 高性能可视化库，可以显示非常大的数据集。Vispy 通过 OpenGL 库利用现代图形处理单元（GPU）的计算能力。

在过去的二十年里，视频游戏行业促进了 GPU 的强大功能。GPU 专注于高性能、实时渲染。因此，它们非常适合互动式科学绘图。

Vispy 提供了一个 Pythonic 的面向对象接口，用于 OpenGL，适用于那些了解 OpenGL 或愿意学习 OpenGL 的人。更高级别的图形接口也在写作时开发中，实验版本已经可以使用。这些接口不需要任何 OpenGL 知识。

在本食谱中，我们将简要介绍 OpenGL 的基本概念。在以下两种情况中，你需要了解这些概念：

• 如果你今天想使用 Vispy，在高级绘图接口可用之前

• 如果你想创建自定义的、复杂的、高性能的可视化，这些可视化在 Vispy 中尚未实现

在这里，我们使用 Vispy 的面向对象接口来显示一个数字信号。

准备就绪

Vispy 依赖于 NumPy。需要一个后端库（例如，PyQt4 或 PySide）。

本食谱已在 github.com/vispy/vispy 上可用的 Vispy 开发版本中进行测试。你应该克隆 GitHub 仓库并使用以下命令安装 Vispy：

python setup.py install

本食谱中使用的 API 可能会在未来版本中发生变化。

如何做到…

1. 让我们导入 NumPy，vispy.app（用于显示画布）和 vispy.gloo（面向对象的 OpenGL 接口）：

   In [1]: import numpy as np
           from vispy import app
           from vispy import gloo
   

2. 为了显示一个窗口，我们需要创建一个 画布：

   In [2]: c = app.Canvas(keys='interactive')
   

3. 使用vispy.gloo时，我们需要编写着色器。这些程序使用类似 C 语言的语言编写，运行在 GPU 上，为我们的可视化提供完全的灵活性。在这里，我们创建一个简单的顶点着色器，它直接在画布上显示 2D 数据点（存储在a_position变量中）。接下来我们将在下一节中看到更多细节：

   In [3]: vertex = """
           attribute vec2 a_position;
           void main (void)
           {
               gl_Position = vec4(a_position, 0.0, 1.0);
           }
           """
   

4. 我们需要创建的另一个着色器是片段着色器。它让我们控制像素的颜色。在这里，我们将所有数据点显示为黑色：

   In [4]: fragment = """
           void main()
           {
               gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
           }
           """
   

5. 接下来，我们创建一个OpenGL Program。这个对象包含着色器并将着色器变量链接到 NumPy 数据：

   In [5]: program = gloo.Program(vertex, fragment)
   

6. 我们将a_position变量链接到一个*(1000, 2)的 NumPy 数组，该数组包含 1000 个数据点的坐标。在默认坐标系中，四个画布角的坐标为(+/-1, +/-1)*：

   In [6]: program['a_position'] = np.c_[
                   np.linspace(-1.0, +1.0, 1000),
                   np.random.uniform(-0.5, +0.5, 1000)]
   

7. 我们在窗口调整大小时创建一个回调函数。更新OpenGL 视口可以确保 Vispy 使用整个画布：

   In [7]: @c.connect
           def on_resize(event):
               gloo.set_viewport(0, 0, *event.size)
   

8. 当画布需要刷新时，我们创建一个回调函数。这个on_draw()函数渲染整个场景。首先，我们将窗口清空为白色（每一帧都需要这样做）。然后，我们使用 OpenGL 程序绘制一系列线段：

   In [8]: @c.connect
           def on_draw(event):
               gloo.clear((1,1,1,1))
               program.draw('line_strip')
   

9. 最后，我们显示画布并运行应用程序：

   In [9]: c.show()
           app.run()
   

   下图显示了一个截图：

   

   使用 Vispy 的基本可视化示例

它是如何工作的……

OpenGL 是一种硬件加速的互动可视化开放标准。它广泛应用于视频游戏、工业（计算机辅助设计，或CAD）、虚拟现实和科学应用（医学成像、计算机图形学等）。

OpenGL 是一项成熟的技术，创建于 1990 年代初期。在 2000 年代初，OpenGL 2.0 引入了一个重大的新特性：可以自定义渲染管线的基本步骤。这个管线定义了数据如何在 GPU 上处理，以进行实时渲染。许多 OpenGL 课程和教程讲解的是旧的、固定的管线。然而，Vispy 仅支持现代的可编程管线。

在这里，我们将介绍本食谱中使用的可编程管线的基本概念。OpenGL 比我们在这里能覆盖的要复杂得多。然而，Vispy 为 OpenGL 的最常见功能提供了一个大大简化的 API。

注意

Vispy 基于OpenGL ES 2.0，这是 OpenGL 的一种变体，支持桌面计算机、移动设备和现代网页浏览器（通过WebGL）。现代图形卡可以支持额外的功能。这些功能将在未来版本的 Vispy 中提供。

在给定 OpenGL 程序的渲染管线中，有四个主要元素：

• 数据缓冲区将数值数据存储在 GPU 上。缓冲区的主要类型有顶点缓冲区、索引缓冲区和纹理。

• 变量在着色器中是可用的。主要有四种类型的变量：属性、常量、变化量和纹理采样器。

• 着色器是用一种类似 C 语言的语言编写的 GPU 程序，称为OpenGL 着色语言（GLSL）。着色器的两种主要类型是顶点着色器和片段着色器。

• 图元类型定义了数据点的渲染方式。主要类型有点、线和三角形。

这是渲染管线的工作方式：

1. 数据被发送到 GPU 并存储在缓冲区中。

2. 顶点着色器并行处理数据，并生成多个 4D 点，这些点位于归一化坐标系中（+/-1, +/-1）。第四维是齐次坐标（通常为 1）。

3. 图形图元（点、线和三角形）是通过顶点着色器返回的数据点生成的（图元组装和光栅化）。

4. 片段着色器并行处理所有图元像素，并返回每个像素的颜色作为 RGBA 组件。

在本例中，只有一个 GPU 变量：a_position属性。属性是每个数据点取一个值的变量。常量是全局变量（所有数据点共享），而变化量则用于将值从顶点着色器传递到片段着色器（通过对两个或三个顶点之间的像素进行自动线性插值）。

在vispy.gloo中，Program是通过顶点着色器和片段着色器创建的。然后，可以使用program['varname'] = value语法设置在着色器中声明的变量。当varname是一个属性变量时，值可以是一个 NumPy 二维数组。在这个数组中，每一行包含每个数据点的组成部分。

同样，我们也可以在程序中声明常量和纹理。

最后，program.draw()函数使用指定的图元类型渲染数据。这里，line_strip图元类型告诉 GPU 遍历所有顶点（由顶点缓冲区返回），并绘制从一个点到下一个点的线段。如果有n个点，则会有n-1条线段。

其他图元类型包括点和三角形，有几种方法可以从顶点列表中生成线或三角形。

此外，还可以提供索引缓冲区。索引缓冲区包含指向顶点缓冲区的索引。使用索引缓冲区可以让我们在图元组装阶段多次复用同一个顶点。例如，当使用triangles图元类型渲染一个立方体时（每三个点生成一个三角形），我们可以使用包含八个数据点的顶点缓冲区和包含三十六个索引的索引缓冲区（三个点构成一个三角形，每个面有两个三角形，共六个面）。

还有更多内容……

这里展示的示例非常简单。然而，OpenGL 和 Vispy 提供的方法却特别强大。它使我们能够完全控制渲染管线，并且几乎可以以最优方式利用 GPU 的计算能力。

高性能是通过最小化数据传输到 GPU 的次数来实现的。当显示静态数据（例如散点图）时，可以仅在初始化时将数据发送到 GPU。而渲染动态数据的速度也相当快；数据传输的数量级大约为 1 GBps。

此外，尽量减少 OpenGL 绘制调用的次数至关重要。每次绘制都会产生显著的开销。通过一次性渲染所有相似的原始类型来实现高性能（批量渲染）。即使点的属性不同（例如，不同大小和颜色的点），GPU 在批量渲染时也特别高效。

最后，可以通过着色器在 GPU 上以非常高的性能执行几何或像素变换。当变换在着色器中实现时，GPU 强大的架构（由数百或数千个计算单元组成）得到了充分利用。

在可视化的上下文中，可以在着色器中执行通用计算。与适当的 GPGPU 框架（如 CUDA 或 OpenCL）相比，存在一个主要缺点：在顶点着色器中，给定线程只能访问一个数据点。同样，在片段着色器中，线程只能访问一个像素。然而，某些类型的仿真或可视化效果需要顶点或像素之间的交互。虽然有方法可以缓解这个问题，但这会导致性能下降。

然而，OpenGL 可以与 CUDA/OpenCL 进行互操作。缓冲区可以在 OpenGL 和 GPGPU 框架之间共享。复杂的 CUDA/OpenCL 计算可以实时地在顶点缓冲区或纹理上实现，从而实现高效的数值仿真渲染。

Vispy 用于科学可视化

正如我们在这个示例中看到的，Vispy 要求用户了解 OpenGL 和 GLSL。然而，当前正在开发更高级的图形接口。这些接口将为科学家们带来 GPU 的强大能力，用于高性能交互式可视化。

视觉组件将提供可重用的、响应式的图形组件，如形状、多边形、3D 网格、图形等。这些组件将是完全可定制的，并且可以在不需要了解 OpenGL 的情况下使用。着色器组合系统将允许高级用户以模块化的方式重用 GLSL 代码片段。

视觉组件将被组织在一个 场景图 中，执行基于 GPU 的 变换。

科学绘图接口将会实现。Vispy 还可以作为现有绘图库（如 matplotlib）的高性能后端。

Vispy 还将支持通过 WebGL 在 IPython notebook 中的完全集成。

最终，Vispy 将能够实现多种科学可视化：

• 散点图可以通过点精灵高效渲染，每个数据点使用一个顶点。平移和缩放可以在顶点着色器中实现，从而实现对数百万点的快速交互式可视化。

• 静态或动态（实时）数字信号可以通过折线显示。使用 Anti-Grain Geometry 的 OpenGL 实现可以实现曲线的高质量渲染，这是一个高质量的 2D 渲染库。

• 图表可以通过组合点和线段来显示。

• 3D 网格可以使用三角形和索引缓冲区显示。几何变换和逼真的光照可以在顶点和片段着色器中实现。

• 实时图像流可以有效地通过纹理显示。

• 轴线、网格、刻度、文本和标签可以在片段着色器中高效渲染。

在 Vispy 的画廊中可以找到许多示例。

以下是一些参考资料：

• Vispy 的画廊位于vispy.org/gallery.html

• 由 Nicolas P. Rougier 编写的现代 OpenGL 教程，位于www.loria.fr/~rougier/teaching/opengl/

• Python 中的神经科学硬件加速交互式数据可视化，一篇文章位于journal.frontiersin.org/Journal/10.3389/fninf.2013.00036/full

• Vispy 用户邮件列表位于groups.google.com/forum/#!forum/vispy

• Vispy-dev 邮件列表位于groups.google.com/forum/#!forum/vispy-dev

• Anti-Grain Geometry 库在维基百科上的页面，位于en.wikipedia.org/wiki/Anti-Grain_Geometry

第七章。统计数据分析

本章将涵盖以下主题：

• 使用 pandas 和 matplotlib 探索数据集

• 开始进行统计假设检验——简单的 z 检验

• 开始使用贝叶斯方法

• 使用列联表和卡方检验估计两个变量之间的相关性

• 使用最大似然法拟合数据的概率分布

• 使用核密度估计法非参数地估计概率分布

• 通过马尔可夫链蒙特卡洛方法从后验分布中抽样拟合贝叶斯模型

• 在 IPython 笔记本中使用 R 编程语言分析数据

介绍

在前面的章节中，我们回顾了 Python 中高性能交互式计算的技术方面。现在，我们开始本书的第二部分，展示可以使用 Python 解决的各种科学问题。

在本章中，我们介绍了用于数据分析的统计方法。除了涉及诸如 pandas、statsmodels 和 PyMC 等统计包外，我们还将解释这些方法背后的数学原理基础。因此，如果你有一定的概率论和微积分基础，本章将会最有益。

下一章，第八章，机器学习，与本章密切相关；其基础数学非常相似，但目标略有不同。在本章中，我们展示了如何从现实世界数据中获得洞见，以及如何在不确定性面前做出明智的决策。而在下一章，目标是从数据中学习，即从部分观察中进行归纳和预测结果。

在本引言中，我们将对本章中所涉及的方法进行广泛的、高层次的概述。

什么是统计数据分析？

统计数据分析的目标是从部分和不确定的观察中理解复杂的现实现象。数据中的不确定性导致我们对现象的知识也存在不确定性。理论的主要目标是量化这种不确定性。

在进行统计数据分析时，重要的是区分数学理论与分析后做出的决策。前者是完美严谨的；或许令人惊讶的是，数学家们能够建立一个精确的数学框架来处理不确定性。然而，统计分析得出的实际人类决策中有主观的部分。在决策过程中，理解统计结果背后的风险和不确定性至关重要。

在本章中，我们将看到统计数据分析背后的基本概念、原则和理论，特别是如何在量化风险的情况下做出决策。当然，我们将始终展示如何使用 Python 实现这些方法。

一些术语

在我们开始这些食谱之前，有许多术语需要介绍。这些概念使我们能够在多个维度上对统计技术进行分类。

探索、推断、决策和预测

探索性方法 使我们能够通过基本的统计聚合和交互式可视化对数据集进行初步查看。我们在本书的第一章以及《学习 IPython 用于交互式计算和数据可视化》一书（Packt Publishing）中介绍了这些基本方法。本章的第一条食谱，使用 pandas 和 matplotlib 探索数据集，展示了另一个例子。

统计推断 是通过部分和不确定的观察获取关于未知过程的信息。特别地，估计 是指为描述该过程的数学变量获得近似值。本章的三条食谱涉及统计推断：

• 通过最大似然法拟合概率分布 食谱

• 通过核密度估计非参数地估计概率分布 食谱

• 通过从后验分布中采样来拟合贝叶斯模型，使用马尔可夫链蒙特卡罗方法 食谱

决策理论 使我们能够通过随机观察在可控风险下对未知过程做出决策。以下两条食谱展示了如何做出统计决策：

• 开始进行统计假设检验：简单的 z 检验 食谱

• 通过列联表和卡方检验估计两个变量之间的相关性 食谱

预测 是指从数据中学习，即根据有限的观察预测随机过程的结果。这是下章的主题，第八章，机器学习。

单变量和多变量方法

在大多数情况下，你可以考虑数据中的两个维度：

• 观察值（或 样本，对于机器学习人员）

• 变量（或 特征）

通常，观察值是同一随机过程的独立实现。每个观察值由一个或多个变量组成。大多数时候，变量要么是数字，要么是属于有限集合的元素（即取有限数量的值）。分析的第一步是理解你的观察值和变量是什么。

如果你有一个变量，则你的问题是 单变量。如果你有两个变量，则是 双变量，如果你至少有两个变量，则是 多变量。单变量方法通常较为简单。话虽如此，单变量方法也可以应用于多变量数据，每次使用一个维度。尽管这种方法无法探索变量之间的相互作用，但它通常是一个有趣的初步方法。

频率学派和贝叶斯方法

至少有两种不同的方式来考虑不确定性，从而产生两类用于推断、决策和其他统计问题的方法。它们分别被称为频率派方法和贝叶斯方法。一些人偏好频率派方法，而另一些人则偏好贝叶斯方法。

频率派解释概率为多个独立试验的统计平均（大数法则）。贝叶斯派则解释为信念程度（不需要多个试验）。贝叶斯解释在只考虑单次试验时非常有用。此外，贝叶斯理论还考虑了我们对随机过程的先验知识。随着数据的增加，先验概率分布会被更新为后验分布。

频率派和贝叶斯派方法各有其优缺点。例如，有人可能会说，频率派方法比贝叶斯方法更容易应用，但解释起来更为困难。有关频率派方法的经典误用，见www.refsmmat.com/statistics/。

无论如何，如果你是统计数据分析的初学者，你可能希望在选择立场之前学习这两种方法的基础知识。本章将向你介绍这两种方法。

以下的配方完全是贝叶斯方法：

• 入门贝叶斯方法配方

• 通过马尔可夫链蒙特卡洛方法从后验分布中采样拟合贝叶斯模型配方

Jake Vanderplas 曾写过几篇关于频率派和贝叶斯派的博客文章，并且在文章中提供了 Python 的示例。该系列的第一篇文章可以在jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/找到。

参数化和非参数化推断方法

在许多情况下，你会基于概率模型来进行分析。该模型描述了你的数据是如何生成的。概率模型并没有实际存在；它仅是一个数学对象，指导你进行分析。一个好的模型是有帮助的，而一个不好的模型可能会误导你。

使用参数化方法时，你假设你的模型属于一个已知的概率分布族。该模型有一个或多个数值的参数，你可以对其进行估计。

使用非参数化模型时，你的模型不需要做出这样的假设。这给你带来了更多的灵活性。然而，这些方法通常实现起来更为复杂，且难以解释。

以下配方分别是参数化和非参数化的：

• 使用最大似然法拟合概率分布配方

• 使用核密度估计非参数化估计概率分布配方

本章仅为你提供了 Python 在统计数据分析方面广泛可能性的一个大致概念。你可以找到许多书籍和在线课程，详细讲解统计方法，例如：

• 维基教科书上的统计学内容：en.wikibooks.org/wiki/Statistics

• 免费的统计学教材可以在 stats.stackexchange.com/questions/170/free-statistical-textbooks 找到

使用 pandas 和 matplotlib 探索数据集

在这第一个食谱中，我们将展示如何使用 pandas 对数据集进行初步分析。这通常是在获得数据后进行的第一步。pandas 让我们能够非常轻松地加载数据、探索变量，并使用 matplotlib 创建基本的图表。

我们将查看一个数据集，该数据集包含了四名网球选手直到 2012 年为止的所有 ATP 比赛。在这里，我们将重点关注罗杰·费德勒。

准备工作

从本书的 GitHub 仓库 github.com/ipython-books/cookbook-data 下载 Tennis 数据集，并将其解压到当前目录。

如何做...

1. 我们导入 NumPy、pandas 和 matplotlib：

   In [1]: import numpy as np
           import pandas as pd
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 数据集是一个 CSV 文件，即一个以逗号分隔值的文本文件。pandas 让我们能够用一个函数加载这个文件：

   In [2]: player = 'Roger Federer'
           filename = "data/{name}.csv".format(
                         name=player.replace(' ', '-'))
           df = pd.read_csv(filename)
   

   我们可以通过在 IPython 笔记本中直接显示数据集来进行首次查看：

   In [3]: df
   Out[3]: Int64Index: 1179 entries, 0 to 1178
           Data columns (total 70 columns):
           year                        1179  non-null values
           tournament                  1179  non-null values
           ...
           player2 total points total  1027  non-null values
           dtypes: float64(49), int64(2), object(19)
   

3. 数据集有很多列。每一行对应罗杰·费德勒的一场比赛。我们来添加一个布尔变量，表示他是否赢得了比赛。tail() 方法显示列的最后几行：

   In [4]: df['win'] = df['winner'] == player
           df['win'].tail()
   Out[4]: 1174    False
           1175     True
           1176     True
           1177     True
           1178    False
           Name: win, dtype: bool
   

4. df['win'] 是一个 Series 对象。它与 NumPy 数组非常相似，只是每个值都有一个索引（这里是比赛索引）。这个对象具有一些标准的统计函数。例如，让我们查看获胜比赛的比例：

   In [5]: print(("{player} has won {vic:.0f}% "
                  "of his ATP matches.").format(
                   player=player, vic=100*df['win'].mean()))
   Roger Federer has won 82% of his ATP matches.
   

5. 现在，我们将观察一些变量随时间的变化。df['start date'] 字段包含了比赛的开始日期，格式为字符串。我们可以使用 pd.to_datetime() 函数将其转换为日期类型：

   In [6]: date = pd.to_datetime(df['start date'])
   

6. 我们现在正在查看每场比赛中的双误比例（考虑到在更长时间的比赛中，双误通常更多！）。这个数字是球员心理状态的一个指标，反映了他的自信水平、在发球时的冒险精神以及其他一些参数。

   In [7]: df['dblfaults'] = (df['player1 double faults'] / 
                              df['player1 total points total'])
   

7. 我们可以使用 head() 和 tail() 方法查看列的开头和结尾，并使用 describe() 获取摘要统计数据。特别地，我们需要注意，有些行包含 NaN 值（即，并不是所有比赛的双误数据都有记录）。

   In [8]: df['dblfaults'].tail()
   Out[8]: 1174    0.018116
           1175    0.000000
           1176    0.000000
           1177    0.011561
           1178         NaN
           Name: dblfaults, dtype: float64
   In [9]: df['dblfaults'].describe()
   Out[9]: count    1027.000000
           mean        0.012129
           std         0.010797
           min         0.000000
           25%         0.004444
           50%         0.010000
           75%         0.018108
           max         0.060606
           dtype: float64
   

8. pandas 中的一个非常强大的功能是groupby()。这个函数允许我们将具有相同列值的行组合在一起。然后，我们可以通过该值对该组进行聚合，计算每个组中的统计数据。例如，下面是我们如何根据比赛场地表面类型来获取胜利比例：

   In [10]: df.groupby('surface')['win'].mean()
   Out[10]: surface
            Indoor: Carpet    0.736842
            Indoor: Clay      0.833333
            Indoor: Hard      0.836283
            Outdoor: Clay     0.779116
            Outdoor: Grass    0.871429
            Outdoor: Hard     0.842324
            Name: win, dtype: float64
   

9. 现在，我们将显示双误差的比例与比赛日期的关系，以及每年的平均值。为此，我们也使用groupby()：

   In [11]: gb = df.groupby('year')
   

10. gb是一个GroupBy实例。它类似于DataFrame对象，但每个组中有多行（每年比赛的所有场次）。我们可以使用mean()操作对这些行进行聚合。我们使用 matplotlib 的plot_date()函数，因为 x 轴包含日期：

    In [12]: plt.plot_date(date, df['dblfaults'], 
                           alpha=.25, lw=0)
             plt.plot_date(gb['start date'].max(), 
                           gb['dblfaults'].mean(), '-', lw=3)
             plt.xlabel('Year')
             plt.ylabel('Proportion of double faults per 
                           match.')
    

    

还有更多...

pandas 是一个用于数据整理和探索性分析的优秀工具。pandas 支持各种格式（文本格式和二进制文件），并允许我们以多种方式操作表格。特别是，groupby()函数非常强大。Wes McKinney 的《Python for Data Analysis》一书对这个库进行了更为详细的讲解。

我们在这里所讲述的仅仅是数据分析过程中的第一步。我们需要更高级的统计方法来获得关于基础现象的可靠信息，做出决策和预测，等等。这个内容将在接下来的章节中讨论。

此外，更复杂的数据集需要更复杂的分析方法。例如，数字记录、图像、声音和视频在应用统计技术之前需要特定的信号处理处理。这些问题将在后续章节中讨论。

开始统计假设检验——一个简单的 z 检验

统计假设检验使我们能够在数据不完全的情况下做出决策。根据定义，这些决策是有不确定性的。统计学家已经开发了严格的方法来评估这种风险。然而，决策过程中总是涉及一些主观性。理论只是帮助我们在不确定的世界中做出决策的工具。

在这里，我们介绍统计假设检验背后的最基本思想。我们将跟随一个极其简单的例子：抛硬币。更准确地说，我们将展示如何进行z 检验，并简要解释其背后的数学思想。这种方法（也称为频率主义方法）尽管在科学中被广泛使用，但也受到了许多批评。稍后我们将展示一种基于贝叶斯理论的更现代的方法。理解这两种方法非常有帮助，因为许多研究和出版物仍然采用频率主义方法。

准备工作

你需要具备基本的概率论知识（随机变量、分布、期望、方差、中心极限定理等）才能理解此方法。

如何实现...

许多频率主义假设检验方法大致包括以下步骤：

1. 写下假设，特别是零假设，它是我们想要证明的假设的对立面（以一定的置信度）。

2. 计算检验统计量，这是一个依赖于检验类型、模型、假设和数据的数学公式。

3. 使用计算得出的值来接受假设、拒绝假设或无法得出结论。

在这个实验中，我们抛掷硬币n次，并观察到h次正面朝上。我们想知道硬币是否公平（零假设）。这个例子非常简单，但对于教学目的非常有用。此外，它是许多更复杂方法的基础。

我们用B(q)表示伯努利分布，其中未知参数为q。您可以访问en.wikipedia.org/wiki/Bernoulli_distribution获取更多信息。

伯努利变量是：

• 0（反面）出现的概率为1-q

• 1（正面朝上）出现的概率为q

以下是进行简单统计z-检验所需的步骤：

1. 假设我们抛掷了n=100次硬币，得到了h=61次正面朝上。我们选择显著性水平为 0.05：硬币是否公平？我们的零假设是：硬币是公平的（q = 1/2）：

   In [1]: import numpy as np
           import scipy.stats as st
           import scipy.special as sp
   In [2]: n = 100  # number of coin flips
           h = 61  # number of heads
           q = .5  # null-hypothesis of fair coin
   

2. 让我们计算z 分数，它由以下公式定义（xbar是分布的估计平均值）。我们将在下一部分*它是如何工作的...*中解释此公式。

   In [3]: xbar = float(h)/n
           z = (xbar - q) * np.sqrt(n / (q*(1-q))); z
   Out[3]: 2.1999999999999997
   

3. 现在，根据 z 分数，我们可以按以下方式计算 p 值：

   In [4]: pval = 2 * (1 - st.norm.cdf(z)); pval
   Out[4]: 0.02780689502699718
   

4. 该 p 值小于 0.05，因此我们拒绝零假设，并得出结论：硬币可能不公平。

它是如何工作的...

投币实验被建模为一系列独立的随机变量，，遵循伯努利分布B(q)。每个*x[i]*代表一次投币实验。经过我们的实验后，我们获得这些变量的实际值（样本）。有时为了区分随机变量（概率对象）和实际值（样本），会使用不同的符号表示。

以下公式给出了样本均值（这里是正面朝上的比例）：

知道分布*B(q)*的期望值和方差，我们计算：

z 检验是的标准化版本（我们去除其均值，并除以标准差，因此我们得到一个均值为 0，标准差为 1 的变量）：

在零假设下，获得大于某个数量z[0]的 z 检验的概率是多少？这个概率称为（双尾）p 值。根据中心极限定理，当n较大时，z 检验大致服从标准正态分布N(0,1)，因此我们得到：

下图展示了 z 分数和 p 值：

z 分数和 p 值的示意图

在这个公式中， 是标准正态分布的累积分布函数。在 SciPy 中，我们可以通过 scipy.stats.norm.cdf 获取它。因此，给定从数据中计算得出的 z 检验，我们计算 p 值：在原假设下，观察到比所观察到的检验值更极端的 z 检验的概率。

如果 p 值小于 5%（这是一个常用的显著性水平，出于任意和历史原因），我们得出结论：

• 原假设为假，因此我们得出结论：硬币是不公平的。

• 原假设为真，如果我们得到了这些值，那就是运气不好。我们无法得出结论。

在这个框架中，我们无法对这两个选项进行歧义消解，但通常选择第一个选项。我们达到了频率主义统计的极限，尽管有一些方法可以缓解这个问题（例如，通过进行几项独立研究并查看它们的所有结论）。

还有更多内容...

存在许多遵循这种模式的统计检验。回顾所有这些检验远超出本书的范围，但你可以查看 en.wikipedia.org/wiki/Statistical_hypothesis_testing 中的参考资料。

由于 p 值不容易解释，它可能导致错误的结论，即使是在同行评审的科学出版物中也是如此。关于该主题的深入探讨，请参见 www.refsmmat.com/statistics/。

另见

• 入门贝叶斯方法部分

入门贝叶斯方法

在上一个例子中，我们使用了频率主义方法对不完全数据进行假设检验。在这里，我们将看到一种基于贝叶斯理论的替代方法。主要思想是认为未知参数是随机变量，就像描述实验的变量一样。关于这些参数的先验知识被整合到模型中。随着数据的不断增加，这些知识会被更新。

频率主义者和贝叶斯主义者对概率的解释不同。频率主义者将概率解释为样本数量趋于无穷大时频率的极限。而贝叶斯主义者则将其解释为一种信念；随着观察到的数据越来越多，这种信念会不断更新。

这里，我们以贝叶斯方法重新审视之前的抛硬币例子。这个例子足够简单，可以进行分析处理。通常，如我们将在本章后面看到的，无法得到解析结果，数值方法变得至关重要。

准备开始

这是一个数学密集的过程。建议具备基本概率论（随机变量、分布、贝叶斯公式）和微积分（导数、积分）的知识。我们使用与前一个方法相同的符号。

如何做...

让q表示获得正面的概率。而在之前的公式中，q只是一个固定的数字，在这里我们认为它是一个随机变量。最初，这个变量遵循一种被称为先验概率分布的分布。它表示我们在开始抛硬币之前，对q的知识。我们将在每次试验后更新这个分布（即后验分布）。

1. 首先，我们假设q是区间[0, 1]上的一个均匀随机变量。这是我们的先验分布：对于所有q，P(q)=1。

2. 然后，我们抛硬币n次。我们用x[i]表示第i次抛掷的结果（0 代表反面，1 代表正面）。

3. 知道观察到的结果*x[i]*后，q的概率分布是多少？贝叶斯定理使我们能够解析地计算出后验分布（数学细节见下一部分）：

4. 我们根据之前的数学公式定义后验分布。我们注意到，这个表达式是*(n+1)*倍的概率质量函数（PMF）二项分布的形式，二项分布可以直接通过scipy.stats获得。（有关二项分布的更多信息，请参见en.wikipedia.org/wiki/Binomial_distribution。）

   In [1]: import numpy as np
           import scipy.stats as st
           import matplotlib.pyplot as plt
           %matplotlib inline
   In [2]: posterior = lambda n, h, q: ((n+1) * 
                                        st.binom(n, q).pmf(h))
   

5. 让我们为观察到h=61次正面朝上的结果和n=100次抛掷总数，绘制这个分布：

   In [3]: n = 100
           h = 61
           q = np.linspace(0., 1., 1000)
           d = posterior(n, h, q)
   In [4]: plt.plot(q, d, '-k')
           plt.ylim(0, d.max()+1)
   

   

   这个曲线表示我们在观察到 61 次正面朝上的结果后，对参数q的信念。

它是如何工作的...

在这一部分，我们解释了贝叶斯定理，并给出了这个例子背后的数学细节。

贝叶斯定理

数据科学中有一个非常通用的思想，那就是用数学模型来解释数据。这通过一个单向过程模型 → 数据来形式化。

一旦这个过程被形式化，数据科学家的任务就是利用数据恢复关于模型的信息。换句话说，我们希望逆转原始过程并得到数据 → 模型。

在一个概率设置中，直接过程表示为条件概率分布 P(data|model)。这是在模型完全指定的情况下，观察到数据的概率。

类似地，逆向过程是P(model|data)。它在知道观察结果（我们拥有的数据）的基础上，给我们关于模型（我们所寻找的）的信息。

贝叶斯定理是一个通用框架的核心，用于逆转模型 → 数据的概率过程。它可以表述如下：

这个方程为我们提供了关于模型的信息，前提是我们知道观测数据。贝叶斯方程广泛应用于信号处理、统计学、机器学习、逆问题以及许多其他科学领域。

在贝叶斯方程中，P(model) 反映了我们关于模型的先验知识。同时，P(data) 是数据的分布。通常它表示为 P(data|model)P(model) 的积分。

总之，贝叶斯方程为我们提供了数据推断的总体路线图：

1. 为直接过程 model → data 指定一个数学模型（P(data|model) 项）。

2. 为模型指定一个先验概率分布 (P(model) 项)。

3. 执行解析或数值计算以求解这个方程。

后验分布的计算

在这个例子中，我们通过以下方程（直接源自贝叶斯定理）找到了后验分布：

由于 x[i] 是独立的，我们得到（h 为正面次数）：

此外，我们可以解析计算以下积分（使用分部积分法和归纳法）：

最后，我们得到：

最大后验估计

我们可以从后验分布中得到一个点估计。例如，最大后验估计 (MAP) 是通过考虑后验分布的 最大 值来作为 q 的估计。我们可以通过解析方法或数值方法找到这个最大值。有关 MAP 的更多信息，请参阅 en.wikipedia.org/wiki/Maximum_a_posteriori_estimation。

在这里，我们可以通过对 q 求导直接推导出后验分布，得到这个估计（假设 1 h * n-1*）：

当 q = h/n 时，该表达式为零。这就是参数 q 的 MAP 估计。这个值恰好是实验中获得的正面比例。

还有更多...

在这个例子中，我们展示了贝叶斯理论中的一些基本概念，并通过一个简单的例子进行说明。我们能够解析推导出后验分布在现实应用中并不常见。尽管如此，这个例子仍然具有启发性，因为它解释了我们接下来将看到的复杂数值方法背后的核心数学思想。

可信区间

后验分布表示在给定观察值的情况下 q 的合理取值。我们可以利用它推导出 可信区间，它很可能包含实际值。可信区间是贝叶斯统计中的类比于频率统计中的置信区间。有关可信区间的更多信息，请参阅 en.wikipedia.org/wiki/Credible_interval。

共轭分布

在这道菜谱中，先验分布和后验分布是共轭的，意味着它们属于同一个分布族（即贝塔分布）。因此，我们能够解析计算后验分布。你可以在en.wikipedia.org/wiki/Conjugate_prior找到有关共轭分布的更多细节。

非信息性（客观）先验分布

我们选择了均匀分布作为未知参数q的先验分布。这是一个简单的选择，它使得计算变得可处理。它反映了一个直观的事实，即我们在先验上并不偏向任何特定的值。然而，也有一些严格的选择完全无信息性先验的方法（参见en.wikipedia.org/wiki/Prior_probability#Uninformative_priors）。一个例子是 Jeffreys 先验，基于这样的思想：先验分布不应依赖于参数化的选择。更多关于 Jeffreys 先验的信息，请参阅en.wikipedia.org/wiki/Jeffreys_prior。在我们的例子中，Jeffreys 先验为：

另见

• 通过马尔科夫链蒙特卡罗方法从后验分布中抽样拟合贝叶斯模型菜谱

使用列联表和卡方检验估计两个变量之间的相关性

而单变量方法处理的是单一变量的观测数据，多变量方法则考虑包含多个特征的观测数据。多变量数据集允许研究变量之间的关系，特别是它们之间的相关性或独立性。

在本菜谱中，我们将查看本章第一道菜谱中的相同网球数据集。采用频率学派方法，我们将估计发球得分数与网球选手获胜的点数比例之间的相关性。

准备工作

从本书的 GitHub 仓库下载网球数据集，链接为github.com/ipython-books/cookbook-data，并将其解压到当前目录。

如何操作...

1. 让我们导入 NumPy、pandas、SciPy.stats 和 matplotlib：

   In [1]: import numpy as np
           import pandas as pd
           import scipy.stats as st
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们加载对应于罗杰·费德勒的数据集：

   In [2]: player = 'Roger Federer'
           filename = "data/{name}.csv".format(
                         name=player.replace(' ', '-'))
           df = pd.read_csv(filename)
   

3. 每行对应一场比赛，70 列包含该比赛中许多选手的特征：

   In [3]: print("Number of columns: " + str(len(df.columns)))
           df[df.columns[:4]].tail()
   Number of columns: 70
                 year                  tournament  start date
           1174  2012  Australian Open, Australia  16.01.2012
           1175  2012                 Doha, Qatar  02.01.2012
           1176  2012                 Doha, Qatar  02.01.2012
           1177  2012                 Doha, Qatar  02.01.2012
           1178  2012                 Doha, Qatar  02.01.2012
   

4. 在这里，我们仅关注获胜点数的比例和（相对的）发球得分数：

   In [4]: npoints = df['player1 total points total']
           points = df['player1 total points won'] / npoints
           aces = df['player1 aces'] / npoints
   In [5]: plt.plot(points, aces, '.')
           plt.xlabel('% of points won')
           plt.ylabel('% of aces')
           plt.xlim(0., 1.)
           plt.ylim(0.)
   

   

   如果这两个变量是独立的，我们就不会在点云中看到任何趋势。在这个图上，稍微有点难以看出。让我们使用 pandas 计算一个相关系数。

5. 我们创建一个新的DataFrame对象，仅包含这些字段（请注意，这一步不是强制性的）。我们还删除了缺失某个字段的行（使用dropna()）：

   In [6]: df_bis = pd.DataFrame({'points': points,
                                  'aces': aces}).dropna()
           df_bis.tail()
   Out[6]:           aces    points
           1173  0.024390  0.585366
           1174  0.039855  0.471014
           1175  0.046512  0.639535
           1176  0.020202  0.606061
           1177  0.069364  0.531792
   

6. 让我们计算比赛中 A 球的相对数量与赢得的点数之间的皮尔逊相关系数：

   In [7]: df_bis.corr()
   Out[7]:             aces    points
           aces    1.000000  0.255457
           points  0.255457  1.000000
   

   约为 0.26 的相关性似乎表明我们的两个变量之间存在正相关关系。换句话说，在一场比赛中 A 球越多，球员赢得的点数就越多（这并不令人惊讶！）。

7. 现在，为了确定变量之间是否存在统计上显著的相关性，我们使用卡方检验来检验列联表中变量的独立性。

8. 首先，我们将变量二值化。在这里，如果球员在比赛中发球 A 球比平常多，值为True，否则为False：

   In [8]: df_bis['result'] = df_bis['points'] > \
                              df_bis['points'].median()
           df_bis['manyaces'] = df_bis['aces'] > \
                              df_bis['aces'].median()
   

9. 然后，我们创建一个列联表，其中包含所有四种可能性（真和真，真和假，依此类推）的频率：

   In [9]: pd.crosstab(df_bis['result'], df_bis['manyaces'])
   Out[9]: manyaces  False  True 
           result                
           False       300    214
           True        214    299
   

10. 最后，我们计算卡方检验统计量和相关的 P 值。零假设是变量之间的独立性。SciPy 在scipy.stats.chi2_contingency中实现了这个测试，返回了几个对象。我们对第二个结果感兴趣，即 P 值：

    In [10]: st.chi2_contingency(_)
    Out[10]: (27.809858855369555,
              1.3384233799633629e-07,
              1L,
              array([[ 257.25024343,  256.74975657],
                     [ 256.74975657,  256.25024343]]))
    

    P 值远低于 0.05，因此我们拒绝零假设，并得出结论：在一场比赛中赢得的点数与赢得的 A 球比例之间存在统计学上显著的相关性（对于罗杰·费德勒！）。

提示

如常，相关性并不意味着因果关系。在这里，外部因素很可能影响两个变量。有关更多细节，请参阅en.wikipedia.org/wiki/Correlation_does_not_imply_causation。

工作原理...

我们在这里提供了一些有关本文中使用的统计概念的细节。

皮尔逊相关系数

皮尔逊相关系数衡量了两个随机变量X和Y之间的线性相关性。它是协方差的归一化版本：

可以通过将这个公式中的期望值替换为样本均值，方差替换为样本方差来估计。关于其推断的更多细节可以在en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient找到。

列联表和卡方检验

列联表包含所有组合结果的频率O[ij]，当存在多个随机变量可以取有限数量的值时。在独立性的零假设下，我们可以基于边际和（每行的总和）计算期望频率E[ij]。卡方统计量的定义如下：

当观察足够多时，这个变量大致遵循卡方分布（正态变量平方和的分布）。一旦我们得到了 p 值，就像开始统计假设检验 - 一个简单的 z 检验中所解释的那样，我们可以拒绝或接受独立性的零假设。然后，我们可以得出（或不得出）变量之间存在显著相关性的结论。

还有更多...

还有许多其他类型的卡方检验，即测试统计量遵循卡方分布的测试。这些测试广泛用于测试分布的拟合度，或测试变量的独立性。更多信息可以在以下页面找到：

• SciPy 文档中的 Chi2 测试可在docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html找到

• 在en.wikipedia.org/wiki/Contingency_table介绍的列联表

• 在en.wikipedia.org/wiki/Pearson's_chi-squared_test介绍的卡方检验

参见

• 开始统计假设检验 - 一个简单的 z 检验食谱

使用最大似然方法将概率分布拟合到数据

解释数据集的一个好方法是对其应用概率模型。找到一个合适的模型可能是一项工作。选择模型后，有必要将其与数据进行比较。这就是统计估计的内容。在这个食谱中，我们对心脏移植后存活时间（1967-1974 年研究）的数据集应用最大似然方法。

准备工作

与本章中通常一样，建议具有概率论和实分析背景。此外，您需要 statsmodels 包来检索测试数据集。有关 statsmodels 的更多信息，请参考statsmodels.sourceforge.net。在 Anaconda 上，您可以使用conda install statsmodels命令安装 statsmodel。

如何做...

1. statsmodels 是一个用于进行统计数据分析的 Python 包。它还包含我们在尝试新方法时可以使用的真实数据集。在这里，我们加载heart数据集：

   In [1]: import numpy as np
           import scipy.stats as st
           import statsmodels.datasets as ds
           import matplotlib.pyplot as plt
           %matplotlib inline
   In [2]: data = ds.heart.load_pandas().data
   

2. 让我们来看看这个DataFrame：

   In [3]: data.tail()
   Out[3]:     survival  censors   age
           64        14        1  40.3
           65       167        0  26.7
           66       110        0  23.7
           67        13        0  28.9
           68         1        0  35.2
   

   这个数据集包含被审查和未被审查的数据：0 的审查意味着患者在研究结束时仍然存活，因此我们不知道确切的存活时间。我们只知道患者至少存活了指定的天数。为简单起见，我们只保留未被审查的数据（这样我们就引入了对未能在移植后存活很长时间的患者的偏见）：

   In [4]: data = data[data.censors==1]
           survival = data.survival
   

3. 让我们通过绘制原始存活数据和直方图来以图形方式查看数据：

   In [5]: plt.subplot(121)
           plt.plot(sorted(survival)[::-1], 'o')
           plt.xlabel('Patient')
           plt.ylabel('Survival time (days)')
           plt.subplot(122)
           plt.hist(survival, bins=15)
           plt.xlabel('Survival time (days)')
           plt.ylabel('Number of patients')
   

   

4. 我们观察到直方图正在迅速下降。幸运的是，今天的生存率要高得多（5 年后的生存率约为 70%）。让我们尝试将指数分布拟合到数据中（关于指数分布的更多信息可以参见en.wikipedia.org/wiki/Exponential_distribution）。根据这个模型，S（生存天数）是一个带有参数的指数随机变量，观察值*s[i]*是从这个分布中抽样得到的。令样本均值为：指数分布的似然函数如下，按照定义（证明见下节）：

   最大似然估计对于速率参数的定义是：值，它最大化了似然函数。换句话说，这是最大化观察到数据的概率的参数，假设这些观察值是从指数分布中抽样得到的。

   在这里，可以证明当时，似然函数取得最大值，这就是速率参数的最大似然估计。让我们通过数值方法计算这个参数：

   In [6]: smean = survival.mean()
           rate = 1./smean
   

5. 为了将拟合的指数分布与数据进行比较，我们首先需要为 x 轴（天数）生成线性间隔的值：

   In [7]: smax = survival.max()
           days = np.linspace(0., smax, 1000)
           dt = smax / 999\.  # bin size: interval between two
                             # consecutive values in `days`
   

   我们可以使用 SciPy 获得指数分布的概率密度函数。参数是尺度，即估计速率的倒数。

   In [8]: dist_exp = st.expon.pdf(days, scale=1./rate)
   

6. 现在，让我们绘制直方图和得到的分布。我们需要将理论分布重新缩放到直方图上（这取决于箱子大小和数据点的总数）：

   In [9]: nbins = 30
           plt.hist(survival, nbins)
           plt.plot(days, dist_exp*len(survival)*smax/nbins,
                    '-r', lw=3)
   

   

   拟合结果远非完美。我们能够找到最大似然估计的解析公式。在更复杂的情况下，这并不总是可能的。因此，我们可能需要求助于数值方法。SciPy 实际上集成了针对大量分布的数值最大似然程序。在这里，我们使用这种其他方法来估计指数分布的参数。

   In [10]: dist = st.expon
            args = dist.fit(survival); args
   Out[10]: (0.99999999994836486, 222.28880590143666)
   

7. 我们可以使用这些参数进行Kolmogorov-Smirnov 检验，该检验评估分布相对于数据的拟合优度。这个检验是基于数据的经验分布函数与参考分布的累积分布函数（CDF）之间的距离。

   In [11]: st.kstest(survival, dist.cdf, args)
   Out[11]: (0.36199685486406347, 8.6470960143358866e-06)
   

   第二个输出值是 p 值。在这里，它非常低：零假设（即观察到的数据来自具有最大似然速率参数的指数分布）可以被高置信度地拒绝。我们试试另一个分布，Birnbaum-Sanders 分布，它通常用于建模故障时间。（关于 Birnbaum-Sanders 分布的更多信息，请访问en.wikipedia.org/wiki/Birnbaum-Saunders_distribution。）

   In [12]: dist = st.fatiguelife
            args = dist.fit(survival)
            st.kstest(survival, dist.cdf, args)
   Out[12]: (0.18773446101946889, 0.073211497000863268)
   

   这次，p 值是 0.07，因此我们在 5%的置信水平下不会拒绝零假设。绘制结果分布时，我们观察到比指数分布更好的拟合：

   In [13]: dist_fl = dist.pdf(days, *args)
            nbins = 30
            plt.hist(survival, nbins)
            plt.plot(days, dist_exp*len(survival)*smax/nbins,
                     '-r', lw=3, label='exp')
            plt.plot(days, dist_fl*len(survival)*smax/nbins,
                     '-g', lw=3, label='BS')
            plt.xlabel("Survival time (days)")
            plt.ylabel("Number of patients")
            plt.legend()
   

   

它是如何工作的...

在这里，我们给出了计算过程，推导出指数分布的速率参数的最大似然估计：

这里， 是样本均值。在更复杂的情况下，我们将需要数值优化方法，其中的原理是使用标准的数值优化算法来最大化似然函数（请参阅第九章，数值优化）。

为了找到该函数的最大值，我们需要计算它关于 的导数：

因此，这个导数的根是

还有更多内容...

这里列出了一些参考资料：

• 最大似然法在维基百科上的介绍，详情请见en.wikipedia.org/wiki/Maximum_likelihood

• Kolmogorov-Smirnov 检验在维基百科上的介绍，详情请见en.wikipedia.org/wiki/Kolmogorov-Smirnov_test

• 拟合优度在en.wikipedia.org/wiki/Goodness_of_fit上有详细介绍

最大似然法是参数化的：模型属于一个预先指定的参数分布族。在下一个方法中，我们将看到一种基于核的方法，它是非参数化的。

另见

• 通过核密度估计非参数地估计概率分布 的方法

使用核密度估计法非参数地估计概率分布

在之前的方法中，我们应用了参数估计方法。我们有一个统计模型（指数分布）来描述我们的数据，并且我们估计了一个参数（分布的速率）。非参数估计处理那些不属于已知分布族的统计模型。这样，参数空间就是无限维的，而不是有限维的（也就是说，我们估计的是函数而不是数值）。

在这里，我们使用 核密度估计 (KDE) 来估算空间分布的概率密度。我们查看了 1848 到 2013 年期间热带气旋的地理位置，数据由 NOAA（美国国家海洋和大气管理局）提供。

准备工作

从本书的 GitHub 仓库下载 Storms 数据集 github.com/ipython-books/cookbook-data，并将其解压到当前目录。数据来自 www.ncdc.noaa.gov/ibtracs/ind…。

你还需要 matplotlib 的工具包 basemap，可以通过 matplotlib.org/basemap/ 获取。使用 Anaconda，你可以通过 conda install basemap 安装它。Windows 用户还可以在 www.lfd.uci.edu/~gohlke/pyt… 找到安装程序。

如何操作...

1. 让我们导入常用的包。使用高斯核的核密度估计在 SciPy.stats 中有实现：

   In [1]: import numpy as np
           import pandas as pd
           import scipy.stats as st
           import matplotlib.pyplot as plt
           from mpl_toolkits.basemap import Basemap
           %matplotlib inline
   

2. 让我们使用 pandas 打开数据：

   In [2]: df = pd.read_csv(
                      "data/Allstorms.ibtracs_wmo.v03r05.csv")
   

3. 数据集包含了自 1848 年以来大多数风暴的信息。单个风暴可能在多个连续的日子中出现多次。

   In [3]: df[df.columns[[0,1,3,8,9]]].head()
   Out[3]:       Serial_Num  Season Basin  Latitude  Longitude
           0  1848011S09080    1848    SI      -8.6       79.8
           1  1848011S09080    1848    SI      -9.0       78.9
           2  1848011S09080    1848    SI     -10.4       73.2
           3  1848011S09080    1848    SI     -12.8       69.9
           4  1848011S09080    1848    SI     -13.9       68.9
   

4. 我们使用 pandas 的 groupby() 函数获取每个风暴的平均位置：

   In [4]: dfs = df.groupby('Serial_Num')
           pos = dfs[['Latitude', 'Longitude']].mean()
           y, x = pos.values.T
           pos.head()
   Out[4]:                 Latitude  Longitude
           Serial_Num                         
           1848011S09080 -15.918182  71.854545
           1848011S15057 -24.116667  52.016667
           1848061S12075 -20.528571  65.342857
           1851080S15063 -17.325000  55.400000
           1851080S21060 -23.633333  60.200000
   

5. 我们使用 basemap 在地图上显示风暴。这个工具包让我们能够轻松地将地理坐标投影到地图上。

   In [5]: m = Basemap(projection='mill', llcrnrlat=-65,
                       urcrnrlat=85, llcrnrlon=-180,
                       urcrnrlon=180)
           x0, y0 = m(-180, -65)
           x1, y1 = m(180, 85)
           m.drawcoastlines()
           m.fillcontinents(color='#dbc8b2')
           xm, ym = m(x, y)
           m.plot(xm, ym, '.r', alpha=.1)
   

   

6. 为了执行核密度估计，我们将风暴的 x 和 y 坐标堆叠成一个形状为 (2, N) 的数组：

   In [6]: h = np.vstack((xm, ym))
   In [7]: kde = st.gaussian_kde(h)
   

7. gaussian_kde() 函数返回了一个 Python 函数。为了在地图上查看结果，我们需要在覆盖整个地图的二维网格上评估该函数。我们通过 meshgrid() 创建此网格，并将 x 和 y 值传递给 kde 函数。kde 接受一个形状为 (2, N) 的数组作为输入，因此我们需要调整数组的形状：

   In [8]: k = 50
           tx, ty = np.meshgrid(np.linspace(x0, x1, 2*k),
                                np.linspace(y0, y1, k))
           v = kde(np.vstack((tx.ravel(), 
                              ty.ravel()))).reshape((k, 2*k))
   

8. 最后，我们使用 imshow() 显示估算的密度：

   In [9]: m.drawcoastlines()
           m.fillcontinents(color='#dbc8b2')
           xm, ym = m(x, y)
           m.imshow(v, origin='lower', extent=[x0,x1,y0,y1],
                    cmap=plt.get_cmap('Reds'))
   

   

它是如何工作的...

一组 n 点 {x[i]} 的 核密度估计器 表示为：

在这里，h>0 是一个缩放参数（带宽），K(u) 是 核函数，它是一个对称函数，其积分为 1。这个估算器与经典的直方图进行比较，其中核是一个 顶帽 函数（一个取值在 {0,1} 中的矩形函数），但是这些块将位于规则的网格上，而不是数据点上。关于核密度估计器的更多信息，请参考 en.wikipedia.org/wiki/Kernel_density_estimation。

可以选择多个核。这里，我们选择了 高斯核，因此 KDE 是以所有数据点为中心的高斯函数的叠加，它是密度的估算。

带宽的选择并非 trivial（简单的）；它在一个过低的值（小偏差，高方差：过拟合）和一个过高的值（高偏差，小方差：欠拟合）之间存在权衡。我们将在下一章回到这个重要的偏差-方差权衡概念。有关偏差-方差权衡的更多信息，请参考en.wikipedia.org/wiki/Bias-variance_dilemma。

以下图示意了 KDE。数据集包含四个位于* [0,1] *的点（黑线）。估计的密度是平滑曲线，这里使用了多个带宽值表示。

核密度估计

小贴士

在 statsmodels 和 scikit-learn 中有其他的 KDE 实现。你可以在jakevdp.github.io/blog/2013/12/01/kernel-density-estimation/找到更多信息。

另见

• 使用最大似然法拟合概率分布到数据食谱

通过从后验分布中采样，使用马尔可夫链蒙特卡罗方法拟合贝叶斯模型

在本食谱中，我们展示了一种非常常见且有用的贝叶斯模型后验分布特征化方法。假设你有一些数据，想要获取关于潜在随机现象的信息。在频率学派的方法中，你可以尝试在给定的分布族中拟合一个概率分布，使用类似最大似然法的参数化方法。优化过程将得出最大化观察数据的概率的参数，假设为零假设。

在贝叶斯方法中，你将参数本身视为随机变量。它们的先验分布反映了你对这些参数的初始知识。在观察后，你的知识得到更新，并在参数的后验分布中体现出来。

贝叶斯推断的一个典型目标是特征化后验分布。贝叶斯定理提供了一种分析方法来实现这一目标，但由于模型的复杂性和维度的数量，它在实际问题中通常不切实际。马尔可夫链 蒙特卡罗方法，例如Metropolis-Hastings 算法，提供了一种数值方法来逼近后验分布。

这里，我们介绍了PyMC包，它提供了一种有效且自然的接口，用于在贝叶斯框架中拟合数据的概率模型。我们将使用来自美国国家海洋和大气管理局（NOAA）的数据，研究自 1850 年代以来北大西洋地区风暴的年频率。

本食谱主要受 PyMC 官网教程的启发（请参见*更多内容...*部分的链接）。

准备工作

你可以在该软件包的网站上找到安装 PyMC 的说明。在本示例中，我们将使用 PyMC2。新版本（PyMC3）在编写时仍在开发中，可能会有显著差异。有关 PyMC 的更多信息，请参阅pymc-devs.github.io/pymc/。对于 Anaconda 用户，可以尝试conda install -c https://conda.binstar.org/pymc pymc。Windows 用户也可以在www.lfd.uci.edu/~gohlke/pyt…找到安装程序。

你还需要从书籍的 GitHub 仓库下载风暴数据集，链接为github.com/ipython-books/cookbook-data，并将其解压到当前目录。

如何操作...

1. 让我们导入标准包和 PyMC：

   In [1]: import numpy as np
           import pandas as pd
           import pymc
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 让我们使用 pandas 导入数据：

   In [2]: df = pd.read_csv(
                   "data/Allstorms.ibtracs_wmo.v03r05.csv",
                   delim_whitespace=False)
   

3. 使用 pandas，只需一行代码就能得到北大西洋地区的年风暴数。我们首先选择该海域的风暴（NA），然后按年份（Season）对行进行分组，再计算独特风暴的数量（Serial_Num），因为每个风暴可能跨越几天（使用nunique()方法）：

   In [3]: cnt = df[df['Basin'] == ' NA'].groupby('Season') \
                              ['Serial_Num'].nunique()
           years = cnt.index
           y0, y1 = years[0], years[-1]
           arr = cnt.values
           plt.plot(years, arr, '-ok')
           plt.xlim(y0, y1)
           plt.xlabel("Year")
           plt.ylabel("Number of storms")
   

   

4. 现在，我们定义我们的概率模型。我们假设风暴遵循时间依赖的泊松过程，并且具有一个确定性的速率。我们假设该速率是一个分段常数函数，在切换点switchpoint之前取值early_mean，在切换点之后取值late_mean。这三个未知参数被视为随机变量（我们将在*它是如何工作的…*部分中详细描述）。

   提示

   泊松过程（en.wikipedia.org/wiki/Poisson_process）是一种特殊的点过程，即描述瞬时事件随机发生的随机过程。泊松过程是完全随机的：事件以给定速率独立发生。另见第十三章，随机动力学系统。

   In [4]: switchpoint = pymc.DiscreteUniform('switchpoint',
                                              lower=0, 
                                              upper=len(arr))
           early_mean = pymc.Exponential('early_mean', beta=1)
           late_mean = pymc.Exponential('late_mean', beta=1)
   

5. 我们将分段常数速率定义为一个 Python 函数：

   In [5]: @pymc.deterministic(plot=False)
           def rate(s=switchpoint, e=early_mean, l=late_mean):
               out = np.empty(len(arr))
               out[:s] = e
               out[s:] = l
               return out
   

6. 最后，观察变量是年风暴数。它遵循一个具有随机均值的泊松变量（底层泊松过程的速率）。这是泊松过程的一个已知数学性质。

   In [6]: storms = pymc.Poisson('storms', mu=rate, value=arr, 
                                 observed=True)
   

7. 现在，我们使用 MCMC 方法从后验分布中进行采样，给定观察数据。sample()方法启动拟合的迭代过程：

   In [7]: model = pymc.Model([switchpoint, early_mean, 
                               late_mean,
                               rate, storms])
   In [8]: mcmc = pymc.MCMC(model)
           mcmc.sample(iter=10000, burn=1000, thin=10)
            [----       17%            ] 1774 of 10000 complete
            [-----------100%-----------] 10000 of 10000 complete
   

8. 让我们绘制采样的马尔科夫链。它们的平稳分布对应我们想要表征的后验分布。

   In [9]: plt.subplot(311)
           plt.plot(mcmc.trace('switchpoint')[:])
           plt.ylabel("Switch point")
           plt.subplot(312)
           plt.plot(mcmc.trace('early_mean')[:])
           plt.ylabel("Early mean")
           plt.subplot(313)
           plt.plot(mcmc.trace('late_mean')[:])
           plt.xlabel("Iteration")
           plt.ylabel("Late mean")
   

   

9. 我们还绘制了样本的分布，这些样本对应我们参数的后验分布，数据点已被考虑在内：

   In [10]: plt.subplot(131)
            plt.hist(mcmc.trace('switchpoint')[:] + y0, 15)
            plt.xlabel("Switch point")
            plt.ylabel("Distribution")
            plt.subplot(132)
            plt.hist(mcmc.trace('early_mean')[:], 15)
            plt.xlabel("Early mean")
            plt.subplot(133)
            plt.hist(mcmc.trace('late_mean')[:], 15)
            plt.xlabel("Late mean")
   

   

10. 通过这些分布的样本均值，我们得到三个未知参数的后验估计，包括风暴频率突然增加的年份：

    In [11]: yp = y0 + mcmc.trace('switchpoint')[:].mean()
             em = mcmc.trace('early_mean')[:].mean()
             lm = mcmc.trace('late_mean')[:].mean()
             print((yp, em, lm))
    (1966.681111111111, 8.2843072252292682, 16.728831395584947)
    

11. 现在，我们可以将估计的比率绘制在观察数据之上：

    In [12]: plt.plot(years, arr, '-ok')
             plt.axvline(yp, color='k', ls='--')
             plt.plot([y0, yp], [em, em], '-b', lw=3)
             plt.plot([yp, y1], [lm, lm], '-r', lw=3)
             plt.xlim(y0, y1)
             plt.xlabel("Year")
             plt.ylabel("Number of storms")
    

    

它是如何工作的……

一般思路是定义一个贝叶斯概率模型，并将其拟合到数据中。该模型可能是估计或决策任务的起点。模型本质上由通过有向无环图（DAG）连接的随机或确定性变量描述。A与B相连，如果B完全或部分由A决定。下图展示了本示例中使用的模型图：

In [13]: graph = pymc.graph.graph(model)
         from IPython.display import display_png
         display_png(graph.create_png(), raw=True)

提示

如你所见，PyMC 可以创建模型的图形表示。你需要安装 GraphViz（参考www.graphviz.org）、pydot 和 pyparsing。由于一个不幸的 bug，你可能需要安装 pyparsing 的特定版本：

pip install pyparsing==1.5.7
pip install pydot

随机变量遵循可以通过模型中的固定数字或其他变量进行参数化的分布。参数本身也可以是随机变量，反映了观察之前的知识。这是贝叶斯建模的核心。

分析的目标是将观察结果纳入模型中，以便随着更多数据的可用，更新我们的知识。尽管贝叶斯定理为我们提供了计算后验分布的精确方法，但在现实世界的问题中很少能实际应用。这主要是由于模型的复杂性。为了应对这一问题，已经开发了数值方法。

这里使用的马尔科夫链蒙特卡洛（MCMC）方法使我们能够通过模拟一个具有所需分布作为平衡分布的马尔科夫链，从复杂的分布中采样。梅特罗波利斯-哈斯廷斯算法是这种方法在当前示例中的具体应用。

该算法在 PyMC 的MCMC类中实现。burn参数决定丢弃多少初始迭代次数。这是必要的，因为马尔科夫链需要经过若干次迭代才能收敛到其平衡分布。thin参数对应于在评估分布时跳过的步数，以尽量减少样本的自相关性。你可以在pymc-devs.github.io/pymc/modelfitting.html找到更多信息。

还有更多内容……

这里有一些参考资料：

• 我们大量借鉴的一个优秀的 PyMC 教程可以在pymc-devs.github.io/pymc/tutorial.html找到。

• 由 Cameron Davidson-Pilon 撰写的关于该主题的必读免费电子书，完全使用 IPython 笔记本编写，可通过camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/获取。

• 在en.wikipedia.org/wiki/Markov_chain_Monte_Carlo介绍的马尔可夫链蒙特卡洛方法

• 在en.wikipedia.org/wiki/Metropolis-Hastings_algorithm介绍的 Metropolis-Hastings 算法

另见

• 贝叶斯方法入门配方

在 IPython 笔记本中使用 R 编程语言分析数据

R (www.r-project.org)是一种免费的特定领域编程语言，用于统计学。其语法非常适合统计建模和数据分析。相比之下，Python 的语法通常更适用于通用编程。幸运的是，IPython 使你可以同时享受两者的优势。例如，你可以在普通的 IPython 笔记本中随时插入 R 代码片段。你可以继续使用 Python 和 pandas 进行数据加载和整理，然后切换到 R 来设计和拟合统计模型。使用 R 代替 Python 进行这些任务，不仅仅是语法问题；R 自带一个令人印象深刻的统计工具箱，目前 Python 尚无匹敌。

在这个配方中，我们将展示如何在 IPython 中使用 R，并通过一个简单的数据分析示例来展示 R 的最基本功能。

准备工作

你需要 statsmodels 包来进行此配方的操作。你可以在之前的配方中找到安装说明，使用最大似然方法拟合概率分布到数据。

你还需要安装 R。使用 R 从 IPython 的步骤有三个。首先，安装 R 和 rpy2（R 到 Python 的接口）。当然，这一步只需要做一次。然后，为了在 IPython 会话中使用 R，你需要加载 IPython 的 R 扩展。

1. 从cran.r-project.org/mirrors.html下载适合你操作系统的 R 并进行安装。在 Ubuntu 上，你可以执行sudo apt-get install r-base-dev。

2. 从rpy.sourceforge.net/rpy2.html下载 rpy2 并进行安装。在 Linux 上使用 Anaconda 时，你可以尝试conda install -c https://conda.binstar.org/r rpy2。或者，你也可以执行pip install rpy2。

3. 然后，要在 IPython 笔记本中执行 R 代码，首先执行%load_ext rmagic。

   提示

   rpy2 似乎在 Windows 上不太兼容。我们建议使用 Linux 或 Mac OS X。

如何操作...

在这里，我们将使用以下工作流程：首先，从 Python 加载数据。然后，使用 R 设计和拟合模型，并在 IPython 笔记本中绘制一些图表。我们也可以从 R 加载数据，或者使用 Python 的 statsmodels 包设计和拟合统计模型，等等。特别地，我们在这里做的分析本可以完全在 Python 中进行，而无需依赖 R 语言。这个示例仅展示了 R 的基础知识，并说明了 R 和 Python 如何在 IPython 会话中协同工作。

1. 让我们使用 statsmodels 包加载 longley 数据集。这个数据集包含了美国从 1947 年到 1962 年的一些经济指标。我们还加载了 IPython 的 R 扩展：

   In [1]: import statsmodels.datasets as sd
   In [2]: data = sd.longley.load_pandas()
   In [3]: %load_ext rmagic
   

2. 我们将 x 和 y 定义为外生（自变量）和内生（因变量）变量。内生变量量化了国家的总就业量。

   In [4]: data.endog_name, data.exog_name
   Out[4]: ('TOTEMP', ['GNPDEFL', 'GNP', 'UNEMP',
                       'ARMED', 'POP', 'YEAR'])
   In [5]: y, x = data.endog, data.exog
   

3. 为了方便，我们将内生变量添加到 x 数据框中：

   In [6]: x['TOTEMP'] = y
   In [7]: x
   Out[7]:     GNPDEFL     GNP  UNEMP     POP  YEAR  TOTEMP
           0      83.0  234289   2356  107608  1947   60323
           1      88.5  259426   2325  108632  1948   61122
           2      88.2  258054   3682  109773  1949   60171
           ...
           13    114.2  502601   3931  125368  1960   69564
           14    115.7  518173   4806  127852  1961   69331
           15    116.9  554894   4007  130081  1962   70551
   

4. 我们将在 R 中绘制一个简单的图表。首先，我们需要将 Python 变量传递给 R。我们可以使用 %R -i var1,var2 魔法。然后，我们可以调用 R 的 plot() 命令：

   In [8]: gnp = x['GNP']
           totemp = x['TOTEMP']
   In [9]: %R -i totemp,gnp plot(gnp, totemp)
   

   

5. 现在，数据已传递给 R，我们可以对数据进行线性回归。lm() 函数让我们进行线性回归。在这里，我们希望将 totemp（总就业）表示为国家 GNP 的函数：

   In [10]: %%R
            # Least-squares regression
            fit <- lm(totemp ~ gnp);
            # Display the coefficients of the fit.
            print(fit$coefficients)
            plot(gnp, totemp)  # Plot the data points.
            abline(fit)  # And plot the linear regression.
    (Intercept)          gnp 
   5.184359e+04 3.475229e-02
   

   

它是如何工作的...

%R 魔法中的 -i 和 -o 选项允许我们在 IPython 和 R 之间传递变量。变量名称需要用逗号分隔。你可以在 rpy.sourceforge.net/ 的文档中找到有关 %R 魔法的更多信息。

在 R 中，波浪符号 (~) 表示因变量对一个或多个自变量的依赖关系。lm() 函数允许我们对数据拟合一个简单的线性回归模型。在这里，totemp 被表示为 gnp 的函数：

在这里，b（截距）和 a 是线性回归模型的系数。这两个值由 fit$coefficients 在 R 中返回，其中 fit 是拟合的模型。

当然，我们的数据点并不完全满足这个关系。系数的选择是为了最小化这个线性预测与实际值之间的误差。这通常通过最小化以下最小二乘误差来完成：

数据点是 (gnp[i], totemp[i])。由 lm() 返回的系数 a 和 b 使得这个和最小化：它们最适合数据。

还有更多...

回归分析是一个重要的统计概念，我们将在下一章中详细讨论。以下是一些参考资料：

• 维基百科上的回归分析，链接：en.wikipedia.org/wiki/Regression_analysis

• 维基百科上的最小二乘法，网址为 en.wikipedia.org/wiki/Linear…

R 是一个非常适合高级统计分析的平台。虽然 Python 有一些统计包，如 pandas 和 statsmodels，能够实现许多常见功能，但目前 R 提供的统计工具箱数量仍然无可匹敌。然而，Python 在统计学以外有着更广泛的应用范围，并且是一个出色的通用编程语言，配备了大量不同的包。

由于 IPython 支持多种语言，你不必在这些语言之间做出选择。你可以继续使用 Python，在需要 Python 仍未涵盖的高度特定统计功能时，切换到 R。

这里有一些关于 R 的参考资料：

• R 介绍，网址为 cran.r-project.org/doc/manuals/R-intro.html

• R 教程，网址为 www.cyclismo.org/tutorial/R/

• CRAN（综合 R 存档网络），包含许多 R 包，网址为 cran.r-project.org

• 可以在 nbviewer.ipython.org/github/ipython/ipython/blob/master/examples/Builtin%20Extensions/R%20Magics.ipynb 查阅 IPython 和 R 教程

另请参见

• 使用 pandas 和 matplotlib 探索数据集 方案

第八章：机器学习

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

• 开始使用 scikit-learn

• 使用逻辑回归预测谁能在泰坦尼克号上生还

• 使用 K 近邻分类器学习识别手写数字

• 从文本中学习 – 使用朴素贝叶斯进行自然语言处理

• 使用支持向量机进行分类任务

• 使用随机森林选择回归任务中的重要特征

• 使用主成分分析（PCA）降低数据集的维度

• 使用聚类方法发现数据集中的隐藏结构

介绍

在上一章中，我们关注的是通过部分观察获得数据的洞察、理解复杂现象，并在不确定性面前做出明智决策。在这里，我们仍然关注使用统计工具分析和处理数据。然而，目标不一定是理解数据，而是从数据中学习。

从数据中学习类似于我们人类的学习方式。从经验中，我们直觉地学习世界的普遍事实和关系，尽管我们未必完全理解其复杂性。随着计算机计算能力的增强，它们也能从数据中学习。这就是机器学习的核心，机器学习是人工智能、计算机科学、统计学和应用数学中的一个现代而迷人的分支。有关机器学习的更多信息，请参考 en.wikipedia.org/wiki/Machine_learning。

本章是对机器学习一些最基础方法的实践介绍。这些方法是数据科学家日常使用的。我们将使用scikit-learn，一个流行且用户友好的 Python 机器学习包，来运用这些方法。

一些词汇

在本介绍中，我们将解释机器学习的基本定义和概念。

从数据中学习

在机器学习中，大部分数据可以表示为一个数值表格。每一行称为观察值、样本或数据点。每一列称为特征或变量。

设N为行数（或数据点的数量），D为列数（或特征的数量）。数字 D 也被称为数据的维度。原因是我们可以将这个表格视为一个在 D 维空间中的向量集 E（或向量空间）。在这里，一个向量 x 包含 D 个数字 (x[1], ..., x[D])，也叫做分量。这种数学视角非常有用，我们将在本章中持续使用它。

我们通常区分监督学习和无监督学习：

• 监督学习是指我们有一个与数据点x相关联的标签y。目标是从我们的数据中学习从x到y的映射。数据为我们提供了一个有限点集的映射，但我们希望将这个映射泛化到整个集合E。

• 无监督学习是指我们没有任何标签。我们希望做的是发现数据中的某种隐藏结构。

监督学习

从数学角度看，监督学习包括找到一个函数f，将点集E映射到标签集F，并且已知一个有限的关联集*(x, y)，这些数据来自我们的数据集。这就是泛化的意义所在：在观察了对(x[i], y[i])后，给定一个新的x*，我们能够通过将函数f应用到x上来找到相应的y。关于监督学习的更多信息，请参阅en.wikipedia.org/wiki/Supervised_learning。

通常做法是将数据集分为两个子集：训练集和测试集。我们在训练集上学习函数f，并在测试集上进行测试。这对于评估模型的预测能力至关重要。如果在同一数据集上训练和测试模型，我们的模型可能无法很好地泛化。这就是过拟合的基本概念，我们将在本章后面详细讲解。

我们通常区分分类和回归，它们是监督学习的两种特定实例。

分类是指我们的标签y只能取有限的值（类别）。例如：

• 手写数字识别：x是一个包含手写数字的图像；y是一个 0 到 9 之间的数字

• 垃圾邮件过滤：x是电子邮件，y是 1 或 0，取决于该电子邮件是否是垃圾邮件

回归是指我们的标签y可以取任何实数（连续值）。例如：

• 股票市场数据预测

• 销售预测

• 从图片中检测一个人的年龄

分类任务将我们的空间E划分为不同的区域（也称为划分），每个区域与标签y的某个特定值相关联。回归任务则产生一个数学模型，将一个实数与空间E中的任何点x关联。这一差异可以通过下图来说明：

分类与回归的区别

分类与回归可以结合使用。例如，在probit 模型中，尽管因变量是二元的（分类），但该变量属于某一类别的概率也可以通过回归来建模。我们将在关于逻辑回归的示例中看到这一点。关于 probit 模型的更多信息，请参阅en.wikipedia.org/wiki/Probit_model。

无监督学习

广义而言，非监督学习帮助我们发现数据中的系统性结构。这比监督学习更难理解，因为通常没有明确的问题和答案。欲了解更多非监督学习的信息，请参考en.wikipedia.org/wiki/Unsupervised_learning。

这里有一些与非监督学习相关的重要术语：

• 聚类：将相似的点聚集在簇内

• 密度估计：估计一个概率密度函数，用以解释数据点的分布

• 降维：通过将数据点投影到低维空间，获取高维数据点的简单表示（特别是用于数据可视化）

• 流形学习：找到包含数据点的低维流形（也称为非线性降维）

特征选择和特征提取

在监督学习的背景下，当我们的数据包含许多特征时，有时需要选择其中的一个子集。我们想保留的特征是那些与我们问题最相关的特征。这就是特征选择的问题。

此外，我们可能希望通过对原始数据集应用复杂的变换来提取新特征。这就是特征提取。例如，在计算机视觉中，直接在像素上训练分类器通常不是最有效的方法。我们可能希望提取相关的兴趣点或进行适当的数学变换。这些步骤取决于我们的数据集和我们希望回答的问题。

例如，在使用学习模型之前，通常需要对数据进行预处理。特征缩放（或数据归一化）是一个常见的预处理步骤，其中特征会线性地重新缩放，以适应范围*[-1,1]或[0,1]*。

特征提取和特征选择涉及领域专业知识、直觉和数学方法的平衡组合。这些早期步骤至关重要，可能比学习步骤本身还要重要。原因是与我们问题相关的少数维度通常隐藏在数据集的高维度中。我们需要揭示感兴趣的低维结构，以提高学习模型的效率。

在本章中，我们将看到一些特征选择和特征提取方法。特定于信号、图像或声音的方法将在第十章，信号处理，和第十一章，图像和音频处理中介绍。

这里有一些进一步的参考资料：

• 在 scikit-learn 中的特征选择，文档说明请参见scikit-learn.org/stable/modules/feature_selection.html

• 维基百科上的特征选择 en.wikipedia.org/wiki/Feature_selection

过拟合、欠拟合与偏差-方差权衡

机器学习中的一个核心概念是过拟合与欠拟合之间的权衡。一个模型可能能够准确地表示我们的数据。然而，如果它过于准确，它可能无法很好地推广到未见过的数据。例如，在面部识别中，过于精确的模型可能无法识别当天发型不同的人。原因是我们的模型可能会在训练数据中学习到无关的特征。相反，一个训练不足的模型也无法很好地推广。例如，它可能无法正确识别双胞胎。有关过拟合的更多信息，请参考 en.wikipedia.org/wiki/Overfitting。

一个减少过拟合的常见方法是向模型中添加结构，例如，通过正则化。这种方法在训练过程中倾向于选择更简单的模型（奥卡姆剃刀原则）。你可以在 en.wikipedia.org/wiki/Regularization_%28mathematics%29 找到更多信息。

偏差-方差困境与过拟合和欠拟合问题密切相关。一个模型的偏差量化了它在多个训练集上的准确性。方差量化了模型对训练集小变化的敏感性。一个鲁棒的模型不会对小的变化过于敏感。这个困境涉及到同时最小化偏差和方差；我们希望模型既精确又鲁棒。简单的模型通常不太准确，但更鲁棒；复杂的模型往往更准确，但鲁棒性差。有关偏差-方差困境的更多信息，请参考 en.wikipedia.org/wiki/Bias-variance_dilemma。

这个权衡的重要性无法被过分强调。这个问题贯穿了整个机器学习领域。我们将在本章中看到具体的例子。

模型选择

正如我们在本章中所看到的，存在许多监督学习和无监督学习的算法。例如，本章将讨论一些著名的分类器，包括逻辑回归、最近邻、朴素贝叶斯和支持向量机。还有许多其他算法我们无法在这里讨论。

没有任何模型能在所有情况下都表现得比其他模型更好。一个模型可能在某个数据集上表现良好，而在另一个数据集上表现差。这就是模型选择的问题。

我们将看到一些系统化的方法来评估模型在特定数据集上的质量（特别是交叉验证）。实际上，机器学习并不是一门“精确的科学”，因为它通常涉及试错过程。我们需要尝试不同的模型，并通过经验选择表现最好的那个。

尽管如此，理解学习模型的细节使我们能够直观地了解哪个模型最适合当前问题。

这里有一些关于这个问题的参考资料：

• 维基百科上的模型选择，详见 en.wikipedia.org/wiki/Model_selection

• scikit-learn 文档中的模型评估，详见 scikit-learn.org/stable/modules/model_evaluation.html

• 关于如何选择分类器的博客文章，详见 blog.echen.me/2011/04/27/choosing-a-machine-learning-classifier/

机器学习参考资料

这里有一些优秀的、数学内容较多的机器学习教科书：

• 模式识别与机器学习，Christopher M. Bishop，(2006)，Springer

• 机器学习——一种概率视角，Kevin P. Murphy，(2012)，MIT Press

• 统计学习的元素，Trevor Hastie，Robert Tibshirani，Jerome Friedman，(2009)，Springer

这里有一些更适合没有强数学背景的程序员的书籍：

• 黑客的机器学习，Drew Conway，John Myles White，(2012)，O'Reilly Media

• 机器学习实战，Peter Harrington，(2012)，Manning Publications Co.

你可以在线找到更多的参考资料。

本章中我们未能涵盖的机器学习方法的重要类别包括神经网络和深度学习。深度学习是机器学习中一个非常活跃的研究领域。许多最先进的成果目前都是通过使用深度学习方法实现的。有关深度学习的更多信息，请参见 en.wikipedia.org/wiki/Deep_learning。

开始使用 scikit-learn

在这个食谱中，我们介绍了机器学习 scikit-learn 包的基础知识 (scikit-learn.org)。这个包是我们在本章中将要使用的主要工具。它简洁的 API 使得定义、训练和测试模型变得非常容易。而且，scikit-learn 专为速度和（相对）大数据而设计。

我们将在这里展示一个非常基础的线性回归示例，应用于曲线拟合的背景下。这个玩具示例将帮助我们说明关键概念，如线性模型、过拟合、欠拟合、正则化和交叉验证。

准备工作

你可以在主文档中找到安装 scikit-learn 的所有指令。更多信息，请参考 scikit-learn.org/stable/install.html。使用 anaconda 时，你可以在终端中输入 conda install scikit-learn。

如何操作...

我们将生成一个一维数据集，使用一个简单的模型（包括一些噪声），并尝试拟合一个函数到这些数据上。通过这个函数，我们可以对新的数据点进行预测。这是一个曲线拟合回归问题。

1. 首先，让我们进行所有必要的导入：

   In [1]: import numpy as np
           import scipy.stats as st
           import sklearn.linear_model as lm
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们现在定义一个确定性的非线性函数，作为我们生成模型的基础：

   In [2]: f = lambda x: np.exp(3 * x)
   

3. 我们生成在 [0,2] 范围内的曲线值：

   In [3]: x_tr = np.linspace(0., 2, 200)
           y_tr = f(x_tr)
   

4. 现在，让我们生成在 [0,1] 范围内的数据点。我们使用函数 f 并加入一些高斯噪声：

   In [4]: x = np.array([0, .1, .2, .5, .8, .9, 1])
           y = f(x) + np.random.randn(len(x))
   

5. 让我们在 [0,1] 范围内绘制我们的数据点：

   In [5]: plt.plot(x_tr[:100], y_tr[:100], '--k')
           plt.plot(x, y, 'ok', ms=10)
   

   

   在图像中，虚线曲线表示生成模型。

6. 现在，我们使用 scikit-learn 将线性模型拟合到数据上。这个过程有三个步骤。首先，我们创建模型（LinearRegression 类的一个实例）。然后，我们将模型拟合到我们的数据上。最后，我们从训练好的模型中预测值。

   In [6]: # We create the model.
           lr = lm.LinearRegression()
           # We train the model on our training dataset.
           lr.fit(x[:, np.newaxis], y)
           # Now, we predict points with our trained model.
           y_lr = lr.predict(x_tr[:, np.newaxis])
   

   注意

   我们需要将x和x_tr转换为列向量，因为在 scikit-learn 中，一般约定观测值是行，特征是列。在这里，我们有七个观测值，每个观测值有一个特征。

7. 我们现在绘制训练后的线性模型结果。这里我们得到了一条绿色的回归线：

   In [7]: plt.plot(x_tr, y_tr, '--k')
           plt.plot(x_tr, y_lr, 'g')
           plt.plot(x, y, 'ok', ms=10)
           plt.xlim(0, 1)
           plt.ylim(y.min()-1, y.max()+1)
           plt.title("Linear regression")
   

   

8. 线性拟合在这里并不适用，因为数据点是根据非线性模型（指数曲线）生成的。因此，我们现在将拟合一个非线性模型。更准确地说，我们将为数据点拟合一个多项式函数。我们仍然可以使用线性回归来做到这一点，方法是预先计算数据点的指数。这是通过生成范德蒙德矩阵来完成的，使用的是np.vander函数。我们将在*它是如何工作的...*部分解释这个技巧。在下面的代码中，我们执行并绘制拟合结果：

   In [8]: lrp = lm.LinearRegression()
           plt.plot(x_tr, y_tr, '--k')
           for deg in [2, 5]:
               lrp.fit(np.vander(x, deg + 1), y)
               y_lrp = lrp.predict(np.vander(x_tr, deg + 1))
               plt.plot(x_tr, y_lrp,
                        label='degree ' + str(deg))
               plt.legend(loc=2)
               plt.xlim(0, 1.4)
               plt.ylim(-10, 40)
               # Print the model's coefficients.
               print(' '.join(['%.2f' % c for c in 
                               lrp.coef_]))
           plt.plot(x, y, 'ok', ms=10)
           plt.title("Linear regression")
   25.00 -8.57 0.00
   -132.71 296.80 -211.76 72.80 -8.68 0.00
   

   

   我们拟合了两个多项式模型，分别是 2 次和 5 次的多项式。2 次多项式似乎拟合数据点的精度不如 5 次多项式。然而，它似乎更稳健；5 次多项式在预测数据点外的值时表现得非常差（例如，可以查看 x * 1* 部分）。这就是我们所说的过拟合；通过使用一个过于复杂的模型，我们在训练数据集上获得了更好的拟合，但在数据集外的模型却不够稳健。

   注意

   注意 5 次多项式的系数非常大；这通常是过拟合的一个标志。

9. 我们现在将使用一种不同的学习模型，称为岭回归。它与线性回归类似，不同之处在于它防止多项式系数变得过大。这正是前一个例子中发生的情况。通过在损失函数中添加正则化****项，岭回归对基础模型施加了一些结构。我们将在下一节中看到更多细节。

   岭回归模型有一个超参数，表示正则化项的权重。我们可以使用 Ridge 类通过反复试验尝试不同的值。然而，scikit-learn 提供了另一个叫做 RidgeCV 的模型，它包括 交叉验证的参数搜索。实际上，这意味着我们无需手动调整该参数——scikit-learn 会为我们完成。由于 scikit-learn 的模型始终遵循 fit-predict API，我们只需在之前的代码中将 lm.LinearRegression() 替换为 lm.RidgeCV()。我们将在下一节中提供更多细节。

   In [9]: ridge = lm.RidgeCV()
           plt.plot(x_tr, y_tr, '--k')
   
           for deg in [2, 5]:
               ridge.fit(np.vander(x, deg + 1), y);
               y_ridge = ridge.predict(np.vander(x_tr, deg+1))
               plt.plot(x_tr, y_ridge,
                        label='degree ' + str(deg))
               plt.legend(loc=2)
               plt.xlim(0, 1.5)
               plt.ylim(-5, 80)
               # Print the model's coefficients.
               print(' '.join(['%.2f' % c 
                               for c in ridge.coef_]))
   
           plt.plot(x, y, 'ok', ms=10)
           plt.title("Ridge regression")
   11.36 4.61 0.00
   2.84 3.54 4.09 4.14 2.67 0.00
   

   

   这次，5 次多项式似乎比简单的 2 次多项式更加精确（后者现在导致 欠拟合）。岭回归在此处缓解了过拟合问题。观察 5 次多项式的系数比前一个例子中要小得多。

它是如何工作的...

本节我们将解释本食谱中涉及的所有方面。

scikit-learn API

scikit-learn 为监督学习和无监督学习实现了一个简洁且一致的 API。我们的数据点应该存储在一个 (N,D) 矩阵 X 中，其中 N 是观测值的数量，D 是特征的数量。换句话说，每一行都是一个观测值。机器学习任务的第一步是明确矩阵 X 的确切含义。

在监督学习中，我们还需要一个 目标，一个长度为 N 的向量 y，每个观测值对应一个标量值。这个值是连续的或离散的，具体取决于我们是回归问题还是分类问题。

在 scikit-learn 中，模型通过包含 fit() 和 predict() 方法的类来实现。fit() 方法接受数据矩阵 X 作为输入，对于监督学习模型，还接受 y。此方法用于在给定数据上训练模型。

predict() 方法也接受数据点作为输入（作为 (M,D) 矩阵）。它返回训练模型预测的标签或转换后的点。

普通最小二乘回归

普通最小二乘回归 是最简单的回归方法之一。它通过 X[ij] 的线性组合来逼近输出值 y[i]：

这里，w = (w[1], ..., w[D]) 是（未知的）参数向量。另外， 代表模型的输出。我们希望这个向量与数据点 y 尽可能匹配。当然，精确的相等式 通常是不可能成立的（总会有一些噪声和不确定性——模型始终是对现实的理想化）。因此，我们希望最小化这两个向量之间的差异。普通最小二乘回归方法的核心是最小化以下 损失函数：

这些组件的平方和称为L² 范数。它之所以方便，是因为它导致了可微分的损失函数，从而可以计算梯度，并进行常见的优化过程。

使用线性回归的多项式插值

普通最小二乘回归将线性模型拟合到数据上。该模型在数据点*x[i]和参数w[j]*中都是线性的。在我们的例子中，由于数据点是根据非线性生成模型（一个指数函数）生成的，因此我们获得了较差的拟合。

然而，我们仍然可以使用线性回归方法，模型在*w[j]上是线性的，但在x[i]*上是非线性的。为此，我们需要通过使用多项式函数的基来增加数据集的维度。换句话说，我们考虑以下数据点：

这里，D是最大阶数。因此，输入矩阵X是与原始数据点*x[i]*相关的范德蒙德矩阵。有关范德蒙德矩阵的更多信息，请参见en.wikipedia.org/wiki/Vandermonde_matrix。

在这里，很容易看出，在这些新数据点上训练线性模型等同于在原始数据点上训练多项式模型。

岭回归

使用线性回归的多项式插值如果多项式的阶数过大，可能导致过拟合。通过捕捉随机波动（噪声）而不是数据的一般趋势，模型失去了部分预测能力。这对应于多项式系数*w[j]*的发散。

解决这个问题的方法是防止这些系数无限增大。通过岭回归（也称为Tikhonov 正则化），这是通过向损失函数中添加正则化项来实现的。有关 Tikhonov 正则化的更多信息，请参见en.wikipedia.org/wiki/Tikhonov_regularization。

通过最小化这个损失函数，我们不仅最小化了模型与数据之间的误差（第一项，与偏差相关），还最小化了模型系数的大小（第二项，与方差相关）。偏差-方差权衡通过超参数量化，它指定了损失函数中两项之间的相对权重。

在这里，岭回归导致了具有较小系数的多项式，因此得到了更好的拟合。

交叉验证和网格搜索

岭回归模型相较于普通最小二乘法模型的一个缺点是多了一个额外的超参数 。预测的质量取决于该参数的选择。一种可能的做法是手动微调这个参数，但这个过程可能非常繁琐，并且可能会导致过拟合问题。

为了解决这个问题，我们可以使用 网格搜索；我们遍历多个可能的 值，并评估每个可能值下模型的性能。然后，我们选择使性能最好的参数。

我们如何评估具有给定 值的模型性能？一种常见的解决方案是使用 交叉验证。该过程将数据集分为训练集和测试集。我们在训练集上拟合模型，然后在 测试集 上测试其预测性能。通过在与训练集不同的数据集上测试模型，我们可以减少过拟合的风险。

有许多方法可以将初始数据集分成两部分。一个方法是移除 一个 样本，形成训练集，并将这个样本放入测试集中。这被称为 留一法 交叉验证。对于 N 个样本，我们将得到 N 组训练集和测试集。交叉验证性能是所有这些数据集划分的平均性能。

如我们稍后所见，scikit-learn 实现了几个易于使用的函数，用于进行交叉验证和网格搜索。在这个示例中，存在一个名为RidgeCV的特殊估计器，它实现了特定于岭回归模型的交叉验证和网格搜索过程。使用这个类可以自动为我们找到最佳的超参数 。

还有更多内容…

这里有一些关于最小二乘法的参考资料：

• 维基百科上的普通最小二乘法，详情请见 en.wikipedia.org/wiki/Ordinary_least_squares

• 维基百科上的线性最小二乘法，详情请见 en.wikipedia.org/wiki/Linear_least_squares_(mathematics)

这里有一些关于交叉验证和网格搜索的参考资料：

• scikit-learn 文档中的交叉验证，详情请见 scikit-learn.org/stable/modules/cross_validation.html

• scikit-learn 文档中的网格搜索，详情请见 scikit-learn.org/stable/modules/grid_search.html

• 维基百科上的交叉验证，详情请见 en.wikipedia.org/wiki/Cross-validation_%28statistics%29

这里有一些关于 scikit-learn 的参考资料：

• scikit-learn 基础教程，见 scikit-learn.org/stable/tutorial/basic/tutorial.html

• scikit-learn 在 SciPy 2013 年会上的教程，见 github.com/jakevdp/sklearn_scipy2013

另请参见

• 使用支持向量机进行分类任务 的实例

使用逻辑回归预测谁会在泰坦尼克号上生还

在这个实例中，我们将介绍 逻辑回归，一种基本的分类器。我们还将展示如何使用 网格搜索 和 交叉验证。

我们将在 Kaggle 数据集上应用这些技术，该数据集的目标是根据真实数据预测泰坦尼克号上的生还情况。

小贴士

Kaggle (www.kaggle.com/competition…) 主办机器学习比赛，任何人都可以下载数据集，训练模型，并在网站上测试预测结果。最佳模型的作者甚至可能获得奖品！这是一个很好的入门机器学习的方式。

准备就绪

从本书的 GitHub 仓库下载 Titanic 数据集，链接：github.com/ipython-books/cookbook-data。

数据集来源于 www.kaggle.com/c/titanic-g…。

如何做...

1. 我们导入标准库：

   In [1]: import numpy as np
           import pandas as pd
           import sklearn
           import sklearn.linear_model as lm
           import sklearn.cross_validation as cv
           import sklearn.grid_search as gs
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们使用 pandas 加载训练集和测试集：

   In [2]: train = pd.read_csv('data/titanic_train.csv')
           test = pd.read_csv('data/titanic_test.csv')
   In [3]: train[train.columns[[2,4,5,1]]].head()
   Out[3]:    
      Pclass     Sex  Age  Survived
   0       3    male   22         0
   1       1  female   38         1
   2       3  female   26         1
   3       1  female   35         1
   4       3    male   35         0
   

3. 为了简化示例，我们只保留少数几个字段，并将 sex 字段转换为二进制变量，以便它能被 NumPy 和 scikit-learn 正确处理。最后，我们移除包含 NaN 值的行：

   In [4]: data = train[['Sex', 'Age', 'Pclass', 'Survived']].copy()
           data['Sex'] = data['Sex'] == 'female'
           data = data.dropna()
   

4. 现在，我们将 DataFrame 对象转换为 NumPy 数组，以便将其传递给 scikit-learn：

   In [5]: data_np = data.astype(np.int32).values
           X = data_np[:,:-1]
           y = data_np[:,-1]
   

5. 让我们看看男性和女性乘客的生还情况，按年龄分类：

   In [6]: # We define a few boolean vectors.
           female = X[:,0] == 1
           survived = y == 1
           # This vector contains the age of the passengers.
           age = X[:,1]
           # We compute a few histograms.
           bins_ = np.arange(0, 81, 5)
           S = {'male': np.histogram(age[survived & ~female], 
                                     bins=bins_)[0],
                'female': np.histogram(age[survived & female], 
                                       bins=bins_)[0]}
           D = {'male': np.histogram(age[~survived & ~female], 
                                     bins=bins_)[0],
                'female': np.histogram(age[~survived & 
                                           female], 
                                       bins=bins_)[0]}
   In [7]: # We now plot the data.
           bins = bins_[:-1]
           for i, sex, color in zip((0, 1),
                                    ('male', 'female'),
                                    ('#3345d0', '#cc3dc0')):
               plt.subplot(121 + i)
               plt.bar(bins, S[sex], bottom=D[sex],
                       color=color,
                       width=5, label='survived')
               plt.bar(bins, D[sex], color='k', width=5,
                       label='died')
               plt.xlim(0, 80)
               plt.grid(None)
               plt.title(sex + " survival")
               plt.xlabel("Age (years)")
               plt.legend()
   

   

6. 让我们尝试训练一个 LogisticRegression 分类器，以预测人们基于性别、年龄和舱位的生还情况。我们首先需要创建一个训练集和一个测试集：

   In [8]: # We split X and y into train and test datasets.
           (X_train, X_test, y_train, 
           y_test) = cv.train_test_split(X, y, test_size=.05)
   In [9]: # We instanciate the classifier.
           logreg = lm.LogisticRegression()
   

7. 我们训练模型，并在测试集上获取预测值：

   In [10]: logreg.fit(X_train, y_train)
            y_predicted = logreg.predict(X_test)
   

   下图显示了实际结果和预测结果：

   In [11]: plt.imshow(np.vstack((y_test, y_predicted)),
                       interpolation='none', cmap='bone')
            plt.xticks([]); plt.yticks([])
            plt.title(("Actual and predicted survival "
                       "outcomes on the test set"))
   

   

   在这个截图中，第一行显示了测试集中几个人的生还情况（白色表示生还，黑色表示未生还）。第二行显示了模型的预测值。

8. 为了评估模型的性能，我们使用 cross_val_score() 函数计算交叉验证分数。默认情况下，该函数使用三折分层交叉验证过程，但可以通过 cv 关键字参数进行修改：

   In [12]: cv.cross_val_score(logreg, X, y)
   Out[12]: array([ 0.78661088,  0.78991597,  0.78059072])
   

   这个函数返回每一对训练集和测试集的预测分数（我们在*如何运作…*中提供了更多细节）。

9. LogisticRegression类接受C超参数作为参数。该参数量化了正则化强度。为了找到合适的值，我们可以使用通用的GridSearchCV类进行网格搜索。它接受一个估算器作为输入以及一个包含参数值的字典。这个新的估算器使用交叉验证来选择最佳参数：

   In [13]: grid = gs.GridSearchCV(logreg, 
                               {'C': np.logspace(-5, 5, 50)})
            grid.fit(X_train, y_train)
            grid.best_params_
   Out[13]: {'C':  5.35}
   

10. 这是最佳估算器的性能：

    In [14]: cv.cross_val_score(grid.best_estimator_, X, y)
    Out[14]: array([ 0.78661088,  0.79831933,  0.78481013])
    

    在使用C超参数通过网格搜索选择后，性能略有提高。

它是如何工作的...

逻辑回归不是一个回归模型，而是一个分类模型。然而，它与线性回归有着密切的关系。该模型通过将sigmoid 函数（更准确地说是逻辑函数）应用于变量的线性组合，预测一个二元变量为 1 的概率。sigmoid 函数的方程为：

下图展示了一个逻辑函数：

逻辑函数

如果需要得到一个二元变量，我们可以将值四舍五入为最接近的整数。

参数w是在学习步骤中通过优化过程得到的。

还有更多...

以下是一些参考资料：

• 维基百科上的逻辑回归，网址：en.wikipedia.org/wiki/Logistic_regression

• scikit-learn 文档中的逻辑回归，网址：scikit-learn.org/stable/modules/linear_model.html#logistic-regression

另见

• 开始使用 scikit-learn的教程

• 使用 K 近邻分类器识别手写数字的教程

• 使用支持向量机进行分类任务的教程

使用 K 近邻分类器识别手写数字

在本教程中，我们将展示如何使用K 近邻（K-NN）分类器来识别手写数字。这个分类器是一个简单但强大的模型，适用于复杂的、高度非线性的数据集，如图像。稍后我们将在本教程中详细解释它的工作原理。

操作步骤...

1. 我们导入模块：

   In [1]: import numpy as np
           import sklearn
           import sklearn.datasets as ds
           import sklearn.cross_validation as cv
           import sklearn.neighbors as nb
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 让我们加载digits数据集，这是 scikit-learn 的datasets模块的一部分。这个数据集包含已手动标记的手写数字：

   In [2]: digits = ds.load_digits()
           X = digits.data
           y = digits.target
           print((X.min(), X.max()))
           print(X.shape)
   0.0 16.0
   (1797L, 64L)
   

   在矩阵X中，每一行包含8 * 8=64个像素（灰度值，范围在 0 到 16 之间）。使用行主序排列。

3. 让我们显示一些图像及其标签：

   In [3]: nrows, ncols = 2, 5
           plt.gray()
           for i in range(ncols * nrows):
               ax = plt.subplot(nrows, ncols, i + 1)
               ax.matshow(digits.images[i,...])
               plt.xticks([]); plt.yticks([])
               plt.title(digits.target[i])
   

   

4. 现在，让我们对数据拟合一个 K 近邻分类器：

   In [4]: (X_train, X_test, y_train, 
            y_test) = cv.train_test_split(X, y, test_size=.25)
   In [5]: knc = nb.KNeighborsClassifier()
   In [6]: knc.fit(X_train, y_train);
   

5. 让我们在测试数据集上评估训练后的分类器的得分：

   In [7]: knc.score(X_test, y_test)
   Out[7]: 0.98888888888888893
   

6. 现在，让我们看看我们的分类器能否识别手写数字！

   In [8]: # Let's draw a 1.
           one = np.zeros((8, 8))
           one[1:-1, 4] = 16  # The image values are 
                              # in [0,16].
           one[2, 3] = 16
   In [9]: plt.imshow(one, interpolation='none')
           plt.grid(False)
           plt.xticks(); plt.yticks()
           plt.title("One")
   

   

   我们的模型能识别这个数字吗？让我们来看看：

   In [10]: knc.predict(one.ravel())
   Out[10]: array([1])
   

   干得好！

它是如何工作的...

这个例子展示了如何处理 scikit-learn 中的图像。图像是一个 2D 的*(N, M)矩阵，有NM*个特征。在组合数据矩阵时，这个矩阵需要被展平；每一行都是一个完整的图像。

K 最近邻算法的思想如下：在特征空间中给定一个新点，找到训练集中距离最近的K个点，并将它们中大多数点的标签赋给该新点。

距离通常是欧氏距离，但也可以使用其他距离。

下图展示了在玩具数据集上使用 15 最近邻分类器获得的空间划分（带有三个标签）：

K 最近邻空间划分

数字K是模型的一个超参数。如果K太小，模型泛化能力不好（高方差）。特别是，它对异常值非常敏感。相反，如果K太大，模型的精度会变差。在极端情况下，如果K等于总数据点数，模型将始终预测相同的值，而不考虑输入（高偏差）。有一些启发式方法来选择这个超参数（请参阅下一节）。

需要注意的是，K 最近邻算法并没有学习任何模型；分类器只是存储所有数据点，并将任何新的目标点与它们进行比较。这是基于实例的学习的一个例子。这与其他分类器（如逻辑回归模型）形成对比，后者明确地在训练数据上学习一个简单的数学模型。

K 最近邻方法在复杂的分类问题上表现良好，这些问题具有不规则的决策边界。但是，对于大型训练数据集，它可能计算密集，因为必须为测试计算大量距离。可以使用专用的基于树的数据结构（如 K-D 树或球树）来加速最近邻的搜索。

K 最近邻方法可以用于分类（如本例）以及回归问题。该模型分配最近邻居的目标值的平均值。在这两种情况下，可以使用不同的加权策略。

还有更多内容…

下面是一些参考资料：

• scikit-learn 文档中的 K-NN 算法，网址为scikit-learn.org/stable/modules/neighbors.html

• 维基百科上的 K-NN 算法，网址为en.wikipedia.org/wiki/K-nearest_neighbors_algorithm

• 博客文章介绍如何选择 K 超参数，网址为datasciencelab.wordpress.com/2013/12/27/finding-the-k-in-k-means-clustering/

• 维基百科上的基于实例的学习，网址为en.wikipedia.org/wiki/Instance-based_learning

另请参阅

• 用逻辑回归预测谁将在泰坦尼克号上幸存食谱

• 使用支持向量机进行分类任务食谱

从文本中学习 – 自然语言处理中的朴素贝叶斯

在这个食谱中，我们展示了如何使用 scikit-learn 处理文本数据。处理文本需要仔细的预处理和特征提取。同时，处理高稀疏矩阵也很常见。

我们将学习识别在公共讨论中发布的评论是否被认为是侮辱某个参与者的。我们将使用来自 Impermium 的标注数据集，该数据集是在 Kaggle 竞赛中发布的。

准备工作

从书籍的 GitHub 仓库下载Troll数据集，网址是github.com/ipython-books/cookbook-data。

该数据集可以从 Kaggle 下载，网址是www.kaggle.com/c/detecting…。

如何做到...

1. 让我们导入我们的库：

   In [1]: import numpy as np
           import pandas as pd
           import sklearn
           import sklearn.cross_validation as cv
           import sklearn.grid_search as gs
           import sklearn.feature_extraction.text as text
           import sklearn.naive_bayes as nb
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 让我们用 pandas 打开 CSV 文件：

   In [2]: df = pd.read_csv("data/troll.csv")
   

3. 每一行是一个评论。我们将考虑两个列：评论是否为侮辱性（1）还是非侮辱性（0），以及评论的 Unicode 编码内容：

   In [3]: df[['Insult', 'Comment']].tail()
         Insult                                            Comment
   3942       1  "you are both morons and that is..."
   3943       0  "Many toolbars include spell check...
   3944       0  "@LambeauOrWrigley\xa0\xa0@K.Moss\xa0\n...
   3945       0  "How about Felix? He is sure turning into...
   3946       0  "You're all upset, defending this hipster...
   

4. 现在，我们将定义特征矩阵X和标签y：

   In [4]: y = df['Insult']
   

   从文本中获得特征矩阵并不是一件简单的事。scikit-learn 只能处理数值矩阵。那么，我们如何将文本转换成数值矩阵呢？一种经典的解决方案是，首先提取一个词汇表，即语料库中使用的所有词汇的列表。然后，我们为每个样本计算每个词汇的出现频率。我们最终得到一个稀疏矩阵，这是一个大矩阵，主要包含零。这里，我们只用两行代码来完成。我们将在*它是如何工作的…*部分提供更多细节。

   注意

   这里的一般规则是，当我们的特征是分类变量时（即，某个词的出现，属于固定集合的n种颜色中的一种，等等），我们应该通过考虑每个类别项的二元特征来向量化它。例如，color特征如果是red、green或blue，我们应该考虑三个二元特征：color_red、color_green和color_blue。我们将在*更多内容…*部分提供更多参考资料。

   In [5]: tf = text.TfidfVectorizer()
           X = tf.fit_transform(df['Comment'])
           print(X.shape)
   (3947, 16469)
   

5. 共有 3947 条评论和 16469 个不同的词汇。让我们估算一下这个特征矩阵的稀疏性：

   In [6]: print(("Each sample has ~{0:.2f}% non-zero"
                  "features.").format(
              100 * X.nnz / float(X.shape[0] * X.shape[1])))
   Each sample has ~0.15% non-zero features.
   

6. 现在，我们像往常一样训练分类器。我们首先将数据划分为训练集和测试集：

   In [7]: (X_train, X_test, y_train, 
            y_test) = cv.train_test_split(X, y,
                                          test_size=.2)
   

7. 我们使用伯努利朴素贝叶斯分类器，并对参数进行网格搜索：

   In [8]: bnb = gs.GridSearchCV(nb.BernoulliNB(), 
                                 param_grid={
                           'alpha': np.logspace(-2., 2., 50)})
           bnb.fit(X_train, y_train)
   

8. 让我们检查一下这个分类器在测试数据集上的表现：

   In [9]: bnb.score(X_test, y_test)
   Out[9]: 0.76455696202531642
   

9. 让我们来看一下与最大系数对应的词汇（我们在侮辱性评论中经常找到的词汇）：

   In [10]: # We first get the words corresponding 
            # to each feature.
            names = np.asarray(tf.get_feature_names())
            # Next, we display the 50 words with the largest
            # coefficients.
            print(','.join(names[np.argsort(
                bnb.best_estimator_.coef_[0,:])[::-1][:50]]))
   you,are,your,to,the,and,of,that,is,it,in,like,on,have,for,not,re,just,an,with,so,all,***,***be,get,***,***up,this,what,xa0,don,***,***go,no,do,can,but,***,***or,as,if,***,***who,know,about,because,here,***,***me,was
   

10. 最后，让我们在几个测试句子上测试我们的估算器：

    In [11]: print(bnb.predict(tf.transform([
                 "I totally agree with you.",
                 "You are so stupid.",
                 "I love you."
                 ])))
    [0 1 1]
    

    结果不错，但我们可能能够做得更好。

它是如何工作的...

scikit-learn 实现了多个实用函数，可以从文本数据中获取稀疏特征矩阵。像CountVectorizer()这样的向量化器从语料库中提取词汇（fit），并基于该词汇构建语料库的稀疏表示（transform）。每个样本由词汇的词频表示。训练后的实例还包含属性和方法，用于将特征索引映射到对应的词汇（get_feature_names()）及反向映射（vocabulary_）。

也可以提取 N-gram。这些是连续出现的词对或元组（ngram_range关键字）。

词频可以以不同的方式加权。这里，我们使用了tf-idf，即词频-逆文档频率。这个量度反映了一个词对于语料库的重要性。评论中频繁出现的词语有较高的权重，除非它们出现在大多数评论中（这意味着它们是常见词汇，例如“the”和“and”会被这种技术过滤掉）。

朴素贝叶斯算法是基于特征之间独立性假设的贝叶斯方法。这一强假设极大简化了计算，从而使分类器非常快速且效果不错。

还有更多内容……

以下是一些参考资料：

• scikit-learn 文档中的文本特征提取，见scikit-learn.org/stable/modules/feature_extraction.html#text-feature-extraction

• 维基百科上的词频-逆文档频率，见en.wikipedia.org/wiki/tf-idf

• scikit-learn 文档中的向量化器，见scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.DictVectorizer.html

• 维基百科上的朴素贝叶斯分类器，见en.wikipedia.org/wiki/Naive_Bayes_classifier

• scikit-learn 文档中的朴素贝叶斯，见scikit-learn.org/stable/modules/naive_bayes.html

• Impermium Kaggle 挑战，见blog.kaggle.com/2012/09/26/impermium-andreas-blog/

• scikit-learn 文档中的文档分类示例，见scikit-learn.org/stable/auto_examples/document_classification_20newsgroups.html

小贴士

除了 scikit-learn，它在文本处理方面有很好的支持，我们还应该提到 NLTK（可在www.nltk.org获取），这是一个 Python 中的自然语言工具包。

另见

• 使用逻辑回归预测谁会在泰坦尼克号上生还的教程

• 使用 K 最近邻分类器学习识别手写数字的教程

• 使用支持向量机进行分类任务的教程

使用支持向量机进行分类任务

在本教程中，我们介绍了支持向量机，或SVM。这些强大的模型可以用于分类和回归。在这里，我们展示了如何在一个简单的分类任务中使用线性和非线性 SVM。

如何操作...

1. 让我们导入所需的包：

   In [1]: import numpy as np
           import pandas as pd
           import sklearn
           import sklearn.datasets as ds
           import sklearn.cross_validation as cv
           import sklearn.grid_search as gs
           import sklearn.svm as svm
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们生成二维点，并根据坐标的线性操作分配二进制标签：

   In [2]: X = np.random.randn(200, 2)
           y = X[:, 0] + X[:, 1] > 1
   

3. 现在我们拟合一个线性支持向量分类器（SVC）。该分类器尝试用一个线性边界（这里是线，但更一般地是超平面）将两组点分开：

   In [3]: # We train the classifier.
           est = svm.LinearSVC()
           est.fit(X, y)
   

4. 我们定义了一个函数，用于显示训练好的分类器的边界和决策函数：

   In [4]: # We generate a grid in the square [-3,3 ]².
           xx, yy = np.meshgrid(np.linspace(-3, 3, 500),
                                np.linspace(-3, 3, 500))
           # This function takes a SVM estimator as input.
           def plot_decision_function(est):
               # We evaluate the decision function on the
               # grid.
               Z = est.decision_function(np.c_[xx.ravel(), 
                                               yy.ravel()])
               Z = Z.reshape(xx.shape)
               cmap = plt.cm.Blues
               # We display the decision function on the grid.
               plt.imshow(Z,
                          extent=(xx.min(), xx.max(), 
                                  yy.min(), yy.max()),
                          aspect='auto', origin='lower', 
                          cmap=cmap)
               # We display the boundaries.
               plt.contour(xx, yy, Z, levels=[0],                
                           linewidths=2,
                           colors='k')
               # We display the points with their true labels.
               plt.scatter(X[:, 0], X[:, 1], s=30, c=.5+.5*y, 
                           lw=1, cmap=cmap, vmin=0, vmax=1)
               plt.axhline(0, color='k', ls='--')
               plt.axvline(0, color='k', ls='--')
               plt.xticks(())
               plt.yticks(())
               plt.axis([-3, 3, -3, 3])
   

5. 让我们看看线性 SVC 的分类结果：

   In [5]: plot_decision_function(est)
           plt.title("Linearly separable, linear SVC")
   

   

   线性 SVC 尝试用一条线将点分开，它在这里做得相当不错。

6. 现在，我们使用XOR函数修改标签。如果坐标的符号不同，则该点的标签为 1。这个分类任务是不可线性分割的。因此，线性 SVC 完全失败：

   In [6]: y = np.logical_xor(X[:, 0]>0, X[:, 1]>0)
           # We train the classifier.
           est = gs.GridSearchCV(svm.LinearSVC(), 
                             {'C': np.logspace(-3., 3., 10)})
           est.fit(X, y)
           print("Score: {0:.1f}".format(
                         v.cross_val_score(est, X, y).mean()))
           # Plot the decision function.
           plot_decision_function(est)
           plt.title("XOR, linear SVC")
   Score: 0.6
   

   

7. 幸运的是，可以通过使用非线性核函数来使用非线性 SVC。核函数指定将点进行非线性变换，映射到更高维的空间中。假设该空间中的变换点更容易线性可分。默认情况下，scikit-learn 中的SVC分类器使用径向基函数（RBF）核函数：

   In [7]: y = np.logical_xor(X[:, 0]>0, X[:, 1]>0)
           est = gs.GridSearchCV(svm.SVC(), 
                        {'C': np.logspace(-3., 3., 10),
                         'gamma': np.logspace(-3., 3., 10)})
           est.fit(X, y)
           print("Score: {0:.3f}".format(
                 cv.cross_val_score(est, X, y).mean()))
           plot_decision_function(est.best_estimator_)
           plt.title("XOR, non-linear SVC")
   Score: 0.975
   

   

   这次，非线性 SVC 成功地将这些非线性可分的点进行了分类。

如何操作...

二分类线性 SVC 尝试找到一个超平面（定义为线性方程），以最好地分离两组点（根据它们的标签分组）。另外，还有一个约束条件是，这个分离超平面需要尽可能远离点。这种方法在存在这样的超平面时效果最好。否则，这种方法可能完全失败，就像我们在XOR示例中看到的那样。XOR被认为是一个非线性可分的操作。

scikit-learn 中的 SVM 类有一个C超参数。这个超参数在训练样本的误分类与决策面简洁性之间进行权衡。较低的C值使得决策面平滑，而较高的C值则试图正确分类所有训练样本。这是另一个超参数量化偏差-方差权衡的例子。可以通过交叉验证和网格搜索来选择这个超参数。

线性 SVC 也可以扩展到多分类问题。多分类 SVC 在 scikit-learn 中直接实现。

非线性 SVC 通过考虑从原始空间到更高维空间的非线性转换来工作。这个非线性转换可以增加类别之间的线性可分性。在实践中，所有点积都被核所替代。

非线性 SVC

有几种广泛使用的非线性核。默认情况下，SVC 使用高斯径向基函数：

在这里，是模型的超参数，可以通过网格搜索和交叉验证来选择。

函数不需要显式计算。这就是核技巧；只需要知道核函数*k(x, x')即可。给定核函数k(x, x')*所对应的函数的存在是由函数分析中的数学定理所保证的。

还有更多……

这里有一些关于支持向量机的参考资料：

• 维基百科上的异或，访问地址：en.wikipedia.org/wiki/Exclusive_or

• 维基百科上的支持向量机，访问地址：en.wikipedia.org/wiki/Support_vector_machine

• scikit-learn 文档中的 SVM，访问地址：scikit-learn.org/stable/modules/svm.html

• 维基百科上的核技巧，访问地址：en.wikipedia.org/wiki/Kernel_method

• 关于核技巧的说明，访问地址：www.eric-kim.net/eric-kim-ne…

• 一个关于非线性 SVM 的示例，访问地址：scikit-learn.org/0.11/auto_examples/svm/plot_svm_nonlinear.html（这个示例启发了这个配方）

另见：

• 使用逻辑回归预测谁会在泰坦尼克号上生存的配方

• 使用 K 近邻分类器识别手写数字的配方

使用随机森林选择回归的重要特征

决策树经常被用来表示工作流程或算法。它们还构成了一种非参数监督学习方法。在训练集上学习一个将观察值映射到目标值的树，并且可以对新观察值做出预测。

随机森林是决策树的集成。通过训练多个决策树并将其结果聚合，形成的模型比任何单一的决策树更具表现力。这个概念就是集成学习的目的。

有许多类型的集成方法。随机森林是自助聚合（也称为bagging）的一种实现，其中模型是在从训练集中随机抽取的子集上训练的。

随机森林能提供每个特征在分类或回归任务中的重要性。在本实例中，我们将使用一个经典数据集，利用该数据集中的各类房屋邻里指标来找出波士顿房价的最具影响力特征。

如何实现...

1. 我们导入所需的包：

   In [1]: import numpy as np
           import sklearn as sk
           import sklearn.datasets as skd
           import sklearn.ensemble as ske
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 我们加载波士顿数据集：

   In [2]: data = skd.load_boston()
   

   数据集的详细信息可以在 data['DESCR'] 中找到。以下是一些特征的描述：

   • CRIM：各镇的人均犯罪率

   • NOX：氮氧化物浓度（每千万分之一）

   • RM：每套住房的平均房间数

   • AGE：1940 年前建造的自有住房单位的比例

   • DIS：到波士顿五个就业中心的加权距离

   • PTRATIO：各镇的师生比

   • LSTAT：低社会经济地位人口的比例

   • MEDV：以 $1000 为单位的自有住房中位数价格

   目标值是 MEDV。

3. 我们创建一个 RandomForestRegressor 模型：

   In [3]: reg = ske.RandomForestRegressor()
   

4. 我们从这个数据集中获取样本和目标值：

   In [4]: X = data['data']
           y = data['target']
   

5. 让我们来拟合模型：

   In [5]: reg.fit(X, y)
   

6. 特征的重要性可以通过 reg.feature_importances_ 查看。我们按重要性递减顺序对其进行排序：

   In [6]: fet_ind = np.argsort(reg.feature_importances_) \
                                                       [::-1]
           fet_imp = reg.feature_importances_[fet_ind]
   

7. 最后，我们绘制特征重要性的直方图：

   In [7]: ax = plt.subplot(111)
           plt.bar(np.arange(len(fet_imp)), 
                   fet_imp, width=1, lw=2)
           plt.grid(False)
           ax.set_xticks(np.arange(len(fet_imp))+.5)
           ax.set_xticklabels(data['feature_names'][fet_ind])
           plt.xlim(0, len(fet_imp))
   

   

   我们发现 LSTAT（低社会经济地位人口比例）和 RM（每套住房的房间数）是决定房价的最重要特征。作为说明，这里是房价与 LSTAT 的散点图：

   In [8]: plt.scatter(X[:,-1], y)
           plt.xlabel('LSTAT indicator')
           plt.ylabel('Value of houses (k$)')
   

   

工作原理...

训练决策树可以使用多种算法。scikit-learn 使用 CART（分类与回归树）算法。该算法利用在每个节点上产生最大信息增益的特征和阈值构建二叉树。终端节点给出输入值的预测结果。

决策树易于理解。它们还可以通过 pydot 这一 Python 包进行可视化，用于绘制图形和树形结构。当我们希望了解决策树究竟学到了什么时，这非常有用（白盒模型）；每个节点上的观测条件可以通过布尔逻辑简洁地表达。

然而，决策树可能会受到过拟合的影响，特别是在树过深时，且可能不稳定。此外，尤其是在使用贪婪算法训练时，全局收敛到最优模型并没有保证。这些问题可以通过使用决策树集成方法来缓解，特别是随机森林。

在随机森林中，多个决策树在训练数据集的自助样本上进行训练（通过有放回的随机抽样）。预测结果通过各个树的预测结果平均值来得出（自助聚合或袋装法）。此外，在每个节点上会随机选择特征的子集（随机子空间方法）。这些方法能使得整体模型优于单棵树。

还有更多...

以下是一些参考资料：

• scikit-learn 文档中的集成学习，链接：scikit-learn.org/stable/modules/ensemble.html

• RandomForestRegressor的 API 参考，链接：scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html

• 维基百科上的随机森林，链接：en.wikipedia.org/wiki/Random_forest

• 维基百科上的决策树学习，链接：en.wikipedia.org/wiki/Decision_tree_learning

• 维基百科上的自助法聚合，链接：en.wikipedia.org/wiki/Bootstrap_aggregating

• 维基百科上的随机子空间方法，链接：en.wikipedia.org/wiki/Random_subspace_method

• 维基百科上的集成学习，链接：en.wikipedia.org/wiki/Ensemble_learning

另见

• 使用支持向量机进行分类任务的食谱

通过主成分分析（PCA）来减少数据集的维度

在之前的食谱中，我们介绍了有监督学习方法；我们的数据点带有离散或连续的标签，算法能够学习从点到标签的映射关系。

从这个食谱开始，我们将介绍无监督学习方法。这些方法在运行有监督学习算法之前可能会有所帮助，能够为数据提供初步的洞察。

假设我们的数据由没有标签的点*x[i]*组成。目标是发现这些点集合中的某种隐藏结构。数据点通常具有内在的低维性：少量特征就足以准确描述数据。然而，这些特征可能被许多与问题无关的特征所掩盖。降维可以帮助我们找到这些结构。这些知识可以显著提升后续有监督学习算法的性能。

无监督学习的另一个有用应用是数据可视化；高维数据集很难在二维或三维中进行可视化。将数据点投影到子空间或子流形上，可以得到更有趣的可视化效果。

在本食谱中，我们将演示一种基本的无监督线性方法——主成分分析（PCA）。该算法可以让我们将数据点线性地投影到低维子空间上。在主成分方向上，这些向量构成了该低维子空间的基，数据点的方差最大。

我们将使用经典的鸢尾花数据集作为示例。这个数据集包含了 150 朵鸢尾花的花瓣和萼片的宽度和长度。这些花属于三种类别之一：Iris setosa、Iris virginica和Iris versicolor。我们可以在这个数据集中获取到类别信息（有标签数据）。然而，由于我们感兴趣的是展示一种无监督学习方法，我们将只使用不带标签的数据矩阵。

如何操作...

1. 我们导入 NumPy、matplotlib 和 scikit-learn：

   In [1]: import numpy as np
           import sklearn
           import sklearn.decomposition as dec
           import sklearn.datasets as ds
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 鸢尾花数据集可以在 scikit-learn 的datasets模块中找到：

   In [2]: iris = ds.load_iris()
           X = iris.data
           y = iris.target
           print(X.shape)
   (150L, 4L)
   

3. 每一行包含与花朵形态学相关的四个参数。让我们展示前两个维度。颜色反映了花朵的鸢尾花种类（标签，介于 0 和 2 之间）：

   In [3]: plt.scatter(X[:,0], X[:,1], c=y,
                       s=30, cmap=plt.cm.rainbow)
   

   

   提示

   如果你正在阅读这本书的印刷版，可能无法区分颜色。你可以在本书的网站上找到彩色图像。

4. 我们现在对数据集应用 PCA 以获得变换后的矩阵。这个操作可以通过 scikit-learn 中的一行代码完成：我们实例化一个PCA模型并调用fit_transform()方法。这个函数计算主成分并将数据投影到这些主成分上：

   In [4]: X_bis = dec.PCA().fit_transform(X)
   

5. 我们现在展示相同的数据集，但采用新的坐标系统（或者等效地，是初始数据集的线性变换版本）：

   In [5]: plt.scatter(X_bis[:,0], X_bis[:,1], c=y,
                       s=30, cmap=plt.cm.rainbow)
   

   

   属于同一类别的点现在被分到了一起，即使PCA估计器没有使用标签。PCA 能够找到一个最大化方差的投影，这里对应的是一个类别分得很好的投影。

6. scikit.decomposition模块包含了经典PCA估计器的几个变体：ProbabilisticPCA、SparsePCA、RandomizedPCA、KernelPCA等。作为示例，让我们来看一下KernelPCA，PCA 的非线性版本：

   In [6]: X_ter = dec.KernelPCA(kernel='rbf'). \
                                          fit_transform(X)
           plt.scatter(X_ter[:,0], X_ter[:,1], c=y, s=30,
                       cmap=plt.cm.rainbow)
   

   

它是如何工作的...

让我们来看看 PCA 背后的数学思想。这种方法基于一种叫做奇异值分解（SVD）的矩阵分解：

这里，X是*(N,D)的数据矩阵，U和V是正交矩阵，而是一个(N,D)*的对角矩阵。

PCA 将X转换为由以下定义的X'：

的对角线元素是X的奇异值。根据惯例，它们通常按降序排列。U的列是称为X的左奇异向量的正交归一化向量。因此，*X'*的列是左奇异向量乘以奇异值。

最终，PCA 将初始的一组观测值（这些观测值可能包含相关变量）转换为一组线性无关的变量，称为主成分。

第一个新特性（或第一个主成分）是将所有原始特征进行转换，使得数据点的离散度（方差）在该方向上最大。在随后的主成分中，方差逐渐减小。换句话说，PCA 为我们提供了一种数据的新表示，其中新特征按照它们对数据点变异性的解释程度进行排序。

还有更多...

这里有一些进一步的参考：

• 维基百科上的鸢尾花数据集，访问 en.wikipedia.org/wiki/Iris_flower_data_set

• 维基百科上的 PCA，访问 en.wikipedia.org/wiki/Principal_component_analysis

• 维基百科上的 SVD 分解，访问 en.wikipedia.org/wiki/Singular_value_decomposition

• 鸢尾花数据集示例，访问 scikit-learn.org/stable/auto_examples/datasets/plot_iris_dataset.html

• scikit-learn 文档中的分解方法，访问 scikit-learn.org/stable/modules/decomposition.html

• 使用 scikit-learn 的无监督学习教程，访问 scikit-learn.org/dev/tutorial/statistical_inference/unsupervised_learning.html

另见

• 使用聚类检测数据集中的隐藏结构的配方

使用聚类检测数据集中的隐藏结构

无监督学习的一个重要部分是聚类问题。其目标是以完全无监督的方式将相似的点聚集在一起。聚类是一个困难的问题，因为簇（或组）的定义并不总是明确的。在大多数数据集中，声称两个点应该属于同一簇可能是依赖于上下文，甚至可能是主观的。

聚类算法有很多种。我们将在这个配方中看到其中的几个，并将它们应用于一个玩具示例。

如何操作...

1. 让我们导入库：

   In [1]: from itertools import permutations
           import numpy as np
           import sklearn
           import sklearn.decomposition as dec
           import sklearn.cluster as clu
           import sklearn.datasets as ds
           import sklearn.grid_search as gs
           import matplotlib.pyplot as plt
           %matplotlib inline
   

2. 让我们生成一个包含三个簇的随机数据集：

   In [2]: X, y = ds.make_blobs(n_samples=200, n_features=2,  
                                centers=3)
   

3. 我们需要一些函数来重新标记并展示聚类算法的结果：

   In [3]: def relabel(cl):
               """Relabel a clustering with three clusters
               to match the original classes."""
               if np.max(cl) != 2:
                   return cl
               perms = np.array(list(permutations((0, 1, 2))))
               i = np.argmin([np.sum(np.abs(perm[cl] - y))
                              for perm in perms])
               p = perms[i]
               return p[cl]
   In [4]: def display_clustering(labels, title):
               """Plot the data points with the cluster   
               colors."""
               # We relabel the classes when there are 3 
               # clusters.
               labels = relabel(labels)
               # Display the points with the true labels on 
               # the left, and with the clustering labels on
               # the right.
               for i, (c, title) in enumerate(zip(
                       [y, labels], ["True labels", title])):
                   plt.subplot(121 + i)
                   plt.scatter(X[:,0], X[:,1], c=c, s=30, 
                               linewidths=0,
                               cmap=plt.cm.rainbow)
                   plt.xticks([]); plt.yticks([])
                   plt.title(title)
   

4. 现在，我们使用K-means算法对数据集进行聚类，这是一种经典且简单的聚类算法：

   In [5]: km = clu.KMeans()
           km.fit(X)
           display_clustering(km.labels_, "KMeans")
   

   

   小贴士

   如果你正在阅读本书的印刷版，可能无法区分颜色。你可以在本书网站上找到彩色图片。

5. 这个算法在初始化时需要知道簇的数量。然而，一般来说，我们并不一定知道数据集中簇的数量。在这里，让我们尝试使用n_clusters=3（这算是作弊，因为我们恰好知道有 3 个簇！）：

   In [6]: km = clu.KMeans(n_clusters=3)
           km.fit(X)
           display_clustering(km.labels_, "KMeans(3)")
   

   

6. 让我们尝试一些在 scikit-learn 中实现的其他聚类算法。由于 API 非常简洁，尝试不同的方法变得非常容易；只需要改变类名即可：

   In [7]: plt.subplot(231)
           plt.scatter(X[:,0], X[:,1], c=y, s=30,
                       linewidths=0, cmap=plt.cm.rainbow)
           plt.xticks([]); plt.yticks([])
           plt.title("True labels")
           for i, est in enumerate([clu.SpectralClustering(3),
                                    clu.AgglomerativeClustering(3),
                                    clu.MeanShift(),
                                    clu.AffinityPropagation(),
                                    clu.DBSCAN()]):
               est.fit(X)
               c = relabel(est.labels_)
               plt.subplot(232 + i)
               plt.scatter(X[:,0], X[:,1], c=c, s=30,
                           linewidths=0, cmap=plt.cm.rainbow)
               plt.xticks([]); plt.yticks([])
               plt.title(est.__class__.__name__)
   

   

   前两个算法需要输入簇的数量。接下来的两个算法不需要，但它们能够找到正确的簇数：3。最后一个算法没有找到正确的簇数（这就是过度聚类——发现了过多的簇）。

它是如何工作的...

K-means 聚类算法通过将数据点 x[j] 分配到 K 个簇 S[i] 中，从而最小化簇内平方和：

在这里， 是簇 i 的中心（即 S[i] 中所有点的平均值）。

尽管精确地解决这个问题非常困难，但存在一些近似算法。一个流行的算法是 Lloyd's 算法。它的基本过程是从一组初始的 K 均值出发 ，并交替进行两个步骤：

• 在分配步骤中，数据点会被分配到与最近均值相关的簇

• 在更新步骤中，从最后的分配中重新计算均值

该算法收敛到一个解，但并不能保证是最优解。

期望最大化算法可以看作是 K-means 算法的概率版本。它在 scikit-learn 的 mixture 模块中实现。

本方案中使用的其他聚类算法在 scikit-learn 文档中有说明。没有哪种聚类算法在所有情况下都能表现得比其他算法更好，每种算法都有其优缺点。你可以在下一节的参考资料中找到更多细节。

还有更多内容...

以下是一些参考文献：

• 维基百科上的 K-means 聚类算法，链接为 en.wikipedia.org/wiki/K-means_clustering

• 维基百科上的期望最大化算法，链接为 en.wikipedia.org/wiki/Expectation-maximization_algorithm

• scikit-learn 文档中的聚类内容，可在 scikit-learn.org/stable/modules/clustering.html 查看

另请参见

• 使用主成分分析（PCA）降维数据集的方案