NumPy-秘籍中文第二版-全-

NumPy 秘籍中文第二版（全）

原文：NumPy Cookbook - Second Edition

协议：CC BY-NC-SA 4.0

零、前言

第二版增加了有关新的 NumPy 功能和数据分析的两个新章节。我们的 NumPy 用户生活在令人兴奋的时代。与 NumPy 相关的新发展似乎每周，甚至每天都会引起我们的注意。在第一版的时候，创建了 NumFocus，它是 NumPy 可用科学开放代码基金会的缩写。同时宣布了 Numba 项目-使用 LLVM 的 NumPy 感知型动态 Python 编译器。此外，谷歌为他们的云产品 Google App Engine 添加了支持。

将来，我们可以期望改进对 GPU 和 CPU 集群的并发支持。使用 NumPy 数组可以进行类似 OLAP 的查询。这是个好消息，但我们必须不断提醒自己，NumPy 在科学（Python）软件生态系统中并不孤单。有 SciPy，matplotlib（一个非常有用的 Python 绘图库），IPython（一个交互式外壳）和 Scikits。在 Python 生态系统之外，R，C 和 Fortran 等语言非常流行。我们将介绍与这些环境交换数据的细节。

本书涵盖的内容

第 1 章，“使用 IPython”介绍了 IPython，这是一个以其 shell 最著名的工具包。基于 Web 的笔记本是一项激动人心的功能，此处详细介绍。考虑一下 MATLAB 和 Mathematica，但是在您的浏览器中，它是开源的并且免费的。

第 2 章，“高级索引和数组概念”显示，由于强大的索引机制，NumPy 具有非常易于使用的数组。本章介绍一些更高级和更棘手的索引技术。

第 3 章，“掌握常用功能”试图记录每个 NumPy 用户应该知道的最基本的功能。 NumPy 具有许多功能，在本书中甚至无法提及！

第 4 章，“将 NumPy 与世界其他地区联系起来”，现实世界中遇到的编程语言，库和工具的数量令人难以置信。一些软件在云上运行，而另一些则驻留在本地计算机或远程服务器上。能够在这样的环境中适应 NumPy 并将其连接起来与编写独立的 NumPy 代码一样重要。

第 5 章，“音频和图像处理”假定，当您想到 NumPy 时，您可能不会想到声音或图像。阅读本章后，情况将会改变。

第 6 章，“特殊数组和通用函数”介绍了一些非常技术性的主题。本章介绍如何执行字符串操作，忽略非法值以及存储异构数据。

第 7 章，“分析和调试”展示了生产优质软件所需的技能。我们演示了几种方便的性能分析和调试工具。

第 8 章，“质量保证”值得关注，因为它与质量有关。我们使用 NumPy 测试工具讨论了常见的方法和技术，例如单元测试，模拟和 BDD。

第 9 章，“使用 Cython 加速代码”介绍了 Cython，它试图结合 C 的速度和 Python 的优势。我们从 NumPy 的角度向您展示 Cython 的工作原理。

第 10 章，“Scikits 的乐趣”涵盖了 Scikits，这是迷人的科学 Python 生态系统的又一组成部分。快速导览将指导您完成一些最有用的 Scikits 项目。

第 11 章，“最新和最强的 NumPy”展示了第一版中未涵盖的新功能。

第 12 章，“使用 NumPy 进行探索性和预测性数据分析”，介绍了现实的气象数据分析。我已经在第二版中添加了本章。

这本书需要什么

要尝试本书中的代码示例，您将需要最新的 NumPy 版本。这意味着您还需要拥有 NumPy 支持的 Python 版本之一。本书提供了用于安装其他相关包的秘籍。

这本书适合谁？

本书适用于具有 Python 和 NumPy 基本知识的科学家，工程师，程序员或分析人员，他们想要进入下一个层次。此外，还需要对数学和统计学有一定的亲和力，或者至少是有兴趣的。

约定

在本书中，您会发现许多可以区分不同类型信息的文本样式。以下是这些样式的一些示例，并解释了其含义。

文本，数据库表名称，文件夹名称，文件名，文件扩展名，路径名，虚拟 URL，用户输入和 Twitter 句柄中的代码字如下所示：“我们可以通过使用include指令包括其他上下文。”

代码块设置如下：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
import matplotlib.pyplot as plt

def get_indices(high, size):
   #2\. Generate random indices
   return np.random.randint(0, high, size)

当我们希望引起您对代码块特定部分的注意时，相关的行或项目将以粗体显示：

from sklearn.datasets import load_sample_images
import matplotlib.pyplot as plt
import skimage.feature

dataset = load_sample_images()
img = dataset.images[0] 
edges = skimage.feature.canny(img[..., 0])
plt.axis('off')
plt.imshow(edges)
plt.show()

任何命令行输入或输出的编写方式如下：

$ sudo easy_install patsy

新术语和 重要词用粗体显示。您在屏幕上看到的字词，例如在菜单或对话框中，将以如下形式出现：“打印按钮实际上并未打印笔记本。”

注意

警告或重要提示会出现在这样的框中。

提示

提示和技巧如下所示。

一、使用 IPython

在本章中，我们将介绍以下秘籍：

安装 IPython
使用 IPython 作为 Shell
阅读手册页
安装 matplotlib
运行 IPython 笔记本
导出 IPython 笔记本
导入网络笔记本
配置笔记本服务器
探索 SymPy 配置文件

简介

IPython，可从 ipython.org 获得，是一个免费的开源项目，可用于 Linux，Unix，MacOSX，和 Windows。 IPython 作者仅要求您在使用 IPython 的任何科学著作中引用 IPython。 IPython 提供了用于交互式计算的架构。该项目最值得注意的部分是 IPython shell。 IPython 提供了以下组件，其中包括：

交互式 Python Shell（基于终端的 Qt 应用）
一个 Web 笔记本（在 IPython 0.12 和更高版本中可用），支持富媒体和绘图

IPython 与 Python 2.5、2.6、2.7、3.1、3.2、3.3 和 3.4 版本兼容。兼容性取决于 IPython 版本。例如，IPython 2.3.0 需要 Python 2.7 或 3.3+。

您可以在这个页面上在云中尝试 IPython 而不将其安装在系统上。与本地安装的软件相比，会有一些延迟，因此这不如真实的软件好。但是， IPython 交互式外壳程序中可用的大多数功能似乎都可用。 PythonAnywhere 也有一个 Vi(m) 编辑器，如果您喜欢 vi 的话，那显然很棒。您可以从 IPython 会话中保存和编辑文件。

安装 IPython

可以通过多种方式安装 IPython ，具体取决于您的操作系统。对于基于终端的外壳，需要依赖readline。网络笔记本需要tornado和zmq。

除了安装 IPython 外，我们还将安装setuptools，它为您提供easy_install命令。 easy_install命令是 Python 的流行package管理器。一旦拥有easy_install即可安装pip。 pip命令类似于easy_install，并添加了诸如卸载之类的选项。

操作步骤

本节介绍如何在 Windows，MacOSX 和 Linux 上安装 IPython。它还描述了如何使用easy_install和pip或从源代码安装 IPython 及其依赖项：

在 Windows 上安装 IPython 和 Setuptools：IPython 网站上提供了适用于 Python2 或 Python3 的二进制 Windows 安装程序。另见这里。

使用安装程序安装 Setuptools。然后安装pip，如下所示：
```
cd C:\Python27\scripts
python .\easy_install-27-script.py pip
```
在 MacOSX 上安装 IPython：如果需要，请安装 Apple Developer Tools（Xcode）。可以在这个页面上找到 Xcode。请遵循easy_install/pip说明或本节后面提供的从源安装的说明。
在 Linux 上安装 IPython：由于 Linux 发行版太多，因此本节将不详尽：
- 在 Debian 上，键入以下命令：
```
$ su – aptitude install ipython python-setuptools
```
- 在 Fedora 上，魔术命令如下：
```
$ su – yum install ipython python-setuptools-devel
```
- 以下命令将在 Gentoo 上安装 IPython：
```
$ su – emerge ipython
```
- 对于 Ubuntu，安装命令如下：
```
$ sudo apt-get install ipython python-setuptools
```
使用easy_install或pip安装 IPython：使用以下命令，通过easy_install安装 IPython 和本章中的秘籍所需的所有依赖项：
```
$ sudo easy_install ipython pyzmq tornado readline
```
或者，您可以通过在终端中键入以下命令，首先使用easy_install安装pip：
```
$ sudo easy_install pip
```
之后，使用pip安装 IPython：
```
$ sudo pip install ipython pyzmq tornado readline
```
从源代码安装：如果您想使用最新的开发版本，则从源代码安装适合您：
1. 从这个页面下载最新的源存档。
2. 从存档中解压缩源代码：
```
$ tar xzf ipython-<version>.tar.gz
```
3. 相反，如果您已安装 Git，则可以克隆 Git 存储库：
```
$ git clone https://github.com/ipython/ipython.git
```
4. 转到下载源中的根目录：
```
$ cd ipython
```
5. 运行安装脚本。这可能需要您使用sudo运行命令，如下所示：
```
$ sudo python setup.py install
```

工作原理

我们使用几种方法安装了 IPython。这些方法大多数都安装最新的稳定版本，但从源代码安装时除外，这将安装开发版本。

另见

来自 IPython 官方网站的说明

使用 IPython 作为 Shell

科学家和工程师习惯进行实验。 IPython 是科学家根据实验而创建的。 IPython 提供的交互式环境被许多人视为对 MATLAB，Mathematica，Maple 和 R 的直接回答。

以下是 IPython Shell 的功能列表：

制表符补全
历史机制
内联编辑
使用%run调用外部 Python 脚本的功能
调用与操作系统外壳程序交互的魔术函数的能力
访问系统命令
pylab开关
访问 Python 调试器和分析器

操作步骤

本节描述了如何使用 IPython Shell：

pylab：pylab开关会自动导入所有 SciPy，NumPy 和 matplotlib 包。没有此开关，我们将不得不自行导入这些包。

我们需要做的就是在命令行中输入以下指令：

$ ipython --pylab
Type "copyright", "credits" or "license" for more information.

IPython 2.4.1 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

Welcome to pylab, a matplotlib-based Python environment [backend: MacOSX].
For more information, type 'help(pylab)'.

In [1]: quit()
quit() or Ctrl + D quits the IPython shell.

保存会话：我们可能希望能够返回到我们的实验。在 IPython 中，很容易保存会话以供以后使用。这是通过以下命令完成的：

In [1]: %logstart
Activating auto-logging. Current session state plus future input saved.
Filename       : ipython_log.py
Mode           : rotate
Output logging : False
Raw input log  : False
Timestamping   : False
State          : active

可以使用以下命令关闭日志记录：

In [9]: %logoff
Switching logging OFF

执行系统 Shell 命令：您可以在默认的 IPython 配置文件中通过在命令前面加上!符号来执行系统外壳命令。例如，以下输入将获取当前日期：
```
In [1]: !date
```
实际上，任何以!为前缀的行都将发送到系统外壳。我们还可以存储命令输出，如下所示：
```
In [2]: thedate = !date
In [3]: thedate
```
显示历史记录：我们可以使用%hist命令显示命令的历史记录，如下所示：
```
In [1]: a = 2 + 2

In [2]: a
Out[2]: 4

In [3]: %hist
a = 2 + 2
a
%hist
```
这是命令行界面（CLI）环境中的常见功能。我们还可以使用-g开关查询历史记录：
```
In [5]: %hist -g a = 2
 1: a = 2 + 2
```

工作原理

我们看到了许多所谓的魔术功能在起作用。这些功能以%字符开头。如果单独在行中使用魔术函数，则%前缀是可选的。

另见

IPython 作为系统外壳，来自 IPython 官方网站。

阅读手册页

我们可以使用help命令打开有关 NumPy 函数的文档。不必知道函数的名称。我们可以输入几个字符，然后让制表符完成工作。例如，让我们浏览arange()函数的可用信息。

操作步骤

我们可以通过以下两种方式之一浏览可用信息：

调用帮助功能：调用help命令。输入该功能的几个字符，然后按Tab键（请参见以下屏幕截图）：
带问号的查询：另一个选择是在函数名称后添加问号。然后，您当然需要知道函数名称，但是不必键入help命令：
```
In [3]: arange?
```

工作原理

制表符的完成取决于readline，因此您需要确保已安装它。问号为您提供docstrings中的信息。

安装 matplotlib

matplotlib（按惯例所有小写字母）是一个非常有用的 Python 绘图库，对于以下秘籍以及以后的内容，我们将需要它。它取决于 NumPy，但很可能已经安装了 NumPy。

操作步骤

我们将看到如何在 Windows，Linux 和 MacOSX 上安装 matplotlib，以及如何从源代码安装它：

在 Windows 上安装 matplotlib：您可以使用 Enthought 发行版（也称为 Canopy）。

可能需要将msvcp71.dll文件放在您的C:\Windows\system32目录中。您可以从这里获取。
在 Linux 上安装 matplotlib：让我们看看如何在 Linux 的各种发行版中安装 matplotlib：

这是 Debian 和 Ubuntu 上的安装命令：
```
$ sudo apt-get install python-matplotlib
```
- 在 Fedora/Redhat 上的安装命令如下：
```
$ su - yum install python-matplotlib
```
从源代码安装：您可以下载 Sourceforge 的tar.gz版本或从 Git 存储库下载最新的源代码：
```
$ git clone git://github.com/matplotlib/matplotlib.git
```
下载后，请使用以下命令照常构建和安装 matplotlib：
```
$ cd matplotlib
$ sudo python setup.py install
```
在 MacOSX 上安装 matplotlib：从这里获取最新的 DMG 文件。您还可以使用 Mac Ports，Fink 或 Homebrew 包管理器。

另见

matplotlib 官方文档中的说明
在这里中说明了如何在 SciPy 栈中安装。

运行 IPython 笔记本

IPython 具有一项令人兴奋的功能-网络笔记本。所谓的笔记本服务器可以通过 Web 服务笔记本电脑。现在，我们可以启动笔记本服务器并获得基于 Web 的 IPython 环境。该环境具有常规 IPython 环境具有的大多数功能。 IPython 笔记本的功能包括：

显示图像和内嵌图
在文本单元格中使用 HTML 和 Markdown（这是一种简化的类似 HTML 的语言，请参见这里）
笔记本的导入导出

准备

在开始之前，我们应确保已安装所有必需的软件。依赖于tornado和zmq。有关更多信息，请参见本章中的“安装 IPython”秘籍。

操作步骤

运行笔记本：我们可以使用以下命令启动笔记本：
```
$ ipython notebook

[NotebookApp] Using existing profile dir: u'/Users/ivanidris/.ipython/profile_default'
[NotebookApp] The IPython Notebook is running at: http://127.0.0.1:8888
[NotebookApp] Use Control-C to stop this server and shut down all kernels.
```
如您所见，我们正在使用默认配置文件。服务器在本地计算机上的端口 8888 上启动。稍后，您将在本章中学习如何配置这些设置。笔记本在默认浏览器中打开；这也是可配置的（请参见以下屏幕截图）：

IPython 在启动笔记本的目录中列出了所有笔记本。在本示例中，未找到笔记本。可以通过按Ctrl + C停止服务器。
以pylab模式运行笔记本：使用以下命令以pylab模式运行网络笔记本：
```
$ ipython notebook --pylab
```
这将加载SciPy，NumPy和matplotlib模块。
运行带有嵌入式图形的笔记本：我们可以使用以下命令，通过inline指令显示嵌入式 matplotlib 图：
```
$ ipython notebook --pylab inline
```
以下步骤演示了 IPython 笔记本的功能：
1. 单击新笔记本按钮以创建新笔记本。
2. 使用arange()函数创建一个数组。输入以下屏幕快照中所示的命令，然后单击单元格 / 运行：
3. 接下来输入以下命令，然后按Enter。您将在Out [2]中看到输出，如以下屏幕截图所示：
4. 将 sinc()函数应用于数组并绘制结果，如以下屏幕快照所示：

工作原理

内联选项使可以显示内联 matplotlib 图。与pylab模式结合使用时，无需导入 NumPy，SciPy 和 matplotlib 包。

另见

本章中的“安装 IPython”秘籍
笔记本的示例
sinc()函数的文档
plot()函数的文档

导出 IPython 笔记本

有时，您想与朋友或同事交换笔记本。网络笔记本提供了几种导出数据的方法。

操作步骤

可以使用以下选项导出 Web 笔记本：

打印选项：打印按钮实际上并未打印笔记本，但允许您将笔记本导出为 PDF 或 HTML 文档。

下载笔记本：使用“下载”按钮将笔记本下载到您选择的位置。我们可以指定将笔记本下载为.py文件（只是普通的 Python 程序），还是下载为 JSON 格式的.ipynb文件。导出后，我们在前一个秘籍中创建的笔记本如下所示：

{
 "metadata": {
  "name": "Untitled1"
 }, 
 "nbformat": 2, 
 "worksheets": [
  {
    "cells": [
    {
      "cell_type": "code", 
      "collapsed": false, 
      "input": [
        "plot(sinc(a))"
      ], 
      "language": "python", 
      "outputs": [
      {
        "output_type": "pyout", 
        "prompt_number": 3, 
        "text": [
          "[<matplotlib.lines.Line2D at 0x103d9c690>]"
        ]
      }, 
      {
        "output_type": "display_data", 
        "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD9CAYAAABZVQdHAAAABHNCSVQICAgIf...
          mgkAAAAASUVORK5CYII=\n"
      }
      ], 
      "prompt_number": 3
    }
    ]
  }
  ]
}

注意

为简洁起见，一些文本已被省略。该文件不适合编辑甚至阅读，但如果忽略图像表示部分，则可读性很强。有关 JSON 的更多信息，请参见这里。

保存笔记本：使用保存按钮保存笔记本。这会自动以本地 JSON 格式.ipynb导出笔记本。该文件将存储在最初启动 IPython 的目录中。

导入网络笔记本

可以将 Python 脚本作为 Web 笔记本导入。显然，我们也可以导入以前导出的笔记本。

操作步骤

此秘籍向您展示如何将 Python 脚本作为 Web 笔记本导入。

使用以下命令加载 Python 脚本：

% load vectorsum.py

以下屏幕截图显示了从“NumPy 入门指南”加载vectorsum.py后进入笔记本页面的示例：

配置笔记本服务器

公共笔记本服务器需要安全。您应该设置一个密码并使用 SSL 证书来连接它。我们需要证书来通过 HTTPS 提供安全的通信（有关更多信息，请参见这里）。 HTTPS 在互联网上广泛使用的标准 HTTP 协议之上添加了一个安全层。 HTTPS 还对从客户端发送到服务器并返回的数据进行加密。证书颁发机构通常是为网站颁发证书的商业组织。 Web 浏览器具有证书颁发机构的知识，并且可以识别证书。网站管理员需要创建证书并由证书颁发机构签名。

操作步骤

以下步骤描述了如何配置安全的笔记本服务器：

我们可以从 IPython 生成密码。启动一个新的 IPython 会话并输入以下命令：
```
In [1]: from IPython.lib import passwd

In [2]: passwd()
Enter password: 
Verify password:
Out[2]: 'sha1:0e422dfccef2:84cfbcbb3ef95872fb8e23be3999c123f862d856'
```
在第二个输入行，将提示您输入密码。您需要记住该密码。生成一个长字符串。复制此字符串，因为稍后将需要它。

要创建 SSL 证书，您需要在路径中使用openssl命令。

设置openssl命令不是火箭科学，但这可能很棘手。不幸的是，这超出了本书的范围。好的一面是，在线上有很多教程可以帮助您进一步发展。

执行以下命令，以mycert.pem为名称创建证书：

$ openssl req -x509 -nodes -days 365 -newkey rsa:1024 -keyout mycert.pem -out mycert.pem
Generating a 1024 bit RSA private key
......++++++
........................++++++
writing new private key to 'mycert.pem'
-----
You are about to be asked to enter information that will be incorporated
into your certificate request.
What you are about to enter is what is called a Distinguished Name or a DN.
There are quite a few fields but you can leave some blank
For some fields there will be a default value,
If you enter '.', the field will be left blank.
-----
Country Name (2 letter code) [AU]:
State or Province Name (full name) [Some-State]:
Locality Name (eg, city) []:
Organization Name (eg, company) [Internet Widgits Pty Ltd]:
Organizational Unit Name (eg, section) []:
Common Name (eg, YOUR name) []:
Email Address []:

openssl工具提示您填写一些字段。有关更多信息，请查看相关的手册页（手册页的缩写），如下所示：

$ man openssl

使用以下命令为服务器创建特殊的配置文件：
```
$ ipython profile create nbserver
```
编辑配置文件。在此示例中，可以在~/.ipython/profile_nbserver/ipython_notebook_config.py中找到。

配置文件很大，因此我们将省略其中的大多数行。我们至少需要更改的行如下：
```
c.NotebookApp.certfile = u'/absolute/path/to/your/certificate'
c.NotebookApp.password = u'sha1:b...your password'
c.NotebookApp.port = 9999
```
请注意，我们指向的是我们创建的 SSL 证书。我们设置了密码，并将端口更改为 9999。

使用以下命令，启动服务器以检查更改是否有效：

$ ipython notebook --profile=nbserver
[NotebookApp] Using existing profile dir: u'/Users/ivanidris/.ipython/profile_nbserver'
[NotebookApp] The IPython Notebook is running at: https://127.0.0.1:9999
[NotebookApp] Use Control-C to stop this server and shut down all kernels.

服务器在端口 9999 上运行，您需要通过 HTTPS 连接到它。如果一切顺利，您应该会看到一个登录页面。另外，您可能需要在浏览器中接受安全例外。

工作原理

我们为公共服务器创建了一个特殊的配置文件。已经有一些样本概要文件，例如默认概要文件。创建配置文件后，将使用配置文件将profile_<profilename>文件夹添加到.ipython目录。然后可以使用--profile=<profile_name>命令行选项加载配置文件。我们可以使用以下命令列出配置文件：

$ ipython profile list

Available profiles in IPython:
 cluster
 math
 pysh
 python3

 The first request for a bundled profile will copy it
 into your IPython directory (/Users/ivanidris/.ipython),
 where you can customize it.

Available profiles in /Users/ivanidris/.ipython:
 default
 nbserver
 sh

另见

探索 SymPy 配置文件

IPython 有一个示例 SymPy 配置文件。 SymPy 是一个 Python 符号数学库。我们可以简化代数表达式或区分函数，类似于 Mathematica 和 Maple。显然，SymPy 是一款有趣的软件，但对于我们穿越 NumPy 的过程而言并不是必需的。将此视为可选或额外的秘籍。就像甜点一样，可以随意跳过它，尽管您可能会错过本章中最甜的部分。

准备

使用easy_install或pip安装 SymPy：

$ sudo easy_install sympy
$ sudo pip install sympy

操作步骤

以下步骤将帮助您探索 SymPy 配置文件：

查看配置文件，该文件位于~/.ipython/profile_sympy/ipython_config.py。内容如下：

c = get_config()
app = c.InteractiveShellApp

# This can be used at any point in a config file to load a sub config
# and merge it into the current one.
load_subconfig('ipython_config.py', profile='default')

lines = """
from __future__ import division
from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
"""

# You have to make sure that attributes that are containers already
# exist before using them.  Simple assigning a new list will override
# all previous values.

if hasattr(app, 'exec_lines'):
 app.exec_lines.append(lines)
else:
 app.exec_lines = [lines]

# Load the sympy_printing extension to enable nice printing of sympy expr's.
if hasattr(app, 'extensions'):
 app.extensions.append('sympyprinting')
else:
 app.extensions = ['sympyprinting']

此代码完成以下任务：

加载默认配置文件
导入 SymPy 包
定义符号

使用以下命令以 SymPy 配置文件启动 IPython：
```
$ ipython --profile=sympy
```
使用以下屏幕快照中所示的命令扩展代数表达式：

另见

SymPy 主页

二、高级索引和数组概念

在本章中，我们将介绍以下秘籍：

安装 SciPy
安装 PIL
调整图像大小
比较视图和副本
翻转 Lena
花式索引
位置列表索引
布尔值索引
数独的步幅技巧
广播数组

简介

NumPy 以其高效的数组而闻名。之所以成名，部分原因是索引容易。我们将演示使用图像的高级索引技巧。在深入研究索引之前，我们将安装必要的软件 -- SciPy 和 PIL。如果您认为有此需要，请参阅第 1 章“使用 IPython”的“安装 matplotlib”秘籍。

在本章和其他章中，我们将使用以下导入：

import numpy as np 
import matplotlib.pyplot as plt
import scipy

我们还将尽可能为print() Python 函数使用最新的语法。

注意

Python2 是仍然很流行的主要 Python 版本，但与 Python3 不兼容。Python2 直到 2020 年才正式失去支持。主要区别之一是print()函数的语法。本书使用的代码尽可能与 Python2 和 Python3 兼容。

本章中的一些示例涉及图像处理。为此，我们将需要 Python 图像库（PIL），但不要担心；必要时会在本章中提供帮助您安装 PIL 和其他必要 Python 软件的说明和指示。

安装 SciPy

SciPy 是科学的 Python 库，与 NumPy 密切相关。实际上，SciPy 和 NumPy 在很多年前曾经是同一项目。就像 NumPy 一样，SciPy 是一个开放源代码项目，已获得 BSD 许可。在此秘籍中，我们将安装 SciPy。 SciPy 提供高级功能，包括统计，信号处理，线性代数，优化，FFT，ODE 求解器，插值，特殊功能和积分。 NumPy 有一些重叠，但是 NumPy 主要提供数组功能。

准备

在第 1 章，“使用 IPython”中，我们讨论了如何安装setuptools和pip。如有必要，请重新阅读秘籍。

操作步骤

在此秘籍中，我们将完成安装 SciPy 的步骤：

从源安装：如果已安装 Git，则可以使用以下命令克隆 SciPy 存储库：
```
$ git clone https://github.com/scipy/scipy.git

$ python setup.py build
$ python setup.py install --user
```
这会将 SciPy 安装到您的主目录。它需要 Python 2.6 或更高版本。

在构建之前，您还需要安装 SciPy 依赖的以下包：
- BLAS和LAPACK库
- C 和 Fortran 编译器
您可能已经在 NumPy 安装过程中安装了此软件。
在 Linux 上安装 SciPy：大多数 Linux 发行版都包含 SciPy 包。我们将遵循一些流行的 Linux 发行版中的必要步骤（您可能需要以 root 用户身份登录或具有sudo权限）：
- 为了在 RedHat，Fedora 和 CentOS 上安装 SciPy，请从命令行运行以下指令：
```
$ yum install python-scipy
```
- 为了在 Mandriva 上安装 SciPy，请运行以下命令行指令：
```
$ urpmi python-scipy
```
- 为了在 Gentoo 上安装 SciPy，请运行以下命令行指令：
```
$ sudo emerge scipy
```
- 在 Debian 或 Ubuntu 上，我们需要输入以下指令：
```
$ sudo apt-get install python-scipy
```
在 MacOSX 上安装 SciPy：需要 Apple Developer Tools（XCode），因为它包含BLAS和LAPACK库。可以在 App Store 或 Mac 随附的安装 DVD 中找到它。或者您可以从 Apple Developer 的连接网站获取最新版本。确保已安装所有内容，包括所有可选包。

您可能已经为 NumPy 安装了 Fortran 编译器。 gfortran的二进制文件可以在这个链接中找到。
使用easy_install或pip安装 SciPy：您可以使用以下两个命令中的任何一个来安装 SciPy（sudo的需要取决于权限）：
```
$ [sudo] pip install scipy
$ [sudo] easy_install scipy

```** 
```
在 Windows 上安装：如果已经安装 Python，则首选方法是下载并使用二进制发行版。或者，您可以安装 Anaconda 或 Enthought Python 发行版，该发行版与其他科学 Python 包一起提供。
检查安装：使用以下代码检查 SciPy 安装：
```
import scipy
print(scipy.__version__)
print(scipy.__file__)
```
这应该打印正确的 SciPy 版本。

工作原理

大多数包管理器都会为您解决依赖项（如果有）。但是，在某些情况下，您需要手动安装它们。这超出了本书的范围。

另见

如果遇到问题，可以在以下位置寻求帮助：

freenode的#scipy IRC 频道
SciPy 邮件列表

安装 PIL

PIL（Python 图像库）是本章中进行图像处理的先决条件。如果愿意，可以安装 Pillow，它是 PIL 的分支。有些人喜欢 Pillow API；但是，我们不会在本书中介绍其安装。

操作步骤

让我们看看如何安装 PIL：

在 Windows 上安装 PIL：使用 Windows 中的 PIL 可执行文件安装 PIL。
在 Debian 或 Ubuntu 上安装：在 Debian 或 Ubuntu 上，使用以下命令安装 PIL：
```
$ sudo apt-get install python-imaging
```
使用easy_install或pip安装：在编写本书时，似乎 RedHat，Fedora 和 CentOS 的包管理器没有对 PIL 的直接支持。因此，如果您使用的是这些 Linux 发行版之一，请执行此步骤。

使用以下任一命令安装：
```
$ easy_install PIL
$ sudo pip install PIL
```

另见

可在这里找到有关 PILLOW（PIL 的分支）的说明。

调整图像大小

在此秘籍中，我们将把 Lena 的样例图像（在 SciPy 发行版中可用）加载到数组中。顺便说一下，本章不是关于图像操作的。我们将只使用图像数据作为输入。

注意

Lena Soderberg 出现在 1972 年的《花花公子》杂志中。由于历史原因，这些图像之一经常用于图像处理领域。不用担心，该图像完全可以安全工作。

我们将使用repeat()函数调整图像大小。此函数重复一个数组，这意味着在我们的用例中按一定的大小调整图像大小。

准备

此秘籍的前提条件是必须安装 SciPy，matplotlib 和 PIL。看看本章和第 1 章，“使用 IPython”的相应秘籍。

操作步骤

通过以下步骤调整图像大小：

首先，导入SciPy。 SciPy 具有lena()函数。它用于将图像加载到 NumPy 数组中：
```
lena = scipy.misc.lena()
```
从 0.10 版本开始发生了一些重构，因此，如果您使用的是旧版本，则正确的代码如下：
```
lena = scipy.lena()
```
使用numpy.testing包中的assert_equal()函数检查 Lena 数组的形状-这是可选的完整性检查测试：
```
np.testing.assert_equal((LENA_X, LENA_Y), lena.shape)
```
使用repeat()函数调整 Lena 数组的大小。我们在x和y方向上给此函数一个调整大小的因子：
```
resized = lena.repeat(yfactor, axis=0).repeat(xfactor, axis=1)
```

我们将在同一网格的两个子图中绘制 Lena 图像和调整大小后的图像。使用以下代码在子图中绘制 Lena 数组：

plt.subplot(211)
plt.title("Lena")
plt.axis("off")
plt.imshow(lena)

matplotlib subplot()函数创建一个子图。此函数接受一个三位整数作为参数，其中第一位是行数，第二位是列数，最后一位是子图的索引，从 1 开始。imshow()函数显示图像。最后，show()函数显示最终结果。

将调整大小后的数组绘制在另一个子图中并显示它。索引现在为 2：

plt.subplot(212)
plt.title("Resized")
plt.axis("off")
plt.imshow(resized)
plt.show()

以下屏幕截图显示了结果，以及原始图像（第一幅）和调整大小后的图像（第二幅）：

以下是本书代码包中resize_lena.py文件中该秘籍的完整代码：

import scipy.misc
import matplotlib.pyplot as plt
import numpy as np

# This script resizes the Lena image from Scipy.

# Loads the Lena image into an array
lena = scipy.misc.lena()

#Lena's dimensions
LENA_X = 512
LENA_Y = 512

#Check the shape of the Lena array
np.testing.assert_equal((LENA_X, LENA_Y), lena.shape)

# Set the resize factors
yfactor = 2
xfactor = 3

# Resize the Lena array
resized = lena.repeat(yfactor, axis=0).repeat(xfactor, axis=1)

#Check the shape of the resized array
np.testing.assert_equal((yfactor * LENA_Y, xfactor * LENA_Y), resized.shape)

# Plot the Lena array
plt.subplot(211)
plt.title("Lena")
plt.axis("off")
plt.imshow(lena)

#Plot the resized array
plt.subplot(212)
plt.title("Resized")
plt.axis("off")
plt.imshow(resized)
plt.show()

工作原理

repeat()函数重复数组，在这种情况下，这会导致原始图像的大小改变。 subplot() matplotlib 函数创建一个子图。 imshow()函数显示图像。最后，show()函数显示最终结果。

另见

第 1 章“使用 IPython”中的“安装 matplotlib”
本章中的“安装 SciPy”
本章中的“安装 PIL”
这个页面中介绍了repeat()函数。

创建视图和副本

了解何时处理共享数组视图以及何时具有数组数据的副本，这一点很重要。例如，切片将创建一个视图。这意味着，如果您将切片分配给变量，然后更改基础数组，则此变量的值将更改。我们将根据著名的 Lena 图像创建一个数组，复制该数组，创建一个视图，最后修改视图。

准备

前提条件与先前的秘籍相同。

操作步骤

让我们创建 Lena 数组的副本和视图：

创建 Lena 数组的副本：
```
acopy = lena.copy()
```
创建数组的视图：
```
aview = lena.view()
```

使用flat迭代器将视图的所有值设置为0：

aview.flat = 0

最终结果是只有一个图像（与副本相关的图像）显示了花花公子模型。其他图像完全消失：

以下是本教程的代码，显示了本书代码包中copy_view.py文件中数组视图和副本的行为：

import scipy.misc
import matplotlib.pyplot as plt

lena = scipy.misc.lena()
acopy = lena.copy()
aview = lena.view()

# Plot the Lena array
plt.subplot(221)
plt.imshow(lena)

#Plot the copy
plt.subplot(222)
plt.imshow(acopy)

#Plot the view
plt.subplot(223)
plt.imshow(aview)

# Plot the view after changes
aview.flat = 0
plt.subplot(224)
plt.imshow(aview)

plt.show()

工作原理

如您所见，通过在程序结尾处更改视图，我们更改了原始 Lena 数组。这样就产生了三张蓝色（如果您正在查看黑白图像，则为空白）图像-复制的数组不受影响。重要的是要记住，视图不是只读的。

另见

NumPy view()函数的文档位于这里

翻转 Lena

我们将翻转 SciPy Lena 图像-当然，所有这些都是以科学的名义，或者至少是作为演示。除了翻转图像，我们还将对其进行切片并对其应用遮罩。

操作步骤

步骤如下：

使用以下代码围绕垂直轴翻转 Lena 数组：
```
plt.imshow(lena[:,::-1])
```
从图像中切出一部分并将其绘制出来。在这一步中，我们将看一下 Lena 数组的形状。该形状是表示数组大小的元组。以下代码有效地选择了花花公子图片的左上象限：
```
plt.imshow(lena[:lena.shape[0]/2,:lena.shape[1]/2])
```

通过在 Lena 数组中找到所有偶数的值，对图像应用遮罩（这对于演示目的来说是任意的）。复制数组并将偶数值更改为 0。这样会在图像上放置很多蓝点（如果您正在查看黑白图像，则会出现暗点）：

mask = lena % 2 == 0
masked_lena = lena.copy()
masked_lena[mask] = 0

所有这些工作都会产生2 x 2的图像网格，如以下屏幕截图所示：

这是本书代码包中flip_lena.py文件中此秘籍的完整代码：

import scipy.misc
import matplotlib.pyplot as plt

# Load the Lena array
lena = scipy.misc.lena()

# Plot the Lena array
plt.subplot(221)
plt.title('Original')
plt.axis('off')
plt.imshow(lena)

#Plot the flipped array
plt.subplot(222)
plt.title('Flipped')
plt.axis('off')
plt.imshow(lena[:,::-1])

#Plot a slice array
plt.subplot(223)
plt.title('Sliced')
plt.axis('off')
plt.imshow(lena[:lena.shape[0]/2,:lena.shape[1]/2])

# Apply a mask
mask = lena % 2 == 0
masked_lena = lena.copy()
masked_lena[mask] = 0
plt.subplot(224)
plt.title('Masked')
plt.axis('off')
plt.imshow(masked_lena)

plt.show()

另见

第 1 章“使用 IPython”中的“安装 matplotlib”
本章中的“安装 SciPy”
本章中的“安装 PIL”

花式索引

在本教程中，我们将应用花式索引将 Lena 图像的对角线值设置为 0。这将沿着对角线绘制黑线并交叉，这不是因为图像有问题，而仅仅作为练习。花式索引是不涉及整数或切片的索引；这是正常的索引编制。

操作步骤

我们将从第一个对角线开始：

将第一个对角线的值设置为0。

要将对角线值设置为0，我们需要为x和y值定义两个不同的范围：
```
lena[range(xmax), range(ymax)] = 0
```

将另一个对角线的值设置为0。

要设置另一个对角线的值，我们需要一组不同的范围，但是原理保持不变：

lena[range(xmax-1,-1,-1), range(ymax)] = 0

最后，我们得到带有对角线标记的图像，如以下屏幕截图所示：

以下是本书代码集中fancy.py文件中该秘籍的完整代码：

import scipy.misc
import matplotlib.pyplot as plt

# This script demonstrates fancy indexing by setting values
# on the diagonals to 0.

# Load the Lena array
lena = scipy.misc.lena()
xmax = lena.shape[0]
ymax = lena.shape[1]

# Fancy indexing
# Set values on diagonal to 0
# x 0-xmax
# y 0-ymax
lena[range(xmax), range(ymax)] = 0

# Set values on other diagonal to 0
# x xmax-0
# y 0-ymax
lena[range(xmax-1,-1,-1), range(ymax)] = 0

# Plot Lena with diagonal lines set to 0
plt.imshow(lena)
plt.show()

工作原理

我们为x值和y值定义了单独的范围。这些范围用于索引 Lena 数组。花式索引是基于内部 NumPy 迭代器对象执行的。执行以下步骤：

创建迭代器对象。
迭代器对象绑定到数组。
数组元素通过迭代器访问。

另见

花式索引的实现文档

位置列表索引

让我们使用ix_()函数来随机播放 Lena 图像。此函数根据多个序列创建网格。

操作步骤

我们将从随机改组数组索引开始：

使用numpy.random模块的shuffle()函数创建随机索引数组：

def shuffle_indices(size):
   arr = np.arange(size)
   np.random.shuffle(arr)

   return arr

绘制乱序的索引：

plt.imshow(lena[np.ix_(xindices, yindices)])

我们得到的是一张完全打乱的 Lena 图像，如以下屏幕截图所示：

这是本书代码包中ix.py文件中秘籍的完整代码：

import scipy.misc
import matplotlib.pyplot as plt
import numpy as np

# Load the Lena array
lena = scipy.misc.lena()
xmax = lena.shape[0]
ymax = lena.shape[1]

def shuffle_indices(size):
   '''
   Shuffles an array with values 0 - size
   '''
   arr = np.arange(size)
   np.random.shuffle(arr)

   return arr

xindices = shuffle_indices(xmax)
np.testing.assert_equal(len(xindices), xmax)
yindices = shuffle_indices(ymax)
np.testing.assert_equal(len(yindices), ymax)

# Plot Lena
plt.imshow(lena[np.ix_(xindices, yindices)])
plt.show()

另见

ix_()函数的文档页面

布尔值索引

布尔索引是基于布尔数组的索引，属于奇特索引的类别。

操作步骤

我们将这种索引技术应用于图像：

在对角线上带有点的图像。

这在某种程度上类似于本章中的“花式索引”秘籍。这次，我们在图像的对角线上选择模4：

def get_indices(size):
   arr = np.arange(size)
   return arr % 4 == 0

然后，我们只需应用此选择并绘制点：

lena1 = lena.copy() 
xindices = get_indices(lena.shape[0])
yindices = get_indices(lena.shape[1])
lena1[xindices, yindices] = 0
plt.subplot(211)
plt.imshow(lena1)

在最大值的四分之一到四分之三之间选择数组值，并将它们设置为0：

lena2[(lena > lena.max()/4) & (lena < 3 * lena.max()/4)] = 0

带有两个新图像的图看起来类似于以下屏幕截图所示：

这是本书代码包中boolean_indexing.py文件中该秘籍的完整代码：

import scipy.misc
import matplotlib.pyplot as plt
import numpy as np

# Load the Lena array
lena = scipy.misc.lena()

def get_indices(size):
   arr = np.arange(size)
   return arr % 4 == 0

# Plot Lena
lena1 = lena.copy() 
xindices = get_indices(lena.shape[0])
yindices = get_indices(lena.shape[1])
lena1[xindices, yindices] = 0
plt.subplot(211)
plt.imshow(lena1)

lena2 = lena.copy() 
# Between quarter and 3 quarters of the max value
lena2[(lena > lena.max()/4) & (lena < 3 * lena.max()/4)] = 0
plt.subplot(212)
plt.imshow(lena2)

plt.show()

工作原理

由于布尔值索引是一种花式索引，因此它的工作方式基本相同。这意味着索引是在特殊的迭代器对象的帮助下发生的。

另见

“花式索引”

数独的步幅技巧

ndarray 类具有strides字段，它是一个元组，指示通过数组时要在每个维中步进的字节数。让我们对将数独谜题拆分为3 x 3正方形的问题应用一些大步技巧。

注意

对数独的规则进行解释超出了本书的范围。简而言之，数独谜题由3 x 3的正方形组成。这些正方形均包含九个数字。有关更多信息，请参见这里。

操作步骤

应用如下的跨步技巧：

让我们定义sudoku数组。此数组充满了一个实际的已解决的数独难题的内容：

sudoku = np.array([
    [2, 8, 7, 1, 6, 5, 9, 4, 3],
    [9, 5, 4, 7, 3, 2, 1, 6, 8],
    [6, 1, 3, 8, 4, 9, 7, 5, 2],
    [8, 7, 9, 6, 5, 1, 2, 3, 4],
    [4, 2, 1, 3, 9, 8, 6, 7, 5],
    [3, 6, 5, 4, 2, 7, 8, 9, 1],
    [1, 9, 8, 5, 7, 3, 4, 2, 6],
    [5, 4, 2, 9, 1, 6, 3, 8, 7],
    [7, 3, 6, 2, 8, 4, 5, 1, 9]
    ])

ndarray的itemsize字段为我们提供了数组中的字节数。给定itemsize，请计算步幅：
```
strides = sudoku.itemsize * np.array([27, 3, 9, 1])
```

现在我们可以使用np.lib.stride_tricks模块的as_strided()函数将拼图分解成正方形：

squares = np.lib.stride_tricks.as_strided(sudoku, shape=shape, strides=strides)
print(squares)

该代码打印单独的数独正方形，如下所示：

[[[[2 8 7]
    [9 5 4]
    [6 1 3]]

  [[1 6 5]
    [7 3 2]
    [8 4 9]]

  [[9 4 3]
    [1 6 8]
    [7 5 2]]]

 [[[8 7 9]
    [4 2 1]
    [3 6 5]]

  [[6 5 1]
    [3 9 8]
    [4 2 7]]

  [[2 3 4]
    [6 7 5]
    [8 9 1]]]

 [[[1 9 8]
    [5 4 2]
    [7 3 6]]

  [[5 7 3]
    [9 1 6]
    [2 8 4]]

  [[4 2 6]
    [3 8 7]
    [5 1 9]]]]

以下是本书代码包中strides.py文件中此秘籍的完整源代码：

import numpy as np

sudoku = np.array([
   [2, 8, 7, 1, 6, 5, 9, 4, 3],
   [9, 5, 4, 7, 3, 2, 1, 6, 8],
   [6, 1, 3, 8, 4, 9, 7, 5, 2],
   [8, 7, 9, 6, 5, 1, 2, 3, 4],
   [4, 2, 1, 3, 9, 8, 6, 7, 5],
   [3, 6, 5, 4, 2, 7, 8, 9, 1],
   [1, 9, 8, 5, 7, 3, 4, 2, 6],
   [5, 4, 2, 9, 1, 6, 3, 8, 7],
   [7, 3, 6, 2, 8, 4, 5, 1, 9]
   ])

shape = (3, 3, 3, 3)

strides = sudoku.itemsize * np.array([27, 3, 9, 1])

squares = np.lib.stride_tricks.as_strided(sudoku, shape=shape, strides=strides)
print(squares)

工作原理

我们应用了跨步技巧，将数独谜题拆分为3 x 3的正方形。步幅告诉我们通过数独数组时每一步需要跳过的字节数。

另见

strides属性的文档在这里

广播数组

在不知道的情况下，您可能已经广播了数组。简而言之，即使操作数的形状不同，NumPy 也会尝试执行操作。在此秘籍中，我们将一个数组和一个标量相乘。标量被扩展为数组操作数的形状，然后执行乘法。我们将下载一个音频文件并制作一个更安静的新版本。

操作步骤

让我们从读取 WAV 文件开始：

我们将使用标准的 Python 代码下载 Austin Powers 的音频文件。 SciPy 具有 WAV 文件模块，可让您加载声音数据或生成 WAV 文件。如果已安装 SciPy，则我们应该已经有此模块。 read()函数返回data数组和采样率。在此示例中，我们仅关心数据：
```
sample_rate, data = scipy.io.wavfile.read(WAV_FILE)
```
使用 matplotlib 绘制原始 WAV 数据。将子图命名为Original：
```
plt.subplot(2, 1, 1)
plt.title("Original")
plt.plot(data)
```
现在，我们将使用 NumPy 制作更安静的音频样本。这只是通过与常量相乘来创建具有较小值的新数组的问题。这就是广播魔术发生的地方。最后，由于 WAV 格式，我们需要确保与原始数组具有相同的数据类型：
```
newdata = data * 0.2
newdata = newdata.astype(np.uint8)
```

可以将新数组写入新的 WAV 文件，如下所示：

scipy.io.wavfile.write("quiet.wav",
    sample_rate, newdata)

使用 matplotlib 绘制新数据数组：

plt.subplot(2, 1, 2)
plt.title("Quiet")
plt.plot(newdata)

plt.show()

结果是原始 WAV 文件数据和具有较小值的新数组的图，如以下屏幕快照所示：

这是本书代码包中broadcasting.py文件中该秘籍的完整代码：

import scipy.io.wavfile
import matplotlib.pyplot as plt
import urllib2
import numpy as np

# Download audio file
response = urllib2.urlopen('http://www.thesoundarchive.com/austinpowers/smashingbaby.wav')
print(response.info())
WAV_FILE = 'smashingbaby.wav'
filehandle = open(WAV_FILE, 'w')
filehandle.write(response.read())
filehandle.close()
sample_rate, data = scipy.io.wavfile.read(WAV_FILE)
print("Data type", data.dtype, "Shape", data.shape)

# Plot values original audio
plt.subplot(2, 1, 1)
plt.title("Original")
plt.plot(data)

# Create quieter audio
newdata = data * 0.2
newdata = newdata.astype(np.uint8)
print("Data type", newdata.dtype, "Shape", newdata.shape)

# Save quieter audio file
scipy.io.wavfile.write("quiet.wav",
    sample_rate, newdata)

# Plot values quieter file
plt.subplot(2, 1, 2)
plt.title("Quiet")
plt.plot(newdata)

plt.show()

另见

以下链接提供了更多背景信息：

三、掌握常用函数

在本章中，我们将介绍许多常用函数：

sqrt()，log()，arange()，astype()和sum()
ceil()，modf()，where()，ravel()和take()
sort()和outer()
diff()，sign()和eig()
histogram()和polyfit()
compress()和randint()

我们将在以下秘籍中讨论这些功能：

斐波纳契数求和
查找素因数
查找回文数
稳态向量
发现幂律
逢低定期交易
随机模拟交易
用 Eratosthenes 筛子来筛选质数

简介

本章介绍常用的 NumPy 函数。这些是您每天将要使用的函数。显然，用法可能与您不同。 NumPy 函数太多，以至于几乎不可能全部了解，但是本章中的函数是我们应该熟悉的最低要求。

斐波纳契数求和

在此秘籍中，我们将求和值不超过 400 万的斐波纳契数列中的偶数项。斐波那契数列是从零开始的整数序列，其中每个数字都是前两个数字的和，但（当然）前两个数字除外，零和一（0、1、1、2、3、5、8、13、21、34、55、89 ...）。

该序列由斐波那契（Fibonacci）在 1202 年发布，最初不包含零。实际上，早在几个世纪以前，印度数学家就已经知道了它。斐波那契数在数学，计算机科学和生物学中都有应用。

注意

有关更多信息，请阅读 Wikipedia 关于斐波那契数字的文章。

该秘籍使用基于黄金比例的公式，这是一个无理数，具有与pi相当的特殊性质。黄金比例由以下公式给出：

我们将使用sqrt()，log()，arange()，astype()和sum()函数。斐波那契数列的递归关系具有以下解，涉及黄金比率：

操作步骤

以下是本书代码包中sum_fibonacci.py文件中此秘籍的完整代码：

import numpy as np

#Each new term in the Fibonacci sequence is generated by adding the previous two terms.
#By starting with 1 and 2, the first 10 terms will be:

#1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

#By considering the terms in the Fibonacci sequence whose values do not exceed four million,
#find the sum of the even-valued terms.

#1\. Calculate phi
phi = (1 + np.sqrt(5))/2
print("Phi", phi)

#2\. Find the index below 4 million
n = np.log(4 * 10 ** 6 * np.sqrt(5) + 0.5)/np.log(phi)
print(n)

#3\. Create an array of 1-n
n = np.arange(1, n)
print(n)

#4\. Compute Fibonacci numbers
fib = (phi**n - (-1/phi)**n)/np.sqrt(5)
print("First 9 Fibonacci Numbers", fib[:9])

#5\. Convert to integers
# optional
fib = fib.astype(int)
print("Integers", fib)

#6\. Select even-valued terms
eventerms = fib[fib % 2 == 0]
print(eventerms)

#7\. Sum the selected terms
print(eventerms.sum())

的第一件事是计算黄金分割率，也称为黄金分割或黄金平均值。

使用sqrt()函数计算5的平方根：

phi = (1 + np.sqrt(5))/2
print("Phi", phi)

这印出了中庸之道：

Phi 1.61803398875

接下来，在秘籍中，我们需要找到低于 400 万的斐波那契数的指数。维基百科页面中提供了一个公式，我们将使用该公式进行计算。我们需要做的就是使用log()函数转换对数。我们不需要将结果四舍五入为最接近的整数。在秘籍的下一步中，这将自动为我们完成：
```
n = np.log(4 * 10 ** 6 * np.sqrt(5) + 0.5)/np.log(phi)
print(n)
```
n的值如下：
```
33.2629480359
```
arange()函数是许多人都知道的非常基本的函数。不过，出于完整性考虑，我们将在这里提及：
```
n = np.arange(1, n)
```
我们可以使用一个方便的公式来计算斐波那契数。我们将需要黄金比例和该秘籍中上一步中的数组作为输入参数。打印前九个斐波那契数字以检查结果：
```
fib = (phi**n - (-1/phi)**n)/np.sqrt(5)
print("First 9 Fibonacci Numbers", fib[:9])
```
注意

我本可以进行单元测试而不是打印声明。单元测试是测试一小段代码（例如函数）的测试。秘籍的这种变化是您的练习。

提示

查看第 8 章，“质量保证”，以获取有关如何编写单元测试的指针。

顺便说一下，我们不是从数字 0 开始的。上面的代码给了我们一系列预期的结果：
```
First 9 Fibonacci Numbers [  1\.   1\.   2\.   3\.   5\.   8\.  13\.  21\.  34.]
```
您可以根据需要将此权限插入单元测试。
转换为整数。

此步骤是可选的。我认为最后有一个整数结果是很好的。好的，我实际上想向您展示astype()函数：
```
fib = fib.astype(int)
print("Integers", fib)
```
为简短起见，此代码为我们提供了以下结果：
```
Integers [      1       1       2       3       5       8      13      21      34
 ... snip ... snip ...
 317811  514229  832040 1346269 2178309 3524578]
```
选择偶数项。

此秘籍要求我们现在选择偶数项。如果遵循第 2 章，“高级索引和数组概念”中的“布尔值索引”秘籍，这对您来说应该很容易：
```
eventerms = fib[fib % 2 == 0]
print(eventerms)
```
我们去了：
```
[      2       8      34     144     610    2584   10946   46368  196418  832040 3524578]
```

工作原理

在此秘籍中，我们使用了sqrt()，log()，arange()，astype()和sum()函数。其描述如下：

函数	描述
`sqrt()`	此函数计算数组元素的平方根
`log()`	此函数计算数组元素的自然对数
`arange()`	此函数创建具有指定范围的数组
`astype()`	此函数将数组元素转换为指定的数据类型
`sum()`	此函数计算数组元素的总和

另见

第 2 章，“高级索引和数组概念”中的“布尔值索引”秘籍

查找素因数

素因数是质数，它们精确地除以整数而不会留下余数。对于较大的数字，找到主要因子似乎几乎是不可能的。因此，素因数在密码学中具有应用。但是，使用正确的算法 -- Fermat 因式分解方法和 NumPy -- 对于小数而言，因式分解变得相对容易。想法是将N分解为两个数字，c和d，根据以下等式：

我们可以递归应用因式分解，直到获得所需的素因子。

操作步骤

以下是解决找到最大质数因子 600851475143 的问题所需的全部代码（请参见本书代码包中的fermatfactor.py文件）：

from __future__ import print_function
import numpy as np

#The prime factors of 13195 are 5, 7, 13 and 29.

#What is the largest prime factor of the number 600851475143 ?

N = 600851475143
LIM = 10 ** 6

def factor(n):
   #1\. Create array of trial values
   a = np.ceil(np.sqrt(n))
   lim = min(n, LIM)
   a = np.arange(a, a + lim)
   b2 = a ** 2 - n

   #2\. Check whether b is a square
   fractions = np.modf(np.sqrt(b2))[0]

   #3\. Find 0 fractions
   indices = np.where(fractions == 0)

   #4\. Find the first occurrence of a 0 fraction
   a = np.ravel(np.take(a, indices))[0]
              # Or a = a[indices][0]

   a = int(a)
   b = np.sqrt(a ** 2 - n) 
   b = int(b)
   c = a + b
   d = a - b

   if c == 1 or d == 1:
      return

   print(c, d)
   factor(c)
   factor(d)

factor(N)

该算法要求我们为a尝试一些试验值：

创建试验值的数组。

创建一个 NumPy 数组并消除循环需求是有意义的。但是，应注意不要创建一个在内存需求方面太大的数组。在我的系统上，一百万个元素的数组似乎正好合适：
```
a = np.ceil(np.sqrt(n))
lim = min(n, LIM)
a = np.arange(a, a + lim)
b2 = a ** 2 - n
```
我们使用ceil()函数以元素为单位返回输入的上限。
获取b数组的小数部分。

现在我们应该检查b是否为正方形。使用 NumPy modf()函数获取b数组的分数部分：
```
fractions = np.modf(np.sqrt(b2))[0]
```
查找0分数。

调用where() NumPy 函数以找到零分数的索引，其中小数部分是0：
```
indices = np.where(fractions == 0)
```
找到零分数的第一个出现。

首先，使用上一步中的indices数组调用take() NumPy 函数，以获取零分数的值。现在，使用ravel() NumPy 函数将这个数组变得扁平：
```
a = np.ravel(np.take(a, indices))[0]
```
提示

这条线有些令人费解，但是确实演示了两个有用的功能。写a = a[indices][0]会更简单。

此代码的输出如下：
```
1234169 486847
1471 839
6857 71
```

工作原理

我们使用ceil()，modf()，where()，ravel()和take() NumPy 函数递归地应用了费马分解。这些函数的说明如下：

函数	描述
`ceil()`	计算数组元素的上限
`modf()`	返回浮点数数字的分数和整数部分
`where()`	根据条件返回数组索引
`ravel()`	返回一个扁平数组
`take()`	从数组中获取元素

查找回文数

回文数字在两种方式下的读取相同。由两个 2 位数字的乘积组成的最大回文为9009 = 91 x 99。让我们尝试查找由两个 3 位数字的乘积组成的最大回文。

操作步骤

以下是本书代码包中palindromic.py文件的完整程序：

import numpy as np

#A palindromic number reads the same both ways. 
#The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 x 99.

#Find the largest palindrome made from the product of two 3-digit numbers.

#1\. Create  3-digits numbers array
a = np.arange(100, 1000)
np.testing.assert_equal(100, a[0])
np.testing.assert_equal(999, a[-1])

#2\. Create products array
numbers = np.outer(a, a)
numbers = np.ravel(numbers)
numbers.sort()
np.testing.assert_equal(810000, len(numbers))
np.testing.assert_equal(10000, numbers[0])
np.testing.assert_equal(998001, numbers[-1])

#3\. Find largest palindromic number
for number in numbers[::-1]:
   s = str(numbers[i])

   if s == s[::-1]:
      print(s)
      break

我们将使用最喜欢的 NumPy 函数arange()创建一个数组，以容纳从 100 到 999 的三位数。

创建一个三位数的数字数组。

使用numpy.testing包中的assert_equal()函数检查数组的第一个和最后一个元素：
```
a = np.arange(100, 1000)
np.testing.assert_equal(100, a[0])
np.testing.assert_equal(999, a[-1])
```
创建乘积数组。

现在，我们将创建一个数组，以将三位数数组的元素的所有可能乘积与其自身保持在一起。我们可以使用outer()函数来完成此操作。需要使用ravel()将生成的数组弄平，以便能够轻松地对其进行迭代。在数组上调用sort()方法，以确保数组正确排序。之后，我们可以进行一些检查：
```
numbers = np.outer(a, a)
numbers = np.ravel(numbers)
numbers.sort()
np.testing.assert_equal(810000, len(numbers))
np.testing.assert_equal(10000, numbers[0])
np.testing.assert_equal(998001, numbers[-1])
```

该代码打印 906609，它是回文数。

工作原理

我们看到了outer()函数的作用。此函数返回两个数组的外部乘积。两个向量的外部乘积（一维数字列表）创建一个矩阵。这与内部乘积相反，该乘积返回两个向量的标量数。外部产品用于物理，信号处理和统计。 sort()函数返回数组的排序副本。

检查结果可能是一个好主意。稍微修改一下代码，找出哪两个 3 位数字产生我们的回文码。尝试以 NumPy 方式实现最后一步。

稳态向量

马尔科夫链是一个至少具有两个状态的系统。有关马尔可夫链的详细信息，请参阅这里。时间t的状态取决于时间t-1的状态，仅取决于t-1的状态。系统在这些状态之间随机切换。链没有关于状态的任何记忆。马尔可夫链通常用于对物理，化学，金融和计算机科学中的现象进行建模。例如，Google 的 PageRank 算法使用马尔可夫链对网页进行排名。

我想为股票定义一个马尔可夫链。假设状态为震荡，上涨和下跌的状态。我们可以根据日末收盘价确定稳定状态。

在遥远的未来，或理论上经过无限长的时间之后，我们的马尔可夫链系统的状态将不再改变。这称为稳定状态。动态平衡是一种稳态。对于股票而言，达到稳定状态可能意味着关联公司已变得稳定。 随机矩阵A包含状态转移概率，当应用于稳态时，它会产生相同的状态x。为此的数学符号如下：

解决此问题的另一种方法是特征值和特征向量。特征值和特征向量是线性代数的基本概念，并且在量子力学，机器学习和其他科学中应用。

操作步骤

以下是本书代码包中steady_state_vector.py文件中稳态向量示例的完整代码：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('AAPL', start, today)
close =  [q[4] for q in quotes]

states = np.sign(np.diff(close))

NDIM = 3
SM = np.zeros((NDIM, NDIM))

signs = [-1, 0, 1]
k = 1

for i, signi in enumerate(signs):
   #we start the transition from the state with the specified sign
   start_indices = np.where(states[:-1] == signi)[0]

   N = len(start_indices) + k * NDIM

   # skip since there are no transitions possible
   if N == 0:
      continue

   #find the values of states at the end positions
   end_values = states[start_indices + 1]

   for j, signj in enumerate(signs):
      # number of occurrences of this transition 
      occurrences = len(end_values[end_values == signj])
      SM[i][j] = (occurrences + k)/float(N)

print(SM)
eig_out = np.linalg.eig(SM)
print(eig_out)

idx_vec = np.where(np.abs(eig_out[0] - 1) < 0.1)
print("Index eigenvalue 1", idx_vec)

x = eig_out[1][:,idx_vec].flatten()
print("Steady state vector", x)
print("Check", np.dot(SM, x))

现在我们需要获取数据：

获取一年的数据。

一种实现方法是使用 matplotlib（请参阅第 1 章的“安装 matplotlib”秘籍，如有必要）。我们将检索去年的数据。这是执行此操作的代码：
```
today = date.today()
start = (today.year - 1, today.month, today.day)
quotes = quotes_historical_yahoo('AAPL', start, today)
```
选择收盘价。

现在，我们有了 Yahoo 金融的历史数据。数据表示为元组列表，但我们仅对收盘价感兴趣。

元组中的第一个元素代表日期。其次是开盘价，最高价，最低价和收盘价。最后一个元素是音量。我们可以选择以下收盘价：
```
close =  [q[4] for q in quotes]
```
收盘价是每个元组中的第五个数字。现在我们应该有大约 253 个收盘价的清单。
确定状态。

我们可以通过使用diff() NumPy 函数减去连续天的价格来确定状态。然后，通过差异的符号给出状态。 sign() NumPy 函数返回-1为负数，1为正数，否则返回0。
```
states = np.sign(np.diff(close))
```
将随机矩阵初始化为0值。

对于每个过渡，我们有三个可能的开始状态和三个可能的结束状态。例如，如果我们从启动状态开始，则可以切换到：
- 向上
- 平面
- 下
使用zeros() NumPy 函数初始化随机矩阵：
```
NDIM = 3
SM = np.zeros((NDIM, NDIM))
```
对于每个符号，选择相应的开始状态索引。

现在，代码变得有些混乱。我们将不得不使用实际的循环！我们将遍历所有可能的符号，并选择与每个符号相对应的开始状态索引。使用where() NumPy 函数选择索引。在这里，k是一个平滑常数，我们将在后面讨论：
```
signs = [-1, 0, 1]
k = 1

for i, signi in enumerate(signs):
   #we start the transition from the state with the specified sign
    start_indices = np.where(states[:-1] == signi)[0]
```
平滑和随机矩阵。

现在，我们可以计算每个过渡的出现次数。将其除以给定开始状态的跃迁总数，就可以得出随机矩阵的跃迁概率。顺便说一下，这不是最好的方法，因为它可能过度拟合。

在现实生活中，我们可能有一天收盘价不会发生变化，尽管对于流动性股票市场来说这不太可能。处理零出现的一种方法是应用加法平滑。这个想法是在我们发现的出现次数上增加一个常数，以消除零。以下代码计算随机矩阵的值：
```
N = len(start_indices) + k * NDIM

# skip since there are no transitions possible
if N == 0:
    continue

#find the values of states at the end positions
end_values = states[start_indices + 1]

for j, signj in enumerate(signs):
    # number of occurrences of this transition 
    occurrences = len(end_values[end_values == signj])
    SM[i][j] = (occurrences + k)/float(N)

print(SM)
```
前述代码所做的是基于出现次数和加性平滑计算每个可能过渡的过渡概率。在其中一个测试运行中，我得到了以下矩阵：
```
[[ 0.5047619   0.00952381  0.48571429]
 [ 0.33333333  0.33333333  0.33333333]
 [ 0.33774834  0.00662252  0.65562914]]
```

特征值和特征向量。

要获得特征值和特征向量，我们将需要linalg NumPy 模块和eig()函数：

eig_out = numpy.linalg.eig(SM)
print(eig_out)

eig()函数返回一个包含特征值的数组和另一个包含特征向量的数组：

(array([ 1\.        ,  0.16709381,  0.32663057]), array([[  5.77350269e-01,   7.31108409e-01,   7.90138877e-04],
 [  5.77350269e-01,  -4.65117036e-01,  -9.99813147e-01],
 [  5.77350269e-01,  -4.99145907e-01,   1.93144030e-02]]))

为特征值1选择特征向量。

目前，我们只对特征值1的特征向量感兴趣。实际上，特征值可能不完全是1，因此我们应该建立误差容限。我们可以在0.9和1.1之间找到特征值的索引，如下所示：
```
idx_vec = np.where(np.abs(eig_out[0] - 1) < 0.1)
print("Index eigenvalue 1", idx_vec)

x = eig_out[1][:,idx_vec].flatten()
```
此代码的其余输出如下：
```
Index eigenvalue 1 (array([0]),)
Steady state vector [ 0.57735027  0.57735027  0.57735027]
Check [ 0.57735027  0.57735027  0.57735027]
```

工作原理

我们获得的特征向量的值未标准化。由于我们正在处理概率，因此它们应该合计为一个。在此示例中介绍了diff()，sign()和eig()函数。它们的描述如下：

函数	描述
`diff()`	计算离散差。默认情况下是一阶。
`sign()`	返回数组元素的符号。
`eig()`	返回数组的特征值和特征向量。

另见

第 1 章，“使用 IPython”中的“安装 matplotlib”秘籍

发现幂律

为了这个秘籍目的，假设我们正在经营一家对冲基金。让它沉入；您现在是百分之一的一部分！

幂律出现在很多地方。有关更多信息，请参见这里。在这样的定律中，一个变量等于另一个变量的幂：

例如，帕累托原理是幂律。它指出财富分配不均。这个原则告诉我们，如果我们按照人们的财富进行分组，则分组的规模将成倍地变化。简而言之，富人不多，亿万富翁更少。因此是百分之一

假设在收盘价对数回报中存在幂定律。当然，这是一个很大的假设，但是幂律假设似乎到处都有。

我们不想交易太频繁，因为每笔交易涉及交易成本。假设我们希望根据重大调整（换句话说就是大幅下降）每月进行一次买卖。问题是要确定适当的信号，因为我们要在大约 20 天内每 1 天启动一次交易。

操作步骤

以下是本书代码包中powerlaw.py文件的完整代码：

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
import matplotlib.pyplot as plt

#1\. Get close prices.
today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('IBM', start, today)
close =  np.array([q[4] for q in quotes])

#2\. Get positive log returns.
logreturns = np.diff(np.log(close))
pos = logreturns[logreturns > 0]

#3\. Get frequencies of returns.
counts, rets = np.histogram(pos)
# 0 counts indices
indices0 = np.where(counts != 0)
rets = rets[:-1] + (rets[1] - rets[0])/2
# Could generate divide by 0 warning
freqs = 1.0/counts
freqs = np.take(freqs, indices0)[0]
rets = np.take(rets, indices0)[0]
freqs =  np.log(freqs)

#4\. Fit the frequencies and returns to a line.
p = np.polyfit(rets,freqs, 1)

#5\. Plot the results.
plt.title('Power Law')
plt.plot(rets, freqs, 'o', label='Data')
plt.plot(rets, p[0] * rets + p[1], label='Fit')
plt.xlabel('Log Returns')
plt.ylabel('Log Frequencies')
plt.legend()
plt.grid()
plt.show()

首先，让我们从 Yahoo 金融获取过去一年的历史日末数据。之后，我们提取该时段的收盘价。在上一秘籍中描述了这些步骤：

获得正的对数回报。

现在，计算收盘价的对数回报。有关对数回报中的更多信息，请参考这里。

首先，我们将获取收盘价的对数，然后使用diff() NumPy 函数计算这些值的第一个差异。让我们从对数回报中选择正值：
```
logreturns = np.diff(np.log(close))
pos = logreturns[logreturns > 0]
```

获得回报的频率。

我们需要使用histogram()函数获得回报的频率。返回计数和垃圾箱数组。最后，我们需要记录频率，以获得良好的线性关系：

counts, rets = np.histogram(pos)
# 0 counts indices
indices0 = np.where(counts != 0)
rets = rets[:-1] + (rets[1] - rets[0])/2
# Could generate divide by 0 warning
freqs = 1.0/counts
freqs = np.take(freqs, indices0)[0]
rets = np.take(rets, indices0)[0]
freqs =  np.log(freqs)

拟合频率并返回一条线。

使用polyfit()函数进行线性拟合：
```
p = np.polyfit(rets,freqs, 1)
```

绘制结果。

最后，我们将绘制数据并将其与 matplotlib 线性拟合：

plt.title('Power Law')
plt.plot(rets, freqs, 'o', label='Data')
plt.plot(rets, p[0] * rets + p[1], label='Fit')
plt.xlabel('Log Returns')
plt.ylabel('Log Frequencies')
plt.legend()
plt.grid()
plt.show()

我们得到一个很好的线性拟合，收益率和频率图，如下所示：

工作原理

histogram()函数计算数据集的直方图。它返回直方图值和桶的边界。 polyfit()函数将数据拟合给定阶数的多项式。在这种情况下，我们选择了线性拟合。我们发现了幂律法-您必须谨慎地提出此类主张，但证据看起来很有希望。

另见

第 1 章，“使用 IPython”中的“安装 matplotlib”秘籍
histogram()函数的文档页面
polyfit()函数的文档页面

逢低定期交易

股票价格周期性地下跌和上涨。我们将研究股价对数收益的概率分布，并尝试一个非常简单的策略。该策略基于对均值的回归。这是弗朗西斯·高尔顿爵士最初在遗传学中发现的一个概念。据发现，高大父母的孩子往往比父母矮。矮小父母的孩子往往比父母高。当然，这是一种统计现象，没有考虑基本因素和趋势，例如营养改善。均值回归也与股票市场有关。但是，它不提供任何保证。如果公司开始生产不良产品或进行不良投资，则对均值的回归将无法节省股票。

让我们首先下载股票的历史数据，例如AAPL。接下来，我们计算收盘价的每日对数回报率。我们将跳过这些步骤，因为它们在上一个秘籍中已经完成。

准备

如有必要，安装 matplotlib 和 SciPy。有关相应的秘籍，请参见“另请参见”部分。

操作步骤

以下是本书代码包中periodic.py文件的完整代码：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt

#1\. Get close prices.
today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('AAPL', start, today)
close =  np.array([q[4] for q in quotes])

#2\. Get log returns.
logreturns = np.diff(np.log(close))

#3\. Calculate breakout and pullback
freq = 0.02
breakout = scipy.stats.scoreatpercentile(logreturns, 100 * (1 - freq) )
pullback = scipy.stats.scoreatpercentile(logreturns, 100 * freq)

#4\. Generate buys and sells
buys = np.compress(logreturns < pullback, close)
sells = np.compress(logreturns > breakout, close)
print(buys)
print(sells)
print(len(buys), len(sells))
print(sells.sum() - buys.sum())

#5\. Plot a histogram of the log returns
plt.title('Periodic Trading')
plt.hist(logreturns)
plt.grid()
plt.xlabel('Log Returns')
plt.ylabel('Counts')
plt.show()

现在来了有趣的部分：

计算突破和回调。

假设我们要每年进行五次交易，大约每 50 天进行一次。一种策略是在价格下跌一定百分比时进行买入（回调），而在价格上涨另一百分比时进行卖出（突破）。

通过设置适合我们交易频率的百分比，我们可以匹配相应的对数回报。 SciPy 提供scoreatpercentile()函数，我们将使用：
```
freq = 0.02
breakout = scipy.stats.scoreatpercentile(logreturns, 100 * (1 - freq) )
pullback = scipy.stats.scoreatpercentile(logreturns, 100 * freq)
```
产生买卖。

使用compress() NumPy 函数为我们的收盘价数据生成买卖。该函数根据条件返回元素：
```
buys = np.compress(logreturns < pullback, close)
sells = np.compress(logreturns > breakout, close)
print(buys)
print(sells)
print(len(buys), len(sells))
print(sells.sum() - buys.sum())
```
AAPL和 50 天期间的输出如下：
```
[  77.76375466   76.69249773  102.72        101.2          98.57      ]
[ 74.95502967  76.55980292  74.13759123  80.93512599  98.22      ]
5 5
-52.1387025726
```
因此，如果我们买卖AAPL股票五次，我们将损失 52 美元。当我运行脚本时，经过更正后整个市场都处于恢复模式。您可能不仅要查看AAPL的股价，还可能要查看APL和SPY的比率。 SPY可以用作美国股票市场的代理。
绘制对数回报的直方图。

只是为了好玩，让我们用 matplotlib 绘制对数回报的直方图：
```
plt.title('Periodic Trading')
plt.hist(logreturns)
plt.grid()
plt.xlabel('Log Returns')
plt.ylabel('Counts')
plt.show()
```
直方图如下所示：

工作原理

我们遇到了compress()函数，该函数返回一个数组，其中包含满足给定条件的输入的数组元素。输入数组保持不变。

另见

第 1 章，“使用 IPython”中的“安装 matplotlib”秘籍
第 2 章，“高级索引和数组概念”中的“安装 SciPy”秘籍
本章中的“发现幂律”秘籍
compress()函数文档页面

随机模拟交易

在先前的秘籍中，我们尝试了一种交易想法。但是，我们没有基准可以告诉我们所获得的结果是否良好。在这种情况下，通常以我们应该能够击败随机过程为前提进行随机交易。我们将从交易年度中随机抽出几天来模拟交易。这应该说明使用 NumPy 处理随机数。

准备

如有必要，安装 matplotlib。请参考相应秘籍的“另请参见”部分。

操作步骤

以下是本书代码包中random_periodic.py文件的完整代码：

from __future__ import print_function
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
import matplotlib.pyplot as plt

def get_indices(high, size):
   #2\. Generate random indices
   return np.random.randint(0, high, size)

#1\. Get close prices.
today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('AAPL', start, today)
close =  np.array([q[4] for q in quotes])

nbuys = 5
N = 2000
profits = np.zeros(N)

for i in xrange(N):
   #3\. Simulate trades
   buys = np.take(close, get_indices(len(close), nbuys))
   sells = np.take(close, get_indices(len(close), nbuys))
   profits[i] = sells.sum() - buys.sum()

print("Mean", profits.mean())
print("Std", profits.std())

#4\. Plot a histogram of the profits
plt.title('Simulation')
plt.hist(profits)
plt.xlabel('Profits')
plt.ylabel('Counts')
plt.grid()
plt.show()

首先，我们需要一个数组，其中填充了随机整数：

生成随机索引。

您可以使用randint() NumPy 函数生成随机整数。这将与一个交易年度的随机日期相关联：
```
return np.random.randint(0, high, size)
```

模拟交易。

您可以使用上一步中的随机指数来模拟交易。使用take() NumPy 函数从步骤 1 的数组中提取随机收盘价：

buys = np.take(close, get_indices(len(close), nbuys))
sells = np.take(close, get_indices(len(close), nbuys))
profits[i] = sells.sum() - buys.sum()

绘制大量模拟的利润直方图：
```
plt.title('Simulation')
plt.hist(profits)
plt.xlabel('Profits')
plt.ylabel('Counts')
plt.grid()
plt.show()
```
以下是AAPL的 2,000 个模拟结果的直方图的屏幕截图，一年内进行了五次买卖：

工作原理

我们使用了randint()函数，该函数可以在numpy.random模块中找到。该模块包含更方便的随机生成器，如下表所述：

函数	描述
`rand()`	从`[0,1]`上的均匀分布中创建一个数组，其形状基于大小参数。如果未指定大小，则返回单个浮点数。
`randn()`	从均值`0`和方差`1`的正态分布中采样值。大小参数的作用与`rand()`相同。
`randint()`	返回一个给定下限，可选上限和可选输出形状的整数数组。

另见

第 1 章，“使用 IPython”中的“安装 matplotlib”秘籍

用 Eratosthenes 筛子筛选质数

Eratosthenes 筛子是一种过滤质数的算法。迭代地标识找到的质数的倍数。根据定义，倍数不是质数，可以消除。此筛子对于不到 1000 万的质数有效。现在让我们尝试找到第 10001 个质数。

操作步骤

第一步是创建自然数列表：

创建一个连续整数列表。 NumPy 为此具有arange()函数：
```
a = np.arange(i, i + LIM, 2)
```

筛选出p的倍数。

我们不确定这是否是 Eratosthenes 想要我们做的，但是它有效。在下面的代码中，我们传递 NumPy 数组，并去除除以p时余数为零的所有元素：

a = a[a % p != 0]

以下是此问题的完整代码：

from __future__ import print_function
import numpy as np

LIM = 10 ** 6
N = 10 ** 9
P = 10001
primes = []
p = 2

#By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

#What is the 10 001st prime number?

def sieve_primes(a, p):
   #2\. Sieve out multiples of p
   a = a[a % p != 0]

   return a

for i in xrange(3, N, LIM):
   #1\. Create a list of consecutive integers
   a = np.arange(i, i + LIM, 2)

   while len(primes) < P:
      a = sieve_primes(a, p)
      primes.append(p)

      p = a[0]

print(len(primes), primes[P-1])

四、将 NumPy 与世界的其他地方连接

在本章中，我们将介绍以下秘籍：

使用缓冲区协议
使用数组接口
与 MATLAB 和 Octave 交换数据
安装 RPy2
与 R 交互
安装 JPype
将 NumPy 数组发送到 JPype
安装 Google App Engine
在 Google Cloud 上部署 NumPy 代码
在 PythonAnywhere Web 控制台中运行 NumPy 代码

简介

本章是关于互操作性的。我们必须不断提醒自己，NumPy 在科学（Python）软件生态系统中并不孤单。与 SciPy 和 matplotlib 一起工作非常容易。还存在用于与其他 Python 包互操作性的协议。在 Python 生态系统之外，Java，R，C 和 Fortran 等语言非常流行。我们将详细介绍与这些环境交换数据的细节。

此外，我们还将讨论如何在云上获取 NumPy 代码。这是在快速移动的空间中不断发展的技术。您可以使用许多选项，其中包括 Google App Engine 和 PythonAnywhere。

使用缓冲区协议

基于 C 的 Python 对象具有所谓的缓冲区接口。 Python 对象可以公开其数据以进行直接访问，而无需复制它们。缓冲区协议使我们能够与其他 Python 软件进行通信，例如 Python 图像库（PIL）。

我们将看到一个从 NumPy 数组保存 PIL 图像的示例。

准备

如有必要，请安装 PIL 和 SciPy。有关说明，查阅本秘籍的“另见”部分。

操作步骤

该秘籍的完整代码在本书代码包的buffer.py文件中：

import numpy as np
import Image #from PIL import Image (Python3)
import scipy.misc

lena = scipy.misc.lena()
data = np.zeros((lena.shape[0], lena.shape[1], 4), dtype=np.int8)
data[:,:,3] = lena.copy()
img = Image.frombuffer("RGBA", lena.shape, data, 'raw', "RGBA", 0, 1)
img.save('lena_frombuffer.png')

data[:,:,3] = 255 
data[:,:,0] = 222 
img.save('lena_modified.png')

首先，我们需要一个 NumPy 数组来玩：

在前面的章节中，我们看到了如何加载 Lena 的样例图像。创建一个填充零的数组，并使用图像数据填充 alpha 通道：
```
lena = scipy.misc.lena()
data = np.zeros((lena.shape[0], lena.shape[1], 4), dtype=numpy.int8)
data[:,:,3] = lena.copy()
```

使用 PIL API 将数据另存为 RGBA 图像：

img = Image.frombuffer("RGBA", lena.shape, data, 'raw', "RGBA", 0, 1)
img.save('lena_frombuffer.png')

通过去除图像数据并使图像变为红色来修改数据数组。使用 PIL API 保存图像：
```
data[:,:,3] = 255 
data[:,:,0] = 222 
img.save('lena_modified.png')
```
以下是之前的图片：

注意

在计算机图形中，原点的位置与您从高中数学中知道的通常的直角坐标系不同。原点位于屏幕，画布或图像的左上角，y 轴向下。

PIL 图像对象的数据由于缓冲接口的作用而发生了变化，因此，我们看到以下图像：

工作原理

我们从缓冲区（一个 NumPy 数组）创建了一个 PIL 图像。更改缓冲区后，我们看到更改反映在图像对象中。我们这样做时没有复制 PIL 图像对象；相反，我们直接访问并修改了其数据，以使模型的图片显示红色图像。通过一些简单的更改，代码就可以与其他基于 PIL 的库一起使用，例如 Pillow。

另见

第 2 章，“高级索引和数组概念”中的“安装 PIL”
第 2 章，“高级索引和数组概念”中的“安装 SciPy”
这个页面中介绍了 Python 缓冲区协议。

使用数组接口

数组接口是用于与其他 Python 应用通信的另一种机制。顾名思义，该协议仅适用于类似数组的对象。进行了示范。让我们再次使用 PIL，但不保存文件。

准备

我们将重用先前秘籍中的部分代码，因此前提条件是相似的。在这里，我们将跳过上一秘籍的第一步，并假定它已经为人所知。

操作步骤

该秘籍的代码在本书代码包的array_interface.py文件中：

from __future__ import print_function
import numpy as np
import Image
import scipy.misc

lena = scipy.misc.lena()
data = np.zeros((lena.shape[0], lena.shape[1], 4), dtype=np.int8)
data[:,:,3] = lena.copy()
img = Image.frombuffer("RGBA", lena.shape, data, 'raw', "RGBA", 0, 1)
array_interface = img.__array_interface__
print("Keys", array_interface.keys())
print("Shape", array_interface['shape'])
print("Typestr", array_interface['typestr'])

numpy_array = np.asarray(img)
print("Shape", numpy_array.shape)
print("Data type", numpy_array.dtype)

以下步骤将使我们能够探索数组接口：

PIL Image对象具有__array_interface__属性。让我们检查它的内容。此属性的值是 Python 字典：

array_interface = img.__array_interface__
print("Keys", array_interface.keys()) 
print("Shape", array_interface['shape'])
print("Typestr", array_interface['typestr'])

此代码显示以下信息：

Keys ['shape', 'data', 'typestr']
Shape (512, 512, 4)
Typestr |u1

NumPy ndarray类也具有__array_interface__属性。我们可以使用asarray()函数将 PIL 图像转换成 NumPy 数组：
```
numpy_array = np.asarray(img)
print("Shape", numpy_array.shape)
print("Data type", numpy_array.dtype)
```
数组的形状和数据类型如下：
```
Shape (512, 512, 4)
Data type uint8
```

如您所见，形状没有改变。

工作原理

数组接口或协议使我们可以在类似数组的 Python 对象之间共享数据。 NumPy 和 PIL 都提供了这样的接口。

另见

本章中的“使用缓冲区协议”
数组接口在这个页面中进行了详细描述。

与 MATLAB 和 Octave 交换数据

MATLAB 及其开放源代码 Octave 是流行的数学应用。 scipy.io包具有savemat()函数，该函数允许您将 NumPy 数组存储为.mat文件作为 Python 字典的值。

准备

安装 MATLAB 或 Octave 超出了本书的范围。 Octave 网站上有一些安装的指南。如有必要，检查本秘籍的“另见”部分，来获取安装 SciPy 的说明。

操作步骤

该秘籍的完整代码在本书代码包的octave.py文件中：

import numpy as np
import scipy.io

a = np.arange(7)

scipy.io.savemat("a.mat", {"array": a})

一旦安装了 MATLAB 或 Octave，就需要按照以下步骤存储 NumPy 数组：

创建一个 NumPy 数组，然后调用savemat()将其存储在.mat文件中。此函数有两个参数-文件名和包含变量名和值的字典。
```
a = np.arange(7)
scipy.io.savemat("a.mat", {"array": a})
```

导航到创建文件的目录。加载文件并检查数组：

octave-3.4.0:2> load a.mat
octave-3.4.0:3> array
array =

 0
 1
 2
 3
 4
 5
 6

另见

第 2 章，“高级索引和数组概念”中的“安装 SciPy”
savemat()函数的 SciPy 文档

安装 RPy2

R 是一种流行的脚本语言，用于统计和数据分析。 RPy2 是 R 和 Python 之间的接口。我们将在此秘籍中安装 RPy2。

操作步骤

如果要安装 RPy2，请选择以下选项之一：

使用pip或easy_install进行安装：RPy2 在 PYPI 上可用，因此我们可以使用以下命令进行安装：
```
$ easy_install rpy2
```
另外，我们可以使用以下命令：
```
$ sudo pip install rpy2
$ pip freeze|grep rpy2
rpy2==2.4.2
```

从源代码安装：我们可以从tar.gz源安装 RPy2：

$ tar -xzf <rpy2_package>.tar.gz
$ cd <rpy2_package>
$ python setup.py build install

另见

与 R 交互

RPy2 只能用作从 Python 调用 R，而不能相反。我们将导入一些样本 R 数据集并绘制其中之一的数据。

准备

如有必要，请安装 RPy2。请参阅先前的秘籍。

操作步骤

该秘籍的完整代码在本书代码包的rdatasets.py文件中：

from rpy2.robjects.packages import importr
import numpy as np
import matplotlib.pyplot as plt

datasets = importr('datasets')
mtcars = datasets.__rdata__.fetch('mtcars')['mtcars']

plt.title('R mtcars dataset')
plt.xlabel('wt')
plt.ylabel('mpg')
plt.plot(mtcars)
plt.grid(True)
plt.show()

motorcars数据集在这个页面中进行了描述。让我们从加载此样本 R 数据集开始：

使用 RPy2 importr()函数将数据集加载到数组中。此函数可以导入R包。在此示例中，我们将导入数据集 R 包。从mtcars数据集创建一个 NumPy 数组：
```
datasets = importr('datasets')
mtcars = np.array(datasets.mtcars)
```
使用 matplotlib 绘制数据集：
```
plt.plot(mtcars)
plt.show()
```
数据包含英里每加仑（mpg）和重量（wt）值，单位为千分之一磅。以下屏幕快照显示了数据，它是一个二维数组：

另见

第 1 章“使用 IPython”中的“安装 matplotlib”

安装 JPype

Jython 是用于 Python 和 Java 的默认互操作性解决方案。但是，Jython 在 Java 虚拟机（JVM）上运行。因此，它无法访问主要用 C 语言编写的 NumPy 模块。 JPype 是一个开放源代码项目，试图解决此问题。接口发生在 Python 和 JVM 之间的本机级别上。让我们安装 JPype。

操作步骤

从这里下载 JPype。
打开压缩包，然后运行以下命令：
```
$ python setup.py install
```

将 NumPy 数组发送到 JPype

在此秘籍中，我们将启动 JVM 并向其发送 NumPy 数组。我们将使用标准 Java 调用打印接收到的数组。显然，您将需要安装 Java。

操作步骤

该秘籍的完整代码在本书代码包的hellojpype.py文件中：

import jpype
import numpy as np

#1\. Start the JVM
jpype.startJVM(jpype.getDefaultJVMPath())

#2\. Print hello world
jpype.java.lang.System.out.println("hello world")

#3\. Send a NumPy array
values = np.arange(7)
java_array = jpype.JArray(jpype.JDouble, 1)(values.tolist())

for item in java_array:
   jpype.java.lang.System.out.println(item)

#4\. Shutdown the JVM
jpype.shutdownJVM()

首先，我们需要从 JPype 启动 JVM：

从 JPype 启动 JVM； JPype 可以方便地找到默认的 JVM 路径：
```
jpype.startJVM(jpype.getDefaultJVMPath())
```
仅出于传统原因，让我们打印"hello world"：
```
jpype.java.lang.System.out.println("hello world")
```

创建一个 NumPy 数组，将其转换为 Python 列表，然后将其传递给 JPype。现在很容易打印数组元素：

values = np.arange(7)
java_array = jpype.JArray(jpype.JDouble, 1)(values.tolist())

for item in java_array:
    jpype.java.lang.System.out.println(item)

完成后，让我们关闭 JVM：
```
jpype.shutdownJVM()
```
JPype 中一次只能运行一个 JVM。如果我们忘记关闭 JVM，则可能导致意外错误。程序输出如下：
```
hello world
0.0
1.0
2.0
3.0
4.0
5.0
6.0
JVM activity report     :
 classes loaded       : 31
JVM has been shutdown
```

工作原理

JPype 允许我们启动和关闭 JVM。它为标准 Java API 调用提供了包装器。如本例所示，我们可以传递要由 JArray 包装器转换为 Java 数组的 Python 列表。 JPype 使用 Java 本机接口（JNI），这是本机 C 代码和 Java 之间的桥梁。不幸的是，使用 JNI 会损害性能，因此您必须注意这一事实。

另见

本章中的“安装 JPype”
JPype 主页

安装 Google App Engine

Google App Engine（GAE）使您可以在 Google Cloud 上构建 Web 应用。自 2012 年以来，是 NumPy 的官方支持；您需要一个 Google 帐户才能使用 GAE。

操作步骤

第一步是下载 GAE：

从这里下载适用于您的操作系统的 GAE。

您也可以从此页面下载文档和 GAE Eclipse 插件。 如果使用 Eclipse 开发，则一定要安装它。

开发环境。

GAE 带有一个模拟生产云的开发环境。在撰写本书时，GAE 正式仅支持 Python 2.5 和 2.7。 GAE 将尝试在您的系统上找到 Python；但是，例如，如果您有多个 Python 版本，则可能需要自行设置。您可以在启动器应用的首选项对话框中设置此设置。

SDK 中有两个重要的脚本：
- dev_appserver.py：开发服务器
- appcfg.py：部署在云上
在 Windows 和 Mac 上，有一个 GAE 启动器应用。启动器具有运行和部署按钮，它们执行与上述脚本相同的操作。

在 Google Cloud 上部署 NumPy 代码

部署 GAE 应用非常容易。对于 NumPy，需要额外的配置步骤，但这仅需几分钟。

操作步骤

让我们创建一个新的应用：

使用启动器创建一个新应用（文件 | 新应用）。命名为numpycloud。这将创建一个包含以下文件的同名文件夹：
- app.yaml：YAML 应用配置文件
- favicon.ico：一个图标
- index.yaml：自动生成的文件
- main.py：Web 应用的主要入口点

将 NumPy 添加到库中。

首先，我们需要让 GAE 知道我们要使用 NumPy。将以下行添加到库部分中的app.yaml配置文件中：

- name: NumPy
  version: "1.6.1"

这不是最新的 NumPy 版本，但它是 GAE 当前支持的最新版本。配置文件应具有以下内容：

application: numpycloud
version: 1
runtime: python27
api_version: 1
threadsafe: yes

handlers:
- url: /favicon\.ico
  static_files: favicon.ico
  upload: favicon\.ico

- url: .*
  script: main.app

libraries:
- name: webapp2
  version: "2.5.1"
- name: numpy
  version: "1.6.1"

为了演示我们可以使用 NumPy 代码，让我们修改main.py文件。有一个MainHandler类，带有用于 GET 请求的处理器方法。用以下代码替换此方法：

def get(self):
    self.response.out.write('Hello world!<br/>')
    self.response.out.write('NumPy sum = ' + str(numpy.arange(7).sum()))

最后，我们将提供以下代码：

import webapp2
import numpy

class MainHandler(webapp2.RequestHandler):
    def get(self):
      self.response.out.write('Hello world!<br/>')
      self.response.out.write('NumPy sum = ' + str(numpy.arange(7).sum()))

app = webapp2.WSGIApplication([('/', MainHandler)],
                              debug=True)

如果您单击在 GAE 启动器中浏览按钮（在 Linux 上，以项目根为参数运行dev_appserver.py），则您应该在默认浏览器中看到一个包含以下文字的网页：

Hello world!NumPy sum = 21

工作原理

GAE 是免费的，具体取决于使用了多少资源。您最多可以创建 10 个 Web 应用。 GAE 采用沙盒方法，这意味着 NumPy 暂时无法使用，但现在可以使用，如本秘籍所示。

在 PythonAnywhere Web 控制台中运行 NumPy 代码

在第 1 章，“使用 IPython”中，我们已经看到了运行 PythonAnywhere 控制台的过程，而没有任何权限。此秘籍将需要您有一个帐户，但不要担心-它是免费的，如果您不需要太多资源，至少是免费的。

注册是一个非常简单的过程，此处将不涉及。 NumPy 已经与其他 Python 软件一起安装。有关完整列表，请参见这里。

我们将建立一个简单的脚本，该脚本每分钟从 Google 财经获取价格数据，并使用 NumPy 对价格进行简单的统计。

操作步骤

当我们签名后，我们可以登录并查看 PythonAnywhere 信息中心。

编写代码。此示例的完整代码如下：

from __future__ import print_function
import urllib2
import re
import time
import numpy as np

prices = np.array([])

for i in xrange(3):
   req = urllib2.Request('http://finance.google.com/finance/info?client=ig&q=AAPL')
   req.add_header('User-agent', 'Mozilla/5.0')
   response = urllib2.urlopen(req)
   page = response.read()
   m = re.search('l_cur" : "(.*)"', page)
   prices = np.append(prices, float(m.group(1)))
   avg = prices.mean()
   stddev = prices.std()

   devFactor = 1
   bottom = avg - devFactor * stddev
   top = avg + devFactor * stddev
   timestr = time.strftime("%H:%M:%S", time.gmtime())

   print(timestr, "Average", avg, "-Std", bottom, "+Std", top)
   time.sleep(60)

除我们在其中生长包含价格的 NumPy 数组并计算价格的均值和标准差的位以外，大多数都是标准 Python。如果有股票代号，例如AAPL，则可以使用 URL 从 Google 财经下载 JSON 格式的价格数据。该 URL 当然可以更改。

接下来，我们使用正则表达式解析 JSON 以提取价格。此价格已添加到 NumPy 数组中。我们计算价格的均值和标准差。价格是根据标准差乘以我们指定的某个因素后在时间戳的顶部和底部打印出来的。

上传代码。

在本地计算机上完成代码后，我们可以将脚本上传到 PythonAnywhere。转到仪表板，然后单击文件选项卡。从页面底部的小部件上传脚本。
要运行代码，请单击控制台选项卡，然后单击 Bash 链接。 PythonAnywhere 应该立即为我们创建一个 bash 控制台。现在，我们可以在一个标准差范围内运行AAPL程序，如以下屏幕截图所示：

工作原理

如果您想在远程服务器上运行 NumPy 代码，则 PythonAnywhere 是完美的选择，尤其是当您需要程序在计划的时间执行时。至少对于免费帐户而言，进行交互式工作并不那么方便，因为每当您在 Web 控制台中输入文本时都会有一定的滞后。

但是，正如我们所看到的，可以在本地创建和测试程序，并将其上传到 PythonAnywhere。这也会释放本地计算机上的资源。我们可以做一些花哨的事情，例如根据股价发送电子邮件或安排在交易时间内激活脚本。通过，使用 Google App Engine 也可以做到这一点，但是它是通过 Google 方式完成的，因此您需要了解其 API。

五、音频和图像处理

在本章中，我们将介绍 NumPy 和 SciPy 的基本图像和音频（WAV 文件）处理。在以下秘籍中，我们将使用 NumPy 对声音和图像进行有趣的操作：

将图像加载到内存映射中
添加图像
图像模糊
重复音频片段
产生声音
设计音频过滤器
使用 Sobel 过滤器进行边界检测

简介

尽管本书中的所有章节都很有趣，但在本章中，我们确实会继续努力并专注于获得乐趣。在第 10 章，“Scikits 的乐趣”中，您会发现更多使用scikits-image的图像处理秘籍。不幸的是，本书没有对音频文件的直接支持，因此您确实需要运行代码示例以充分了解其中的秘籍。

将图像加载到内存映射中

建议将大文件加载到内存映射中。内存映射文件仅加载大文件的一小部分。 NumPy 内存映射类似于数组。在此示例中，我们将生成彩色正方形的图像并将其加载到内存映射中。

准备

如有必要，“安装 matplotlib”的“另请参见”部分具有对相应秘籍的引用。

操作步骤

我们将通过初始化数组来开始：

首先，我们需要初始化以下数组：
- 保存图像数据的数组
- 具有正方形中心随机坐标的数组
- 具有平方的随机半径（复数个半径）的数组
- 具有正方形随机颜色的数组
初始化数组：
```
img = np.zeros((N, N), np.uint8)
NSQUARES = 30
centers = np.random.random_integers(0, N, size=(NSQUARES, 2))
radii = np.random.randint(0, N/9, size=NSQUARES)
colors = np.random.randint(100, 255, size=NSQUARES)
```
如您所见，我们正在将第一个数组初始化为零。其他数组使用numpy.random包中的函数初始化，这些函数生成随机整数。
下一步是生成正方形。我们在上一步中使用数组创建正方形。使用clip()函数，我们将确保正方形不会在图像区域外徘徊。 meshgrid()函数为我们提供了正方形的坐标。如果我们给此函数两个大小分别为N和M的数组，它将给我们两个形状为N x M的数组。第一个数组的元素将沿 x 轴重复。第二个数组将沿 y 轴重复其元素。以下示例 IPython 会话应该使这一点更加清楚：

注意
```
In: x = linspace(1, 3, 3)

In: x
Out: array([ 1.,  2.,  3.])
In: y = linspace(1, 2, 2)

In: y
Out: array([ 1.,  2.])

In: meshgrid(x, y)
Out:
[array([[ 1.,  2.,  3.],
 [ 1.,  2.,  3.]]),
 array([[ 1.,  1.,  1.],
 [ 2.,  2.,  2.]])]
```

最后，我们将设置正方形的颜色：

for i in xrange(NSQUARES):
   xindices = range(centers[i][0] - radii[i], centers[i][0] + radii[i])
   xindices = np.clip(xindices, 0, N - 1)
   yindices = range(centers[i][1] - radii[i], centers[i][1] + radii[i])
   yindices = np.clip(yindices, 0, N - 1)

   if len(xindices) == 0 or len(yindices) == 0:
      continue
   coordinates = np.meshgrid(xindices, yindices)
   img[coordinates] = colors[i]

在将图像数据加载到内存映射之前，我们需要使用tofile()函数将其存储在文件中。然后使用memmap()函数将图像文件中的图像数据加载到内存映射中：
```
img.tofile('random_squares.raw')
img_memmap = np.memmap('random_squares.raw', shape=img.shape)
```

为了确认一切正常，我们使用 matplotlib 显示图像：

plt.imshow(img_memmap)
plt.axis('off')
plt.show()

注意，我们没有显示轴。生成图像的示例如下所示：

这是本书代码包中memmap.py文件的完整源代码：

import numpy as np
import matplotlib.pyplot as plt

N = 512
NSQUARES = 30

# Initialize
img = np.zeros((N, N), np.uint8)
centers = np.random.random_integers(0, N, size=(NSQUARES, 2))
radii = np.random.randint(0, N/9, size=NSQUARES)
colors = np.random.randint(100, 255, size=NSQUARES)

# Generate squares
for i in xrange(NSQUARES):
   xindices = range(centers[i][0] - radii[i], centers[i][0] + radii[i])
   xindices = np.clip(xindices, 0, N - 1)
   yindices = range(centers[i][1] - radii[i], centers[i][1] + radii[i])
   yindices = np.clip(yindices, 0, N - 1)

   if len(xindices) == 0 or len(yindices) == 0:
      continue

   coordinates = np.meshgrid(xindices, yindices)
   img[coordinates] = colors[i]

# Load into memory map
img.tofile('random_squares.raw')
img_memmap = np.memmap('random_squares.raw', shape=img.shape)

# Display image
plt.imshow(img_memmap)
plt.axis('off')
plt.show()

工作原理

我们在此秘籍中使用了以下函数：

函数	描述
`zeros()`	此函数给出一个由零填充的数组。
`random_integers()`	此函数返回一个数组，数组中的随机整数值在上限和下限之间。
`randint()`	该函数与`random_integers()`函数相同，除了它使用半开间隔而不是关闭间隔。
`clip()`	该函数在给定最小值和最大值的情况下裁剪数组的值。
`meshgrid()`	此函数从包含 x 坐标的数组和包含 y 坐标的数组返回坐标数组。
`tofile()`	此函数将数组写入文件。
`memmap()`	给定文件名，此函数从文件创建 NumPy 内存映射。（可选）您可以指定数组的形状。
`axis()`	该函数是用于配置绘图轴的 matplotlib 函数。例如，我们可以将其关闭。

另见

第 1 章“使用 IPython”中的“安装 matplotlib”
NumPy 内存映射文档

合成图像

在此秘籍中，我们将结合著名的 Mandelbrot 分形和 Lena 图像。 Mandelbrot 集是数学家 Benoit Mandelbrot 发明的。这些类型的分形由递归公式定义，您可以在其中通过将当前拥有的复数乘以自身并为其添加一个常数来计算序列中的下一个复数。本秘籍将涵盖更多详细信息。

准备

如有必要，请安装 SciPy。“另请参见”部分具有相关秘籍的参考。

操作步骤

首先初始化数组，然后生成和绘制分形，最后将分形与 Lena 图像组合：

使用meshgrid()，zeros()和linspace()函数初始化对应于图像区域中像素的x，y和z数组：

x, y = np.meshgrid(np.linspace(x_min, x_max, SIZE),
                   np.linspace(y_min, y_max, SIZE))
c = x + 1j * y
z = c.copy()
fractal = np.zeros(z.shape, dtype=np.uint8) + MAX_COLOR

如果z是复数，则对于 Mandelbrot 分形具有以下关系：

在此，c是常数复数。这可以在复平面上绘制，水平轴显示实数值，垂直轴显示虚数值。我们将使用所谓的逃逸时间算法绘制分形。该算法以大约 2 个单位的距离扫描原点周围小区域中的点。这些点中的每一个都用作c值，并根据逃避区域所需的迭代次数为其指定颜色。如果所需的迭代次数超过了预定义的次数，则像素将获得默认背景色。有关更多信息，请参见此秘籍中已提及的维基百科文章：

for n in range(ITERATIONS):
    print(n)
    mask = numpy.abs(z) <= 4 
    z[mask] = z[mask] ** 2 +  c[mask]
    fractal[(fractal == MAX_COLOR) & (-mask)] = (MAX_COLOR - 1) * n / ITERATIONS
Plot the fractal with matplotlib:
plt.subplot(211)
plt.imshow(fractal)
plt.title('Mandelbrot')
plt.axis('off')
Use the choose() function to pick a value from the fractal or Lena array:
plt.subplot(212)
plt.imshow(numpy.choose(fractal < lena, [fractal, lena]))
plt.axis('off')
plt.title('Mandelbrot + Lena')

结果图像如下所示：

以下是本书代码集中mandelbrot.py文件中该秘籍的完整代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import lena

ITERATIONS = 10
lena = lena()
SIZE = lena.shape[0]
MAX_COLOR = 255.
x_min, x_max = -2.5, 1
y_min, y_max = -1, 1

# Initialize arrays
x, y = np.meshgrid(np.linspace(x_min, x_max, SIZE),
                   np.linspace(y_min, y_max, SIZE))
c = x + 1j * y
z = c.copy()
fractal = np.zeros(z.shape, dtype=np.uint8) + MAX_COLOR
# Generate fractal
for n in range(ITERATIONS):
    mask = np.abs(z) <= 4 
    z[mask] = z[mask] ** 2 +  c[mask]
    fractal[(fractal == MAX_COLOR) & (-mask)] = (MAX_COLOR - 1) * n / ITERATIONS

# Display the fractal
plt.subplot(211)
plt.imshow(fractal)
plt.title('Mandelbrot')
plt.axis('off')

# Combine with lena
plt.subplot(212)
plt.imshow(np.choose(fractal < lena, [fractal, lena]))
plt.axis('off')
plt.title('Mandelbrot + Lena')

plt.show()

工作原理

在此示例中使用了以下函数：

函数	描述
`linspace()`	此函数返回范围内具有指定间隔的数字
`choose()`	此函数通过根据条件从数组中选择值来创建数组
`meshgrid()`	此函数从包含 x 坐标的数组和包含 y 坐标的数组返回坐标数组

另见

第 1 章，“使用 IPython”中的“安装 matplotlib”秘籍
第 2 章，“高级索引和数组”中的“安装 SciPy”秘籍

图像模糊

我们可以使用高斯过滤器来模糊图像。该过滤器基于正态分布。对应的 SciPy 函数需要标准差作为参数。在此秘籍中，我们还将绘制极地玫瑰和螺旋形。这些数字没有直接关系，但是在这里将它们组合起来似乎更有趣。

操作步骤

我们从初始化极坐标图开始，之后我们将模糊 Lena 图像并使用极坐标进行绘图：

初始化极坐标图：

NFIGURES = 5
k = np.random.random_integers(1, 5, NFIGURES)
a = np.random.random_integers(1, 5, NFIGURES)
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']

要模糊 Lena，请应用标准差为4的高斯过滤器：

plt.subplot(212)
blurred = scipy.ndimage.gaussian_filter(lena, sigma=4)

plt.imshow(blurred)
plt.axis('off')

matplotlib 有一个polar()函数，它以极坐标进行绘制：

theta = np.linspace(0, k[0] * np.pi, 200)
plt.polar(theta, np.sqrt(theta), choice(colors))

for i in xrange(1, NFIGURES):
 theta = np.linspace(0, k[i] * np.pi, 200)
 plt.polar(theta, a[i] * np.cos(k[i] * theta), choice(colors))

这是这本书的代码集中spiral.py文件中该秘籍的完整代码：

import numpy as np
import matplotlib.pyplot as plt
from random import choice
import scipy
import scipy.ndimage
# Initialization
NFIGURES = 5
k = np.random.random_integers(1, 5, NFIGURES)
a = np.random.random_integers(1, 5, NFIGURES)

colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']

lena = scipy.misc.lena()
plt.subplot(211)
plt.imshow(lena)
plt.axis('off')

# Blur Lena
plt.subplot(212)
blurred = scipy.ndimage.gaussian_filter(lena, sigma=4)

plt.imshow(blurred)
plt.axis('off')

# Plot in polar coordinates
theta = np.linspace(0, k[0] * np.pi, 200)
plt.polar(theta, np.sqrt(theta), choice(colors))

for i in xrange(1, NFIGURES):
   theta = np.linspace(0, k[i] * np.pi, 200)
   plt.polar(theta, a[i] * np.cos(k[i] * theta), choice(colors))

plt.axis('off')

plt.show()

工作原理

在本教程中，我们使用了以下函数：

函数	描述
`gaussian_filter()`	此函数应用高斯过滤器
`random_integers()`	此函数返回一个数组，数组中的随机整数值在上限和下限之间
`polar()`	该函数使用极坐标绘制图形

另见

可以在这个页面中找到scipy.ndimage文档。

重复音频片段

正如我们在第 2 章，“高级索引和数组概念”中所看到的那样，我们可以使用 WAV 文件来完成整洁的事情。只需使用urllib2标准 Python 模块下载文件并将其加载到 SciPy 中即可。让我们下载一个 WAV 文件并重复 3 次。我们将跳过在第 2 章，“高级索引和数组概念”中已经看到的一些步骤。

操作步骤

尽管NumPy具有repeat()函数，但在这种情况下，更适合使用tile()函数。函数repeat()的作用是通过重复单个元素而不重复其内容来扩大数组。以下 IPython 会话应阐明这些函数之间的区别：
```
In: x = array([1, 2])

In: x
Out: array([1, 2])

In: repeat(x, 3)
Out: array([1, 1, 1, 2, 2, 2])

In: tile(x, 3)
Out: array([1, 2, 1, 2, 1, 2])
```
现在，有了这些知识，就可以应用tile()函数：
```
repeated = np.tile(data, 3)
```

使用 matplotlib 绘制音频数据：

plt.title("Repeated")
plt.plot(repeated)

原始声音数据和重复数据图如下所示：

这是本书代码包中repeat_audio.py文件中该秘籍的完整代码：

import scipy.io.wavfile
import matplotlib.pyplot as plt
import urllib2
import numpy as np

response = urllib2.urlopen('http://www.thesoundarchive.com/austinpowers/smashingbaby.wav')
print(response.info())
WAV_FILE = 'smashingbaby.wav'
filehandle = open(WAV_FILE, 'w')
filehandle.write(response.read())
filehandle.close()
sample_rate, data = scipy.io.wavfile.read(WAV_FILE)
print("Data type", data.dtype, "Shape", data.shape)

plt.subplot(2, 1, 1)
plt.title("Original")
plt.plot(data)

plt.subplot(2, 1, 2)

# Repeat the audio fragment
repeated = np.tile(data, 3)

# Plot the audio data
plt.title("Repeated")
plt.plot(repeated)
scipy.io.wavfile.write("repeated_yababy.wav",
    sample_rate, repeated)

plt.show()

工作原理

以下是此秘籍中最重要的函数：

函数	描述
`scipy.io.wavfile.read()`	将 WAV 文件读入数组
`numpy.tile()`	重复数组指定次数
`scipy.io.wavfile.write()`	从 NumPy 数组中以指定的采样率创建 WAV 文件

另见

可以在这个页面中找到 scipy.io文档。

产生声音

声音可以用具有一定幅度，频率和相位的正弦波在数学上表示如下。我们可以从这个页面中指定的列表中随机选择符合以下公式的频率：

此处，n是钢琴键的编号。我们将键的编号从 1 到 88。我们将随机选择振幅，持续时间和相位。

操作步骤

首先初始化随机值，然后生成正弦波，编写旋律，最后使用 matplotlib 绘制生成的音频数据：

初始化随机值为：

200-2000之间的幅度
0.01-0.2的持续时间
使用已经提到的公式的频率

0和2 pi之间的相位值：

NTONES = 89
amps = 2000\. * np.random.random((NTONES,)) + 200.
durations = 0.19 * np.random.random((NTONES,)) + 0.01
keys = np.random.random_integers(1, 88, NTONES)
freqs = 440.0 * 2 ** ((keys - 49.)/12.)
phi = 2 * np.pi * np.random.random((NTONES,))

编写generate()函数以生成正弦波：

def generate(freq, amp, duration, phi):
 t = np.linspace(0, duration, duration * RATE)
 data = np.sin(2 * np.pi * freq * t + phi) * amp

 return data.astype(DTYPE)

一旦我们产生了一些音调，我们只需要组成一个连贯的旋律。现在，我们将连接正弦波。这不会产生良好的旋律，但可以作为更多实验的起点：
```
for i in xrange(NTONES):
   newtone = generate(freqs[i], amp=amps[i], duration=durations[i], phi=phi[i])
   tone = np.concatenate((tone, newtone))
```

使用 matplotlib 绘制生成的音频数据：

plt.plot(np.linspace(0, len(tone)/RATE, len(tone)), tone)
plt.show()

生成的音频数据图如下：

可以在此处找到本示例的源代码，该代码来自本书代码包中的tone_generation.py文件：

import scipy.io.wavfile
import numpy as np
import matplotlib.pyplot as plt

RATE = 44100
DTYPE = np.int16

# Generate sine wave
def generate(freq, amp, duration, phi):
 t = np.linspace(0, duration, duration * RATE)
 data = np.sin(2 * np.pi * freq * t + phi) * amp

 return data.astype(DTYPE)

# Initialization
NTONES = 89
amps = 2000\. * np.random.random((NTONES,)) + 200.
durations = 0.19 * np.random.random((NTONES,)) + 0.01
keys = np.random.random_integers(1, 88, NTONES)
freqs = 440.0 * 2 ** ((keys - 49.)/12.)
phi = 2 * np.pi * np.random.random((NTONES,))

tone = np.array([], dtype=DTYPE) 

# Compose 
for i in xrange(NTONES):
   newtone = generate(freqs[i], amp=amps[i], duration=durations[i], phi=phi[i])
   tone = np.concatenate((tone, newtone))

scipy.io.wavfile.write('generated_tone.wav', RATE, tone)

# Plot audio data
plt.plot(np.linspace(0, len(tone)/RATE, len(tone)), tone)
plt.show()

工作原理

我们创建了带有随机生成声音的 WAV 文件。 concatenate()函数用于连接正弦波。

另见

concatenate()函数的文档

设计音频过滤器

我记得在模拟电子课上学习了所有类型的过滤器。然后，我们实际上构造了这些过滤器。可以想象，用软件制作过滤器要比用硬件制作过滤器容易得多。

我们将构建一个过滤器，并将其应用于要下载的音频片段。在本章之前，我们已经完成了一些步骤，因此我们将省略那些部分。

操作步骤

顾名思义，iirdesign()函数允许我们构造几种类型的模拟和数字过滤器。可以在scipy.signal模块中找到。该模块包含信号处理函数的完整列表：

使用scipy.signal模块的iirdesign()函数设计过滤器。 IIR 代表无限冲激响应；有关更多信息，请参见这里。我们将不涉及iirdesign())函数的所有细节。在这个页面中查看文档。简而言之，这些是我们将要设置的参数：
- 频率标准化为 0 到 1
- 最大损失
- 最小衰减
- 过滤器类型：
```
b,a = scipy.signal.iirdesign(wp=0.2, ws=0.1, gstop=60, gpass=1, ftype='butter')
```
此过滤器的配置对应于 Butterworth 过滤器。 Butterworth 过滤器由物理学家 Stephen Butterworth 于 1930 年首次描述。

使用scipy.signal.lfilter()函数应用过滤器。它接受上一步的值作为参数，当然也接受要过滤的数据数组：

filtered = scipy.signal.lfilter(b, a, data)

写入新的音频文件时，请确保其数据类型与原始数据数组相同：

scipy.io.wavfile.write('filtered.wav', sample_rate, filtered.astype(data.dtype))

在绘制原始数据和过滤后的数据之后，我们得到以下图：

音频过滤器的代码列出如下：

import scipy.io.wavfile
import matplotlib.pyplot as plt
import urllib2
import scipy.signal

response =urllib2.urlopen('http://www.thesoundarchive.com/austinpowers/smashingbaby.wav')
print response.info()
WAV_FILE = 'smashingbaby.wav'
filehandle = open(WAV_FILE, 'w')
filehandle.write(response.read())
filehandle.close()
sample_rate, data = scipy.io.wavfile.read(WAV_FILE)
print("Data type", data.dtype, "Shape", data.shape)

plt.subplot(2, 1, 1)
plt.title("Original")
plt.plot(data)

# Design the filter
b,a = scipy.signal.iirdesign(wp=0.2, ws=0.1, gstop=60, gpass=1, ftype='butter')

# Filter
filtered = scipy.signal.lfilter(b, a, data)

# Plot filtered data
plt.subplot(2, 1, 2)
plt.title("Filtered")
plt.plot(filtered)

scipy.io.wavfile.write('filtered.wav', sample_rate, filtered.astype(data.dtype))

plt.show()

工作原理

我们创建并应用了 Butterworth 带通过滤器。引入了以下函数来创建过滤器：

函数	描述
`scipy.signal.iirdesign()`	创建一个 IIR 数字或模拟过滤器。此函数具有广泛的参数列表，该列表在这个页面中记录。
`scipy.signal.lfilter()`	给定一个数字过滤器，对数组进行滤波。

使用 Sobel 过滤器进行边界检测

Sobel 过滤器可以用于图像中的边界检测。边界检测基于对图像强度执行离散差分。由于图像是二维的，因此渐变也有两个分量，除非我们将自身限制为一维。我们将 Sobel 过滤器应用于 Lena 的图片。

操作步骤

在本部分中，您将学习如何应用 Sobel 过滤器来检测 Lena 图像中的边界：

要在 x 方向上应用 Sobel 过滤器，请将轴参数设置为0：
```
sobelx = scipy.ndimage.sobel(lena, axis=0, mode='constant')
```
要在 y 方向上应用 Sobel 过滤器，请将轴参数设置为1：
```
sobely = scipy.ndimage.sobel(lena, axis=1, mode='constant')
```

默认的 Sobel 过滤器仅需要输入数组：

default = scipy.ndimage.sobel(lena)

是原始图像图和所得图像图，显示了使用 Sobel 过滤器进行边界检测：

完整的边界检测代码如下：

import scipy
import scipy.ndimage
import matplotlib.pyplot as plt

lena = scipy.misc.lena()

plt.subplot(221)
plt.imshow(lena)
plt.title('Original')
plt.axis('off')

# Sobel X filter
sobelx = scipy.ndimage.sobel(lena, axis=0, mode='constant')

plt.subplot(222)
plt.imshow(sobelx)
plt.title('Sobel X')
plt.axis('off')

# Sobel Y filter
sobely = scipy.ndimage.sobel(lena, axis=1, mode='constant')

plt.subplot(223)
plt.imshow(sobely)
plt.title('Sobel Y')
plt.axis('off')

# Default Sobel filter
default = scipy.ndimage.sobel(lena)

plt.subplot(224)
plt.imshow(default)
plt.title('Default Filter')
plt.axis('off')

plt.show()

工作原理

我们将 Sobel 过滤器应用于著名模型 Lena 的图片。如本例所示，我们可以指定沿哪个轴进行计算。默认设置为独立于轴。

六、特殊数组和通用函数

在本章中，我们将介绍以下秘籍：

创建通用函数
查找勾股三元组
用chararray执行字符串操作
创建一个遮罩数组
忽略负值和极值
使用recarray函数创建一个得分表

简介

本章是关于特殊数组和通用函数的。这些是您每天可能不会遇到的主题，但是它们仍然很重要，因此在此需要提及。通用函数（Ufuncs）逐个元素或标量地作用于数组。 Ufuncs 接受一组标量作为输入，并产生一组标量作为输出。通用函数通常可以映射到它们的数学对等物上，例如加法，减法，除法，乘法等。这里提到的特殊数组是基本 NumPy 数组对象的所有子类，并提供其他功能。

创建通用函数

我们可以使用frompyfunc() NumPy 函数从 Python 函数创建通用函数。

操作步骤

以下步骤可帮助我们创建通用函数：

定义一个简单的 Python 函数以使输入加倍：
```
def double(a):
    return 2 * a
```

用frompyfunc()创建通用函数。指定输入参数的数目和返回的对象数目（均等于1）：

from __future__ import print_function
import numpy as np

def double(a):
   return 2 * a

ufunc = np.frompyfunc(double, 1, 1)
print("Result", ufunc(np.arange(4)))

该代码在执行时输出以下输出：

Result [0 2 4 6]

工作原理

我们定义了一个 Python 函数，该函数会将接收到的数字加倍。实际上，我们也可以将字符串作为输入，因为这在 Python 中是合法的。我们使用frompyfunc() NumPy 函数从此 Python 函数创建了一个通用函数。通用函数是 NumPy 类，具有特殊功能，例如广播和适用于 NumPy 数组的逐元素处理。实际上，许多 NumPy 函数都是通用函数，但是都是用 C 编写的。

另见

frompyfunc() NumPy 函数的文档

查找勾股三元组

对于本教程，您可能需要阅读有关勾股三元组的维基百科页面。勾股三元组是一组三个自然数，即a < b < c，为此，。

这是勾股三元组的示例：。

勾股三元组与勾股定理密切相关，您可能在中学几何学过的。

勾股三元组代表直角三角形的三个边，因此遵循勾股定理。让我们找到一个分量总数为 1,000 的勾股三元组。我们将使用欧几里得公式进行此操作：

在此示例中，我们将看到通用函数的运行。

操作步骤

欧几里得公式定义了m和n索引。

创建包含以下索引的数组：
```
m = np.arange(33)
n = np.arange(33)
```
第二步是使用欧几里得公式计算勾股三元组的数量a，b和c。使用outer()函数获得笛卡尔积，差和和：
```
a = np.subtract.outer(m ** 2, n ** 2)
b = 2 * np.multiply.outer(m, n)
c = np.add.outer(m ** 2, n ** 2)
```
现在，我们有许多包含a，b和c值的数组。但是，我们仍然需要找到符合问题条件的值。使用where() NumPy 函数查找这些值的索引：
```
idx =  np.where((a + b + c) == 1000)
```

使用numpy.testing模块检查解决方案：

np.testing.assert_equal(a[idx]**2 + b[idx]**2, c[idx]**2)

以下代码来自本书代码包中的triplets.py文件：

from __future__ import print_function
import numpy as np

#A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,
#a ** 2 + b ** 2 = c ** 2
#
#For example, 3 ** 2 + 4 ** 2 = 9 + 16 = 25 = 5 ** 2.
#
#There exists exactly one Pythagorean triplet for which a + b + c = 1000.
#Find the product abc.

#1\. Create m and n arrays
m = np.arange(33)
n = np.arange(33)

#2\. Calculate a, b and c
a = np.subtract.outer(m ** 2, n ** 2)
b = 2 * np.multiply.outer(m, n)
c = np.add.outer(m ** 2, n ** 2)

#3\. Find the index
idx =  np.where((a + b + c) == 1000)

#4\. Check solution
np.testing.assert_equal(a[idx]**2 + b[idx]**2, c[idx]**2)
print(a[idx], b[idx], c[idx])
      # [375] [200] [425]

工作原理

通用函数不是实函数，而是表示函数的对象。工具具有outer()方法，我们已经在实践中看到它。 NumPy 的许多标准通用函数都是用 C 实现的，因此比常规的 Python 代码要快。 Ufuncs 支持逐元素处理和类型转换，这意味着更少的循环。

另见

outer()通用函数的文档

使用`chararray`执行字符串操作

NumPy 具有保存字符串的专用chararray对象。它是ndarray的子类，并具有特殊的字符串方法。我们将从 Python 网站下载文本并使用这些方法。 chararray相对于普通字符串数组的优点如下：

索引时会自动修剪数组元素的空白
字符串末尾的空格也被比较运算符修剪
向量化字符串操作可用，因此不需要循环

操作步骤

让我们创建字符数组：

创建字符数组作为视图：

carray = np.array(html).view(np.chararray)

使用expandtabs()函数将制表符扩展到空格。此函数接受制表符大小作为参数。如果未指定，则值为8：
```
carray = carray.expandtabs(1)
```

使用splitlines()函数将行分割成几行：

carray = carray.splitlines()

以下是此示例的完整代码：

import urllib2
import numpy as np
import re

response = urllib2.urlopen('http://python.org/')
html = response.read()
html = re.sub(r'<.*?>', '', html)
carray = np.array(html).view(np.chararray)
carray = carray.expandtabs(1)
carray = carray.splitlines()
print(carray)

工作原理

我们看到了专门的chararray类在起作用。它提供了一些向量化的字符串操作以及有关空格的便捷行为。

另见

chararray类的文档

创建遮罩数组

遮罩数组可用于忽略丢失或无效的数据项。 numpy.ma模块中的MaskedArray类是ndarray的子类，带有遮罩。我们将使用 Lena 图像作为数据源，并假装其中一些数据已损坏。最后，我们将绘制原始图像，原始图像的对数值，遮罩数组及其对数值。

操作步骤

让我们创建被屏蔽的数组：

要创建一个遮罩数组，我们需要指定一个遮罩。创建一个随机遮罩，其值为0或1：
```
random_mask = np.random.randint(0, 2, size=lena.shape)
```

使用上一步中的遮罩，创建一个遮罩数组：

masked_array = np.ma.array(lena, mask=random_mask)

以下是此遮罩数组教程的完整代码：

from __future__ import print_function
import numpy as np
from scipy.misc import lena
import matplotlib.pyplot as plt

lena = lena()
random_mask = np.random.randint(0, 2, size=lena.shape)

plt.subplot(221)
plt.title("Original")
plt.imshow(lena)
plt.axis('off')

masked_array = np.ma.array(lena, mask=random_mask)
print(masked_array)
plt.subplot(222)
plt.title("Masked")
plt.imshow(masked_array)
plt.axis('off')
plt.subplot(223)
plt.title("Log")
plt.imshow(np.log(lena))
plt.axis('off')

plt.subplot(224)
plt.title("Log Masked")
plt.imshow(np.log(masked_array))
plt.axis('off')

plt.show()

这是显示结果图像的屏幕截图：

工作原理

我们对 NumPy 数组应用了随机的遮罩。这具有忽略对应于遮罩的数据的效果。您可以在numpy.ma 模块中找到一系列遮罩数组操作。在本教程中，我们仅演示了如何创建遮罩数组。

另见

numpy.ma模块的文档

忽略负值和极值

当我们想忽略负值时，例如当取数组值的对数时，屏蔽的数组很有用。遮罩数组的另一个用例是排除极值。这基于极限值的上限和下限。

我们将把这些技术应用于股票价格数据。我们将跳过前面几章已经介绍的下载数据的步骤。

操作步骤

我们将使用包含负数的数组的对数：

创建一个数组，该数组包含可被三除的数字：

triples = np.arange(0, len(close), 3)
print("Triples", triples[:10], "...")

接下来，使用与价格数据数组大小相同的数组创建一个数组：

signs = np.ones(len(close))
print("Signs", signs[:10], "...")

借助您在第 2 章，“高级索引和数组概念”中学习的索引技巧，将每个第三个数字设置为负数。

signs[triples] = -1
print("Signs", signs[:10], "...")

最后，取该数组的对数：

ma_log = np.ma.log(close * signs)
print("Masked logs", ma_log[:10], "...")

这应该为AAPL打印以下输出：

Triples [ 0  3  6  9 12 15 18 21 24 27] ...
Signs [ 1\.  1\.  1\.  1\.  1\.  1\.  1\.  1\.  1\.  1.] ...
Signs [-1\.  1\.  1\. -1\.  1\.  1\. -1\.  1\.  1\. -1.] ...
Masked logs [-- 5.93655586575 5.95094223368 -- 5.97468290742 5.97510711452 --
 6.01674381162 5.97889061623 --] ...

让我们将极值定义为低于平均值的一个标准差，或高于平均值的一个标准差（这仅用于演示目的）。编写以下代码以屏蔽极值：

dev = close.std()
avg = close.mean()
inside = numpy.ma.masked_outside(close, avg - dev, avg + dev)
print("Inside", inside[:10], "...")

此代码显示前十个元素：

Inside [-- -- -- -- -- -- 409.429675172 410.240597855 -- --] ...

绘制原始价格数据，绘制对数后的数据，再次绘制指数，最后绘制基于标准差的遮罩后的数据。以下屏幕截图显示了结果（此运行）：

本教程的完整程序如下：

from __future__ import print_function
import numpy as np
from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import matplotlib.pyplot as plt
def get_close(ticker):
   today = date.today()
   start = (today.year - 1, today.month, today.day)

   quotes = quotes_historical_yahoo(ticker, start, today)

   return np.array([q[4] for q in quotes])

close = get_close('AAPL')

triples = np.arange(0, len(close), 3)
print("Triples", triples[:10], "...")

signs = np.ones(len(close))
print("Signs", signs[:10], "...")

signs[triples] = -1
print("Signs", signs[:10], "...")

ma_log = np.ma.log(close * signs)
print("Masked logs", ma_log[:10], "...")

dev = close.std()
avg = close.mean()
inside = np.ma.masked_outside(close, avg - dev, avg + dev)
print("Inside", inside[:10], "...")

plt.subplot(311)
plt.title("Original")
plt.plot(close)

plt.subplot(312)
plt.title("Log Masked")
plt.plot(np.exp(ma_log))

plt.subplot(313)
plt.title("Not Extreme")
plt.plot(inside)

plt.tight_layout()
plt.show()

工作原理

numpy.ma模块中的函数掩盖了数组元素，我们认为这些元素是非法的。例如，log()和sqrt()函数不允许使用负值。屏蔽值类似于数据库和编程中的NULL或None值。具有屏蔽值的所有操作都将导致屏蔽值。

另见

numpy.ma模块的文档

使用`recarray`函数创建得分表

recarray类是ndarray的子类。这些数组可以像数据库中一样保存记录，具有不同的数据类型。例如，我们可以存储有关员工的记录，其中包含诸如薪水之类的数字数据和诸如员工姓名之类的字符串。

现代经济理论告诉我们，投资归结为优化风险和回报。风险是由对数回报的标准差表示的。另一方面，奖励由对数回报的平均值表示。我们可以拿出相对分数，高分意味着低风险和高回报。这只是理论上的，未经测试，所以不要太在意。我们将计算几只股票的得分，并将它们与股票代号一起使用 NumPy recarray()函数中的表格格式存储。

操作步骤

让我们从创建记录数组开始：

为每个记录创建一个包含符号，标准差得分，平均得分和总得分的记录数组：

weights = np.recarray((len(tickers),), dtype=[('symbol', np.str_, 16), 
    ('stdscore', float), ('mean', float), ('score', float)])

为了简单起见，请根据对数收益在循环中初始化得分：

for i, ticker in enumerate(tickers):
   close = get_close(ticker)
   logrets = np.diff(np.log(close))
   weights[i]['symbol'] = ticker
   weights[i]['mean'] = logrets.mean()
   weights[i]['stdscore'] = 1/logrets.std()
   weights[i]['score'] = 0

如您所见，我们可以使用在上一步中定义的字段名称来访问元素。

现在，我们有一些数字，但是它们很难相互比较。归一化分数，以便我们以后可以将它们合并。在这里，归一化意味着确保分数加起来为：
```
for key in ['mean', 'stdscore']:
    wsum = weights[key].sum()
    weights[key] = weights[key]/wsum
```

总体分数将只是中间分数的平均值。对总分上的记录进行排序以产生排名：

weights['score'] = (weights['stdscore'] + weights['mean'])/2
weights['score'].sort()

The following is the complete code for this example:

from __future__ import print_function
import numpy as np
from matplotlib.finance import quotes_historical_yahoo
from datetime import date

tickers = ['MRK', 'T', 'VZ']

def get_close(ticker):
   today = date.today()
   start = (today.year - 1, today.month, today.day)

   quotes = quotes_historical_yahoo(ticker, start, today)

   return np.array([q[4] for q in quotes])

weights = np.recarray((len(tickers),), dtype=[('symbol', np.str_, 16), 
   ('stdscore', float), ('mean', float), ('score', float)])

for i, ticker in enumerate(tickers):
   close = get_close(ticker)
   logrets = np.diff(np.log(close))
   weights[i]['symbol'] = ticker
   weights[i]['mean'] = logrets.mean()
   weights[i]['stdscore'] = 1/logrets.std()
   weights[i]['score'] = 0

for key in ['mean', 'stdscore']:
   wsum = weights[key].sum()
   weights[key] = weights[key]/wsum

weights['score'] = (weights['stdscore'] + weights['mean'])/2
weights['score'].sort()

for record in weights:
   print("%s,mean=%.4f,stdscore=%.4f,score=%.4f" % (record['symbol'], record['mean'], record['stdscore'], record['score']))

该程序产生以下输出：

MRK,mean=0.8185,stdscore=0.2938,score=0.2177
T,mean=0.0927,stdscore=0.3427,score=0.2262
VZ,mean=0.0888,stdscore=0.3636,score=0.5561

分数已归一化，因此值介于0和1之间，我们尝试从秘籍开始使用定义获得最佳收益和风险组合。根据输出，VZ得分最高，因此是最好的投资。当然，这只是一个 NumPy 演示，数据很少，所以不要认为这是推荐。

工作原理

我们计算了几只股票的得分，并将它们存储在recarray NumPy 对象中。这个数组使我们能够混合不同数据类型的数据，在这种情况下，是股票代码和数字得分。记录数组使我们可以将字段作为数组成员访问，例如arr.field。本教程介绍了记录数组的创建。您可以在numpy.recarray模块中找到更多与记录数组相关的功能。

另见

numpy.recarray模块的文档

七、性能分析和调试

在本章中，我们将介绍以下秘籍：

使用timeit进行性能分析
使用 IPython 进行分析
安装line_profiler
使用line_profiler分析代码
具有cProfile扩展名的性能分析代码
使用 IPython 进行调试
使用PuDB进行调试

简介

调试是从软件中查找和删除错误的行为。分析是指构建程序的概要文件，以便收集有关内存使用或时间复杂度的信息。分析和调试是开发人员生活中必不可少的活动。对于复杂的软件尤其如此。好消息是，许多工具可以为您提供帮助。我们将回顾 NumPy 用户中流行的技术。

使用`timeit`进行性能分析

timeit是一个模块，可用于计时代码段。它是标准 Python 库的一部分。我们将使用几种数组大小对sort() NumPy 函数计时。经典的快速排序和归并排序算法的平均运行时间为O(N log N)，因此我们将尝试将这个模型拟合到结果。

操作步骤

我们将要求数组进行排序：

创建数组以排序包含随机整数值的各种大小：

times = np.array([])

for size in sizes:
    integers = np.random.random_integers (1, 10 ** 6, size)

要测量时间，请创建一个计时器，为其提供执行函数，并指定相关的导入。然后，排序 100 次以获取有关排序时间的数据：
```
def measure():
    timer = timeit.Timer('dosort()', 'from __main__ import dosort')
    return timer.timeit(10 ** 2)
```
通过一次乘以时间来构建测量时间数组：
```
times = np.append(times, measure())
```

将时间拟合为n log n的理论模型。由于我们将数组大小更改为 2 的幂，因此很容易：

fit = np.polyfit(sizes * powersOf2, times, 1)

以下是完整的计时代码：

import numpy as np
import timeit
import matplotlib.pyplot as plt

# This program measures the performance of the NumPy sort function
# and plots time vs array size.
integers = []

def dosort():
   integers.sort()

def measure():
   timer = timeit.Timer('dosort()', 'from __main__ import dosort')

   return timer.timeit(10 ** 2)

powersOf2 = np.arange(0, 19)
sizes = 2 ** powersOf2

times = np.array([])

for size in sizes:
   integers = np.random.random_integers(1, 10 ** 6, size)
   times = np.append(times, measure())

fit = np.polyfit(sizes * powersOf2, times, 1)
print(fit)
plt.title("Sort array sizes vs execution times")
plt.xlabel("Size")
plt.ylabel("(s)")
plt.semilogx(sizes, times, 'ro')
plt.semilogx(sizes, np.polyval(fit, sizes * powersOf2))
plt.grid()
plt.show()

以下屏幕截图显示了运行时间与数组大小的关系图：

工作原理

我们测量了sort() NumPy 函数的平均运行时间。此秘籍中使用了以下函数：

函数	描述
`random_integers()`	给定值和数组大小的范围时，此函数创建一个随机整数数组
`append()`	此函数将值附加到 NumPy 数组
`polyfit()`	此函数将数据拟合为给定阶数的多项式
`polyval()`	此函数计算多项式，并为给定的 `x`值返回相应的值
`semilogx()`	此函数使用对数刻度在 X 轴上绘制数据

另见

timeit的文档

使用 IPython 进行分析

在 IPython 中，我们可以使用timeit来分析代码的小片段。我们可以也分析较大的脚本。我们将展示两种方法。

操作步骤

首先，我们将介绍一个小片段：

以pylab模式启动 IPython：
```
$ ipython --pylab
```
创建一个包含 1000 个介于 0 到 1000 之间的整数值的数组：
```
In [1]: a = arange(1000)
```
测量在数组中搜索“所有问题的答案”（42）所花费的时间。是的，所有问题的答案都是 42。如果您不相信我，请阅读这个页面：
```
In [2]: %timeit searchsorted(a, 42)
100000 loops, best of 3: 7.58 us per loop
```

剖析以下小脚本，该小脚本可以反转包含随机值的大小可变的矩阵。 NumPy 矩阵的.I属性（即大写I）表示该矩阵的逆：

import numpy as np

def invert(n):
  a = np.matrix(np.random.rand(n, n))

  return a.I

sizes = 2 ** np.arange(0, 12)

for n in sizes:
  invert(n)

将此代码计时如下：

In [1]: %run -t invert_matrix.py

IPython CPU timings (estimated):
 User   :       6.08 s.
 System :       0.52 s.
Wall time:      19.26 s.

然后使用p选项对脚本进行配置：

In [2]: %run -p invert_matrix.py

852 function calls in 6.597 CPU seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
       12    3.228    0.269    3.228    0.269 {numpy.linalg.lapack_lite.dgesv}
       24    2.967    0.124    2.967    0.124 {numpy.core.multiarray._fastCopyAndTranspose}
       12    0.156    0.013    0.156    0.013 {method 'rand' of 'mtrand.RandomState' objects}
       12    0.087    0.007    0.087    0.007 {method 'copy' of 'numpy.ndarray' objects}
       12    0.069    0.006    0.069    0.006 {method 'astype' of 'numpy.ndarray' objects}
       12    0.025    0.002    6.304    0.525 linalg.py:404(inv)
       12    0.024    0.002    6.328    0.527 defmatrix.py:808(getI)
        1    0.017    0.017    6.596    6.596 invert_matrix.py:1(<module>)
       24    0.014    0.001    0.014    0.001 {numpy.core.multiarray.zeros}
       12    0.009    0.001    6.580    0.548 invert_matrix.py:3(invert)
       12    0.000    0.000    6.264    0.522 linalg.py:244(solve)
       12    0.000    0.000    0.014    0.001 numeric.py:1875(identity)
        1    0.000    0.000    6.597    6.597 {execfile}
       36    0.000    0.000    0.000    0.000 defmatrix.py:279(__array_finalize__)
       12    0.000    0.000    2.967    0.247 linalg.py:139(_fastCopyAndTranspose)
       24    0.000    0.000    0.087    0.004 defmatrix.py:233(__new__)
       12    0.000    0.000    0.000    0.000 linalg.py:99(_commonType)
       24    0.000    0.000    0.000    0.000 {method '__array_prepare__' of 'numpy.ndarray' objects}
       36    0.000    0.000    0.000    0.000 linalg.py:66(_makearray)
       36    0.000    0.000    0.000    0.000 {numpy.core.multiarray.array}
       12    0.000    0.000    0.000    0.000 {method 'view' of 'numpy.ndarray' objects}
       12    0.000    0.000    0.000    0.000 linalg.py:127(_to_native_byte_order)
        1    0.000    0.000    6.597    6.597 interactiveshell.py:2270(safe_execfile)

工作原理

我们通过分析器运行了上述 NumPy 代码。下表概述了分析器的输出：

函数	描述
`ncalls`	这是调用次数
`tottime`	这是一个函数花费的总时间
`percall`	这是每次通话所花费的时间，计算方法是将总时间除以通话次数
`cumtime`	这是在函数和由函数调用的函数（包括递归调用）上花费的累积时间

另见

IPython 魔术文档

安装`line_profiler`

line_profiler由 NumPy 的开发人员之一创建。此模块对 Python 代码进行逐行分析。我们将在此秘籍中描述必要的安装步骤。

准备

您可能需要安装setuptools。先前的秘籍中对此进行了介绍；如有必要，请参阅“另见”部分。为了安装开发版本，您将需要 Git。安装 Git 超出了本书的范围。

操作步骤

选择适合您的安装选项：

使用以下任一命令将line_profiler与easy_install一起安装：
```
$ easy_install line_profiler
$ pip install line_profiler
```
安装开发版本。

使用 Git 查看源代码：
```
$ git clone https://github.com/rkern/line_profiler
```
签出源代码后，按如下所示构建它：
```
$ python setup.py install
```

另见

第 1 章，“使用 IPython”中的“安装 IPython”

使用`line_profiler`分析代码

现在我们已经安装完毕，可以开始分析。

操作步骤

显然，我们将需要代码来分析：

编写以下代码，以自身乘以大小可变的随机矩阵。此外，线程将休眠几秒钟。使用@profile注解函数以进行概要分析：

import numpy as np
import time

@profile
def multiply(n):
  A = np.random.rand(n, n)
  time.sleep(np.random.randint(0, 2))
  return np.matrix(A) ** 2

for n in 2 ** np.arange(0, 10):
  multiply(n)

使用以下命令运行事件分析器：

$ kernprof.py -l -v mat_mult.py
Wrote profile results to mat_mult.py.lprof
Timer unit: 1e-06 s

File: mat_mult.py
Function: multiply at line 4
Total time: 3.19654 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
 4                                           @profile
 5                                           def multiply(n):
 6        10        13461   1346.1      0.4     A = numpy.random.rand(n, n)
 7        10      3000689 300068.9     93.9     time.sleep(numpy.random.randint(0, 2))
 8        10       182386  18238.6      5.7     return numpy.matrix(A) ** 2

工作原理

@profile装饰器告诉line_profiler要分析哪些函数。下表说明了分析器的输出：

函数	描述
`Line #`	文件中的行号
`Hits`	执行该行的次数
`Time`	执行该行所花费的时间
`Per Hit`	执行该行所花费的平均时间
`% Time`	执行该行所花费的时间相对于执行所有行所花费的时间的百分比
`Line Contents`	该行的内容

另见

Github line_profiler项目页面

`cProfile`扩展和代码性能分析

cProfile是 Python 2.5 中引入的C扩展名。它可以用于确定性分析。确定性分析表示所获得的时间测量是精确的，并且不使用采样。这与统计分析相反，统计分析来自随机样本。我们将使用cProfile对一个小的 NumPy 程序进行分析，该程序会对具有随机值的数组进行转置。

操作步骤

同样，我们需要代码来配置：

编写以下transpose()函数以创建具有随机值的数组并将其转置：

def transpose(n):
  random_values = np.random.random((n, n))
  return random_values.T

运行分析器，并为其提供待分析函数：

cProfile.run('transpose (1000)')

可以在以下片段中找到本教程的完整代码：

import numpy as np
import cProfile

def transpose(n):
   random_values = np.random.random((n, n))
   return random_values.T

cProfile.run('transpose (1000)')

对于1000 x 1000的数组，我们得到以下输出：

4 function calls in 0.029 CPU seconds
 Ordered by: standard name
 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 1    0.001    0.001    0.029    0.029 <string>:1(<module>)
 1    0.000    0.000    0.028    0.028 cprofile_transpose.py:5(transpose)
 1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
 1    0.028    0.028    0.028    0.028 {method 'random_sample' of 'mtrand.RandomState' objects}

输出中的列与 IPython 分析秘籍中看到的列相同。

另见

使用 IPython 进行调试

“如果调试是清除软件错误的过程，则编程必须是放入它们的过程。”

-- 荷兰计算机科学家 Edsger Dijkstra，1972 年图灵奖的获得者

调试是没人真正喜欢，但是掌握这些东西非常重要的东西之一。这可能需要几个小时，并且由于墨菲定律，您很可能没有时间。因此，重要的是要系统地了解您的工具。找到错误并实现修复后，您应该进行单元测试（如果该错误具有来自问题跟踪程序的相关 ID，我通常在末尾附加 ID 来命名测试）。这样，您至少不必再次进行调试。下一章将介绍单元测试。我们将调试以下错误代码。它尝试访问不存在的数组元素：

import numpy as np

a = np.arange(7)
print(a[8])

IPython 调试器充当普通的 Python pdb调试器；它添加了选项卡补全和语法突出显示等功能。

操作步骤

以下步骤说明了典型的调试会话：

启动 IPython Shell。通过发出以下命令在 IPython 中运行错误脚本：

In [1]: %run buggy.py
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
.../site-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
 173             else:
 174                 filename = fname
--> 175             __builtin__.execfile(filename, *where)

.../buggy.py in <module>()
 2
 3 a = numpy.arange(7)
----> 4 print a[8]

IndexError: index out of bounds

现在您的程序崩溃了，启动调试器。在发生错误的行上设置一个断点：
```
In [2]: %debug
> .../buggy.py(4)<module>()
 2 
 3 a = numpy.arange(7)
----> 4 print a[8]
```

使用list命令列出代码，或使用简写l：

ipdb> list
 1 import numpy as np
 2 
 3 a = np.arange(7)
----> 4 print(a[8])

现在，我们可以在调试器当前所在的行上求值任意代码：
```
ipdb> len(a)
7

ipdb> print(a)
[0 1 2 3 4 5 6]
```

调用栈是包含有关正在运行的程序的活动函数的信息的栈。使用bt命令查看调用栈：

ipdb> bt
 .../py3compat.py(175)execfile()
 171             if isinstance(fname, unicode):
 172                 filename = fname.encode(sys.getfilesystemencoding())
 173             else:
 174                 filename = fname
--> 175             __builtin__.execfile(filename, *where)

> .../buggy.py(4)<module>()
 0 print a[8]

向上移动调用栈：

ipdb> u
> .../site-packages/IPython/utils/py3compat.py(175)execfile()
 173             else:
 174                 filename = fname
--> 175             __builtin__.execfile(filename, *where)

下移调用栈：

ipdb> d
> .../buggy.py(4)<module>()
 2
 3 a = np.arange(7)
----> 4 print(a[8])

工作原理

在本教程中，您学习了如何使用 IPython 调试 NumPy 程序。我们设置一个断点并导航调用栈。使用了以下调试器命令：

函数	描述
`list`或 `l`	列出源代码
`bt`	显示调用栈
`u`	向上移动调用栈
`d`	下移调用栈

另见

使用 PuDB 进行调试

PuDB 是基于视觉的，全屏，基于控制台的 Python 调试器，易于安装。 PuDB 支持光标键和 vi 命令。如果需要，我们还可以将此调试器与 IPython 集成。

操作步骤

我们将从安装pudb开始：

要安装pudb，我们只需执行以下命令（或等效的pip命令）：

$ sudo easy_install pudb
$ pip install pudb
$ pip freeze|grep pudb
pudb==2014.1

让我们调试前面示例中的buggy程序。如下所示启动调试器：
```
$ python -m pudb buggy.py
```
以下屏幕截图显示了调试器的用户界面：

屏幕快照在顶部显示了最重要的调试命令。我们还可以看到正在调试的代码，变量，栈和定义的断点。键入q退出大多数菜单。键入n将调试器移至下一行。我们还可以使用光标键或 vi 的j和k键移动，例如，通过键入b设置断点。

另见

PyPi PuDB 页面

八、质量保证

“如果您对计算机撒谎，它将帮助您。”

-- Perry Farrar，ACM 通讯，第 28 卷

在本章中，我们将介绍以下秘籍：

安装 Pyflakes
使用 Pyflakes 执行静态分析
用 Pylint 分析代码
使用 Pychecker 执行静态分析
使用docstrings测试代码
编写单元测试
使用模拟测试代码
以 BDD 方式来测试

简介

与普遍的看法相反，质量保证与其说是发现错误，不如说是发现它们。我们将讨论两种提高代码质量，从而防止出现问题的方法。首先，我们将对已经存在的代码进行静态分析。然后，我们将讨论单元测试；这包括模拟和行为驱动开发（BDD）。

安装 Pyflakes

Pyflakes 是 Python 代码分析包。它可以分析代码并发现潜在的问题，例如：

未使用的导入
未使用的变量

准备

如有必要，请安装pip或easy_install。

操作步骤

选择以下之一来安装pyflakes：

使用pip命令安装 pyflakes：
```
$ sudo pip install pyflakes
```
使用easy_install命令安装 Pyflakes：
```
$ sudo easy_install pyflakes
```
这是在 Linux 上安装此包的两种方法：

Linux 包的名称也为pyflakes。例如，在 RedHat 上执行以下操作：
```
$ sudo yum install pyflakes
```
在 Debian/Ubuntu 上，命令如下：
```
$ sudo apt-get install pyflakes
```

另见

Pyflakes 主页

使用 Pyflakes 执行静态分析

我们将对 NumPy 代码库的一部分执行静态分析。为此，我们将使用 Git 签出代码。然后，我们将使用pyflakes对部分代码进行静态分析。

操作步骤

要检查 NumPy 代码中，我们需要 Git。安装 Git 超出了本书的范围：

用 Git 命令检索代码如下：
```
$ git clone git://github.com/numpy/numpy.git numpy
```
或者，从这里下载源档案。

上一步使用完整的 NumPy 代码创建一个numpy目录。转到此目录，并在其中运行以下命令：

$ pyflakes *.py
pavement.py:71: redefinition of unused 'md5' from line 69
pavement.py:88: redefinition of unused 'GIT_REVISION' from line 86
pavement.py:314: 'virtualenv' imported but unused
pavement.py:315: local variable 'e' is assigned to but never used
pavement.py:380: local variable 'sdir' is assigned to but never used
pavement.py:381: local variable 'bdir' is assigned to but never used
pavement.py:536: local variable 'st' is assigned to but never used
setup.py:21: 're' imported but unused
setup.py:27: redefinition of unused 'builtins' from line 25
setup.py:124: redefinition of unused 'GIT_REVISION' from line 118
setupegg.py:17: 'setup' imported but unused
setupscons.py:61: 'numpy' imported but unused
setupscons.py:64: 'numscons' imported but unused
setupsconsegg.py:6: 'setup' imported but unused

这将对代码样式进行分析，并检查当前目录中所有 Python 脚本中的 PEP-8 违规情况。如果愿意，还可以分析单个文件。

工作原理

正如您所见，分析代码样式并使用 Pyflakes 查找违反 PEP-8 的行为非常简单。另一个优点是速度。但是，Pyflakes 报告的错误类型的数量是有限的。

使用 Pylint 分析代码

Pylint 是另一个由 Logilab 创建的开源静态分析器。 Pylint 比 Pyflakes 更复杂；它允许更多的自定义和代码检查。但是，它比 Pyflakes 慢。有关更多信息，请参见手册。

在本秘籍中，我们再次从 Git 存储库下载 NumPy 代码-为简便起见，省略了此步骤。

准备

您可以从源代码发行版中安装 Pylint。但是，有很多依赖项，因此最好使用easy_install或pip进行安装。安装命令如下：

$ easy_install pylint
$ sudo pip install pylint

操作步骤

我们将再次从 NumPy 代码库的顶部目录进行分析。注意，我们得到了更多的输出。实际上，Pylint 打印了太多文本，因此在这里大部分都必须省略：

$ pylint *.py
No config file found, using default configuration
************* Module pavement
C: 60: Line too long (81/80)
C:139: Line too long (81/80)
...
W: 50: TODO
W:168: XXX: find out which env variable is necessary to avoid the pb with python
W: 71: Reimport 'md5' (imported line 143)
F: 73: Unable to import 'paver'
F: 74: Unable to import 'paver.easy'
C: 79: Invalid name "setup_py" (should match (([A-Z_][A-Z0-9_]*)|(__.*__))$)
F: 86: Unable to import 'numpy.version'
E: 86: No name 'version' in module 'numpy'
C:149: Operator not followed by a space
if sys.platform =="darwin":
 ^^
C:202:prepare_nsis_script: Missing docstring
W:228:bdist_superpack: Redefining name 'options' from outer scope (line 74)
C:231:bdist_superpack.copy_bdist: Missing docstring
W:275:bdist_wininst_nosse: Redefining name 'options' from outer scope (line 74)

工作原理

Pylint 默认输出原始文本；但是我们可以根据需要请求 HTML 输出。消息具有以下格式：

MESSAGE_TYPE: LINE_NUM:[OBJECT:] MESSAGE

消息类型可以是以下之一：

[R]：这意味着建议进行重构
[C]：这意味着存在违反代码风格的情况
[W]：用于警告小问题
[E]：用于错误或潜在的错误
[F]：这表明发生致命错误，阻止了进一步的分析

另见

使用 Pyflakes 执行静态分析

使用 Pychecker 执行静态分析

Pychecker 是一个古老的静态分析工具。它不是十分活跃的开发工具，但它在此提到的速度又足够好。在编写本书时，最新版本是 0.8.19，最近一次更新是在 2011 年。Pychecker 尝试导入每个模块并对其进行处理。然后，它搜索诸如传递不正确数量的参数，使用不存在的方法传递不正确的格式字符串以及其他问题之类的问题。在本秘籍中，我们将再次分析代码，但是这次使用 Pychecker。

操作步骤

从 Sourceforge 下载tar.gz。解压缩源归档文件并运行以下命令：

$ python setup.py install

或者，使用pip安装 Pychecker：

$ sudo pip install http://sourceforge.net/projects/pychecker/files/pychecker/0.8.19/pychecker-0.8.19.tar.gz/download

分析代码，就像先前的秘籍一样。我们需要的命令如下：

$ pychecker *.py
...
Warnings...

...

setup.py:21: Imported module (re) not used
setup.py:27: Module (builtins) re-imported

...

使用文档字符串测试代码

Doctests 是注释字符串，它们嵌入在类似交互式会话的 Python 代码中。这些字符串可用于测试某些假设或仅提供示例。我们需要使用doctest模块来运行这些测试。

让我们写一个简单的示例，该示例应该计算阶乘，但不涵盖所有可能的边界条件。换句话说，某些测试将失败。

操作步骤

用将通过的测试和将失败的另一个测试编写docstring。 docstring文本应类似于在 Python shell 中通常看到的文本：

"""
Test for the factorial of 3 that should pass.
>>> factorial(3)
6

Test for the factorial of 0 that should fail.
>>> factorial(0)
1
"""

编写以下 NumPy 代码：
```
return np.arange(1, n+1).cumprod()[-1]
```
我们希望这段代码有时会故意失败。它将创建一个序列号数组，计算该数组的累积乘积，并返回最后一个元素。

使用doctest模块运行测试：

doctest.testmod()

以下是本书代码包中docstringtest.py文件的完整测试示例代码：

import numpy as np
import doctest

def factorial(n):
   """
   Test for the factorial of 3 that should pass.
   >>> factorial(3)
   6

   Test for the factorial of 0 that should fail.
   >>> factorial(0)
   1
   """
   return np.arange(1, n+1).cumprod()[-1]

doctest.testmod()

我们可以使用-v选项获得详细的输出，如下所示：

$ python docstringtest.py -v
Trying:
 factorial(3)
Expecting:
 6
ok
Trying:
 factorial(0)
Expecting:
 1
**********************************************************************
File "docstringtest.py", line 11, in __main__.factorial
Failed example:
 factorial(0)
Exception raised:
 Traceback (most recent call last):
 File ".../doctest.py", line 1253, in __run
 compileflags, 1) in test.globs
 File "<doctest __main__.factorial[1]>", line 1, in <module>
 factorial(0)
 File "docstringtest.py", line 14, in factorial
 return numpy.arange(1, n+1).cumprod()[-1]
 IndexError: index out of bounds
1 items had no tests:
 __main__
**********************************************************************
1 items had failures:
 1 of   2 in __main__.factorial
2 tests in 2 items.
1 passed and 1 failed.
***Test Failed*** 1 failures.

工作原理

如您所见，我们没有考虑零和负数。实际上，由于数组为空，我们出现了index out of bounds错误。当然，这很容易解决，我们将在下一个教程中进行。

另见

doctest官方文档

编写单元测试

测试驱动开发（TDD）是本世纪软件开发诞生的最好的事情。 TDD 的最重要方面之一是，几乎把重点放在单元测试上。

注意

TDD 方法使用所谓的测试优先方法，在此方法中，我们首先编写一个失败的测试，然后编写相应的代码以通过测试。测试应记录开发人员的意图，但要比功能设计的水平低。一组测试通过降低回归概率来增加置信度，并促进重构。

单元测试是自动测试，通常测试一小段代码，通常是一个函数或方法。 Python 具有用于单元测试的 PyUnit API。作为 NumPy 的用户，我们也可以使用numpy.testing模块中的便捷函数。顾名思义，该模块专用于测试。

操作步骤

让我们编写一些代码进行测试：

首先编写以下factorial()函数：

def factorial(n):
  if n == 0:
    return 1

  if n < 0:
    raise ValueError, "Don't be so negative"

  return np.arange(1, n+1).cumprod()

该代码与前面的秘籍中的代码相同，但是我们添加了一些边界条件检查。

让我们写一个类；此类将包含单元测试。它从unittest模块扩展了TestCase类，是 Python 标准测试的一部分。我们通过调用factorial()函数并运行以下代码来运行测试：

一个正数-幸福的道路！
边界条件等于0

负数，这将导致错误：

class FactorialTest(unittest.TestCase):
    def test_factorial(self):
      #Test for the factorial of 3 that should pass.
      self.assertEqual(6, factorial(3)[-1])
      np.testing.assert_equal(np.array([1, 2, 6]), factorial(3))

    def test_zero(self):
      #Test for the factorial of 0 that should pass.
      self.assertEqual(1, factorial(0))

    def test_negative(self):
      #Test for the factorial of negative numbers that should fail.
      # It should throw a ValueError, but we expect IndexError
      self.assertRaises(IndexError, factorial(-10))

factorial()函数和整个单元测试的代码如下：

import numpy as np
import unittest

def factorial(n):
   if n == 0:
      return 1

   if n < 0:
      raise ValueError, "Don't be so negative"

   return np.arange(1, n+1).cumprod()

class FactorialTest(unittest.TestCase):
   def test_factorial(self):
      #Test for the factorial of 3 that should pass.
      self.assertEqual(6, factorial(3)[-1])
      np.testing.assert_equal(np.array([1, 2, 6]), factorial(3))

   def test_zero(self):
      #Test for the factorial of 0 that should pass.
      self.assertEqual(1, factorial(0))

   def test_negative(self):
      #Test for the factorial of negative numbers that should fail.
      # It should throw a ValueError, but we expect IndexError
      self.assertRaises(IndexError, factorial(-10))

if __name__ == '__main__':
    unittest.main()

负数测试失败，如以下输出所示：

.E.
======================================================================
ERROR: test_negative (__main__.FactorialTest)
----------------------------------------------------------------------
Traceback (most recent call last):
 File "unit_test.py", line 26, in test_negative
 self.assertRaises(IndexError, factorial(-10))
 File "unit_test.py", line 9, in factorial
 raise ValueError, "Don't be so negative"
ValueError: Don't be so negative

----------------------------------------------------------------------
Ran 3 tests in 0.001s

FAILED (errors=1)

工作原理

我们看到了如何使用标准unittest Python 模块实现简单的单元测试。我们编写了一个测试类，该类从unittest模块扩展了TestCase类。以下函数用于执行各种测试：

函数	描述
`numpy.testing.assert_equal()`	测试两个 NumPy 数组是否相等
`unittest.assertEqual()`	测试两个值是否相等
`unittest.assertRaises()`	测试是否引发异常

testing NumPy 包具有许多我们应该了解的测试函数，如下所示：

函数	描述
`assert_almost_equal()`	如果两个数字不等于指定的精度，则此函数引发异常
`assert_approx_equal()`	如果两个数字在一定意义上不相等，则此函数引发异常
`assert_array_almost_equal()`	如果两个数组不等于指定的精度，此函数会引发异常
`assert_array_equal()`	如果两个数组不相等，则此函数引发异常
`assert_array_less()`	如果两个数组的形状不同，并且此函数引发异常，则第一个数组的元素严格小于第二个数组的元素
`assert_raises()`	如果使用定义的参数调用的可调用对象未引发指定的异常，则此函数将失败
`assert_warns()`	如果未抛出指定的警告，则此函数失败
`assert_string_equal()`	此函数断言两个字符串相等

使用模拟测试代码

模拟是用来代替真实对象的对象，目的是测试真实对象的部分行为。如果您看过电影《身体抢夺者》，您可能已经对基本概念有所了解。一般来说，仅在被测试的真实对象的创建成本很高（例如数据库连接）或测试可能产生不良副作用时才有用。例如，我们可能不想写入文件系统或数据库。

在此秘籍中，我们将测试一个核反应堆，当然不是真正的反应堆！此类核反应堆执行阶乘计算，从理论上讲，它可能导致连锁反应，进而导致核灾难。我们将使用mock包通过模拟来模拟阶乘计算。

操作步骤

首先，我们将安装mock包；之后，我们将创建一个模拟并测试一段代码：

要安装mock包，请执行以下命令：
```
$ sudo easy_install mock
```
核反应堆类有一个do_work()方法，该方法调用了我们要模拟的危险的factorial()方法。创建一个模拟，如下所示：
```
reactor.factorial = MagicMock(return_value=6)
```
这样可以确保模拟返回值6。

我们可以通过多种方式检查模拟的行为，然后从中检查真实对象的行为。例如，断言使用正确的参数调用了潜在爆炸性的factorial()方法，如下所示：

reactor.factorial.assert_called_with(3, "mocked")

带有模拟的完整测试代码如下：

from __future__ import print_function
from mock import MagicMock
import numpy as np
import unittest

class NuclearReactor():
   def __init__(self, n):
      self.n = n

   def do_work(self, msg):
      print("Working")

      return self.factorial(self.n, msg)

   def factorial(self, n, msg):
      print(msg)

      if n == 0:
         return 1

      if n < 0:
         raise ValueError, "Core meltdown"

      return np.arange(1, n+1).cumprod()

class NuclearReactorTest(unittest.TestCase):
   def test_called(self):
      reactor = NuclearReactor(3)
      reactor.factorial = MagicMock(return_value=6)
      result = reactor.do_work("mocked")
      self.assertEqual(6, result)
      reactor.factorial.assert_called_with(3, "mocked")

   def test_unmocked(self):
      reactor = NuclearReactor(3)
      reactor.factorial(3, "unmocked")
      np.testing.assert_raises(ValueError)

if __name__ == '__main__':
    unittest.main()

我们将一个字符串传递给factorial()方法，以显示带有模拟的代码不会执行实际的代码。该单元测试的工作方式与上一秘籍中的单元测试相同。这里的第二项测试不测试任何内容。第二个测试的目的只是演示，如果我们在没有模拟的情况下执行真实代码，会发生什么。

测试的输出如下：

Working
.unmocked
.
----------------------------------------------------------------------
Ran 2 tests in 0.000s

OK

工作原理

模拟没有任何行为。他们就像外星人的克隆人，假装是真实的人。只能比外星人傻—外星人克隆人无法告诉您被替换的真实人物的生日。我们需要设置它们以适当的方式进行响应。例如，在此示例中，模拟返回6 。我们可以记录模拟发生了什么，被调用了多少次以及使用了哪些参数。

另见

Mock 包主页

以 BDD 方式来测试

BDD（行为驱动开发）是您可能遇到的另一个热门缩写。在 BDD 中，我们首先根据某些约定和规则定义（英语）被测系统的预期行为。在本秘籍中，我们将看到这些约定的示例。

这种方法背后的想法是，我们可以让可能无法编程或编写测试大部分内容的人员参加。这些人编写的功能采用句子的形式，包括多个步骤。每个步骤或多或少都是我们可以编写的单元测试，例如，使用 NumPy。有许多 Python BDD 框架。在本秘籍中，我们使用 Lettuce 来测试阶乘函数。

操作步骤

在本节中，您将学习如何安装 Lettuce，设置测试以及编写测试规范：

要安装 Lettuce，请运行以下命令之一：

$ pip install lettuce
$ sudo easy_install lettuce

Lettuce 需要特殊的目录结构进行测试。在tests目录中，我们将有一个名为features的目录，其中包含factorial.feature文件，以及steps.py文件中的功能说明和测试代码：
```
./tests:
features

./tests/features:
factorial.feature	steps.py
```

提出业务需求是一项艰巨的工作。以易于测试的方式将其全部写下来更加困难。幸运的是，这些秘籍的要求非常简单-我们只需写下不同的输入值和预期的输出。我们在Given，When和Then部分中有不同的方案，它们对应于不同的测试步骤。为阶乘函数定义以下三种方案：

Feature: Compute factorial

    Scenario: Factorial of 0
      Given I have the number 0 
      When I compute its factorial 
      Then I see the number 1

    Scenario: Factorial of 1
      Given I have the number 1 
      When I compute its factorial 
      Then I see the number 1

    Scenario: Factorial of 3
      Given I have the number 3 
      When I compute its factorial 
      Then I see the number 1, 2, 6

我们将定义与场景步骤相对应的方法。要特别注意用于注释方法的文本。它与业务场景文件中的文本匹配，并且我们使用正则表达式获取输入参数。在前两个方案中，我们匹配数字，在最后一个方案中，我们匹配任何文本。 fromstring() NumPy 函数用于从 NumPy 数组创建字符串，字符串中使用整数数据类型和逗号分隔符。以下代码测试了我们的方案：

from lettuce import *
import numpy as np

@step('I have the number (\d+)')
def have_the_number(step, number):
    world.number = int(number)

@step('I compute its factorial')
def compute_its_factorial(step):
    world.number = factorial(world.number)

@step('I see the number (.*)')
def check_number(step, expected):
    expected = np.fromstring(expected, dtype=int, sep=',')
    np.testing.assert_equal(world.number, expected, \
        "Got %s" % world.number)

def factorial(n):
   if n == 0:
      return 1

   if n < 0:
      raise ValueError, "Core meltdown"

   return np.arange(1, n+1).cumprod()

要运行测试，请转到tests目录，然后键入以下命令：

$ lettuce

Feature: Compute factorial        # features/factorial.feature:1

 Scenario: Factorial of 0        # features/factorial.feature:3
 Given I have the number 0     # features/steps.py:5
 When I compute its factorial  # features/steps.py:9
 Then I see the number 1       # features/steps.py:13

 Scenario: Factorial of 1        # features/factorial.feature:8
 Given I have the number 1     # features/steps.py:5
 When I compute its factorial  # features/steps.py:9
 Then I see the number 1       # features/steps.py:13

 Scenario: Factorial of 3        # features/factorial.feature:13
 Given I have the number 3     # features/steps.py:5
 When I compute its factorial  # features/steps.py:9
 Then I see the number 1, 2, 6 # features/steps.py:13

1 feature (1 passed)
3 scenarios (3 passed)
9 steps (9 passed)

工作原理

我们定义了具有三个方案和相应步骤的函数。我们使用 NumPy 的测试函数来测试不同步骤，并使用fromstring()函数从规格文本创建 NumPy 数组。

另见

Lettuce 文档

九、使用 Cython 加速代码

在本章中，我们将介绍以下秘籍：

安装 Cython
构建 HelloWorld 程序
将 Cython 与 NumPy 结合使用
调用 C 函数
分析 Cython 代码
用 Cython 近似阶乘

简介

Cython 是基于 Python 的相对年轻的编程语言。它允许编码人员将 C 的速度与 Python 的功能混合在一起。与 Python 的区别在于我们可以选择声明静态类型。许多编程语言（例如 C）具有静态类型，这意味着我们必须告诉 C 变量的类型，函数参数和返回值类型。另一个区别是 C 是一种编译语言，而 Python 是一种解释语言。根据经验，可以说 C 比 Python 更快，但灵活性更低。通过 Cython 代码，我们可以生成 C 或 C++ 代码。之后，我们可以将生成的代码编译为 Python 扩展模块。

在本章中，您将学习 Cython。我们将获得一些与 NumPy 一起运行的简单 Cython 程序。另外，我们将介绍 Cython 代码。

安装 Cython

为了使用 Cython，我们需要安装它。 Enthought Canopy，Anaconda 和 Sage 发行版包括 Cython。有关更多信息，请参见这里，这里和这里。我们将不在这里讨论如何安装这些发行版。显然，我们需要一个 C 编译器来编译生成的 C 代码。在某些操作系统（例如 Linux）上，编译器将已经存在。在本秘籍中，我们将假定您已经安装了编译器。

操作步骤

我们可以使用以下任何一种方法来安装 Cython：

通过执行以下步骤从源存档中安装 Cython ：
- 下载源归档文件。
- 打开包装。
- 使用cd命令浏览到目录。
- 运行以下命令：
```
$ python setup.py install
```
使用以下任一命令从 PyPI 存储库安装 Cython：
```
$ easy_install cython
$ sudo pip install cython
```
使用非官方 Windows 安装程序，在 Windows 上安装 Cython。

另见

相关的 Cython 在线文档

构建 HelloWorld 程序

与编程语言的传统一样，我们将从 HelloWorld 示例开始。与 Python 不同，我们需要编译 Cython 代码。我们从.pyx文件开始，然后从中生成 C 代码。可以编译此.c文件，然后将其导入 Python 程序中。

操作步骤

本节介绍如何构建 Cython HelloWorld 程序：

首先，编写一些非常简单的代码以显示Hello World。这只是普通的 Python 代码，但文件具有pyx扩展名：
```
def say_hello():
  print "Hello World!"
```

创建一个名为setup.py的文件来帮助构建 Cython 代码：

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [Extension("hello", ["hello.pyx"])]

setup(
        name = 'Hello world app',
        cmdclass = {'build_ext': build_ext},
        ext_modules = ext_modules
     )

如您所见，我们在上一步中指定了文件，并为应用命名。

使用以下命令进行构建：
```
$ python setup.py build_ext --inplace
```
这将生成 C 代码，将其编译为您的平台，并产生以下输出：
```
running build_ext
cythoning hello.pyx to hello.c
building 'hello' extension
creating build
```
现在，我们可以使用以下语句导入模块：
```
from hello import say_hello
```

工作原理

在此秘籍中，我们创建了一个传统的 HelloWorld 示例。 Cython 是一种编译语言，因此我们需要编译代码。我们编写了一个包含Hello World代码的.pyx文件和一个用于生成和构建 C 代码的setup.py文件。

另见

相关的 Cython 的在线文档

将 Cython 与 NumPy 结合使用

我们可以集成 Cython 和 NumPy 代码，就像可以集成 Cython 和 Python 代码一样。让我们来看一个示例，该示例分析股票的上涨天数（股票收盘价高于前一日的天数）的比率。我们将应用二项式比值置信度的公式。您可以参考这里了解更多信息。以下公式说明该比率的重要性：

式中，p是概率，n是观察数。

操作步骤

本节通过以下步骤介绍如何将 Cython 与 NumPy 结合使用：

编写一个.pyx文件，其中包含一个函数，该函数可计算上升天数的比率和相关的置信度。首先，此函数计算价格之间的差异。然后，它计算出正差的数量，从而得出上升天数的比率。最后，在引言中的维基百科页面上应用置信度公式：
```
import numpy as np

def pos_confidence(numbers):
   diffs = np.diff(numbers)
   n = float(len(diffs))
   p = len(diffs[diffs > 0])/n
   confidence = np.sqrt(p * (1 - p)/ n)

   return (p, confidence)
```

使用上一个示例中的setup.py文件作为模板。更改明显的内容，例如.pyx文件的名称：

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [Extension("binomial_proportion", ["binomial_proportion.pyx"])]

setup(
        name = 'Binomial proportion app',
        cmdclass = {'build_ext': build_ext},
        ext_modules = ext_modules
     )

我们现在可以建立；有关更多详细信息，请参见前面的秘籍。

构建后，通过导入使用上一步中的 Cython 模块。我们将编写一个 Python 程序，使用matplotlib下载股价数据。然后我们将confidence()函数应用于收盘价：

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy
import sys
from binomial_proportion import pos_confidence

#1\. Get close prices.
today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo(sys.argv[1], start, today)
close =  numpy.array([q[4] for q in quotes])
print pos_confidence(close)

AAPL程序的输出如下：

(0.56746031746031744, 0.031209043355655924)

工作原理

我们计算了APL股价上涨的可能性和相应的置信度。我们通过创建 Cython 模块，将 NumPy 代码放入.pyx文件中，并按照上一教程中的步骤进行构建。最后，我们导入并使用了 Cython 模块。

另见

有关 Cython 的在线文档

调用 C 函数

我们可以从 Cython 调用 C 函数。在此示例中，我们调用 C log()函数。此函数仅适用于单个数字。请记住，NumPy log()函数也可以与数组一起使用。我们将计算股票价格的所谓对数回报。

操作步骤

我们首先编写一些 Cython 代码：

首先，从libc命名空间导入 C 的log()函数。然后，将此函数应用于for循环中的数字。最后，使用 NumPy diff()函数在第二步中获取对数值之间的一阶差：
```
from libc.math cimport log
import numpy as np

def logrets(numbers):
   logs = [log(x) for x in numbers] 
   return np.diff(logs)
```
先前的秘籍中已经介绍了架构。我们只需要更改setup.py文件中的某些值。

再次使用 matplotlib 下载股价数据。应用您刚刚在价格上创建的 Cython logrets()函数并绘制结果：

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
from log_returns import logrets
import matplotlib.pyplot as plt

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('AAPL', start, today)
close =  np.array([q[4] for q in quotes])
plt.plot(logrets(close))
plt.title('Logreturns of AAPL for the previous year')
plt.xlabel('Days')
plt.ylabel('Log returns')
plt.grid()
plt.show()

的AAPL对数回报结果图类似于以下屏幕截图所示：

工作原理

我们从 Cython 代码中调用了 C log()函数。该函数与 NumPy 函数一起用于计算股票的对数收益。这样，我们可以创建自己的包含便利函数的专用 API。令人高兴的是，我们的代码应该或多或少地像 Python 代码一样，以与 C 代码差不多的速度执行。

另见

相关的 Cython 在线文档位于

分析 Cython 代码

我们将使用以下公式对 Cython 和 NumPy 代码进行剖析，这些代码试图近似于欧拉常数：

有关更多背景信息，请参见。

操作步骤

本节演示如何通过以下步骤来分析 Cython 代码：

对于e的 NumPy 近似值，请按照下列步骤操作：

首先，我们将创建一个1到n的数组（在我们的示例中n是40）。

然后，我们将计算该数组的累积乘积，该乘积恰好是阶乘。在那之后，我们采取阶乘的倒数。最后，我们从维基百科页面应用公式。最后，我们放入标准配置代码，为我们提供以下程序：

from __future__ import print_function
import numpy as np
import cProfile
import pstats

def approx_e(n=40, display=False):
   # array of [1, 2, ... n-1]
   arr = np.arange(1, n)

   # calculate the factorials and convert to floats
   arr = arr.cumprod().astype(float)

   # reciprocal 1/n
   arr = np.reciprocal(arr)

   if display:
    print(1 + arr.sum())

# Repeat multiple times because NumPy is so fast
def run(repeat=2000):
    for i in range(repeat):
        approx_e()

cProfile.runctx("run()", globals(), locals(), "Profile.prof")

s = pstats.Stats("Profile.prof")
s.strip_dirs().sort_stats("time").print_stats()

approx_e(display=True)

以下代码段显示了e的近似值的分析输出和结果。有关概要分析输出的更多信息，请参见第 7 章，“分析和调试”。

 8004 function calls in 0.016 seconds

 Ordered by: internal time

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 2000    0.007    0.000    0.015    0.000 numpy_approxe.py:6(approx_e)
 2000    0.004    0.000    0.004    0.000 {method 'cumprod' of 'numpy.ndarray' objects}
 2000    0.002    0.000    0.002    0.000 {numpy.core.multiarray.arange}
 2000    0.002    0.000    0.002    0.000 {method 'astype' of 'numpy.ndarray' objects}
 1    0.001    0.001    0.016    0.016 numpy_approxe.py:20(run)
 1    0.000    0.000    0.000    0.000 {range}
 1    0.000    0.000    0.016    0.016 <string>:1(<module>)
 1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

2.71828182846

Cython 代码使用与上一步所示相同的算法，但是实现方式不同。便利函数较少，实际上我们现在需要一个for循环。另外，我们需要为某些变量指定类型。 .pyx文件的代码如下所示：

def approx_e(int n=40, display=False):
    cdef double sum = 0.
    cdef double factorial = 1.
    cdef int k

    for k in xrange(1,n+1):
        factorial *= k
        sum += 1/factorial

    if display:
        print(1 + sum)

以下 Python 程序导入 Cython 模块并进行一些分析：

import pstats
import cProfile
import pyximport
pyximport.install()

import approxe

# Repeat multiple times because Cython is so fast
def run(repeat=2000):
    for i in range(repeat):
        approxe.approx_e()

cProfile.runctx("run()", globals(), locals(), "Profile.prof")

s = pstats.Stats("Profile.prof")
s.strip_dirs().sort_stats("time").print_stats()

approxe.approx_e(display=True)

这是 Cython 代码的分析输出：

 2004 function calls in 0.001 seconds

 Ordered by: internal time

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 2000    0.001    0.000    0.001    0.000 {approxe.approx_e}
 1    0.000    0.000    0.001    0.001 cython_profile.py:9(run)
 1    0.000    0.000    0.000    0.000 {range}
 1    0.000    0.000    0.001    0.001 <string>:1(<module>)
 1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

2.71828182846

工作原理

我们分析了 NumPy 和 Cython 代码。 NumPy 已针对速度进行了优化，因此 NumPy 和 Cython 程序都是高性能程序，我们对此不会感到惊讶。但是，当比较 2,000 次近似代码的总时间时，我们意识到 NumPy 需要 0.016 秒，而 Cython 仅需要 0.001 秒。显然，实际时间取决于硬件，操作系统和其他因素，例如计算机上运行的其他进程。同样，提速取决于代码类型，但我希望您同意，根据经验，Cython 代码会更快。

另见

相关的 Cython 在线文档

Cython 的近似阶乘

最后一个示例是 Cython 的近似阶乘。我们将使用两种近似方法。首先，我们将应用斯特林近似方法。斯特林近似的公式如下：

其次，我们将使用 Ramanujan 的近似值，并使用以下公式：

操作步骤

本节介绍如何使用 Cython 近似阶乘。您可能还记得，在本秘籍中，我们使用在 Cython 中可选的类型。从理论上讲，声明静态类型应加快速度。静态类型化提供了一些有趣的挑战，这些挑战在编写 Python 代码时可能不会遇到，但请不要担心。我们将尝试使其简单：

除了将函数参数和一个局部变量声明为ndarray数组外，我们将编写的 Cython 代码类似于常规的 Python 代码。为了使静态类型起作用，我们需要cimport NumPy。另外，我们必须使用cdef关键字声明局部变量的类型：

import numpy
cimport numpy

def ramanujan_factorial(numpy.ndarray n):
   sqrt_pi = numpy.sqrt(numpy.pi, dtype=numpy.float64)
   cdef numpy.ndarray root = (8 * n + 4) * n + 1 
   root = root * n + 1/30.
   root = root ** (1/6.)
   return sqrt_pi * calc_eton(n) * root

def stirling_factorial(numpy.ndarray n):
    return numpy.sqrt(2 * numpy.pi * n) * calc_eton(n)

def calc_eton(numpy.ndarray n):
    return (n/numpy.e) ** n

如先前的教程所示，构建需要我们创建一个setup.py文件，但是现在我们通过调用get_include()函数来包含与 NumPy 相关的目录。通过此修订，setup.py文件具有以下内容：

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy

ext_modules = [Extension("factorial", ["factorial.pyx"], include_dirs = [numpy.get_include()])] 

setup(
        name = 'Factorial app',
        cmdclass = {'build_ext': build_ext},
        ext_modules = ext_modules
     )

绘制两种近似方法的相对误差。像我们在整本书中所做的那样，将使用 NumPy cumprod()函数计算相对于阶乘值的误差：

from factorial import ramanujan_factorial
from factorial import stirling_factorial
import numpy as np
import matplotlib.pyplot as plt

N = 50
numbers = np.arange(1, N)
factorials = np.cumprod(numbers, dtype=float)

def error(approximations):
   return (factorials - approximations)/factorials

plt.plot(error(ramanujan_factorial(numbers)), 'b-', label='Ramanujan')
plt.plot(error(stirling_factorial(numbers)), 'ro', label='Stirling')
plt.title('Factorial approximation relative errors')
plt.xlabel('n')
plt.ylabel('Relative error')
plt.grid()
plt.legend(loc='best')
plt.show()

下图显示了 Ramanujan 近似值（点）和 Stirling 近似值（线）的相对误差：

工作原理

在此示例中，我们看到了 Cython 静态类型的演示。该秘籍的主要成分如下：

cimport，它导入 C 声明
包含具有get_include()函数的目录
cdef关键字，用于定义局部变量的类型

另见

相关的 Cython 在线文档

十、Scikits 的乐趣

在本章中，我们将介绍以下秘籍：

安装 scikit-learn
加载示例数据集
用 scikit-learn 对道琼斯股票进行聚类
安装 Statsmodels
使用 Statsmodels 执行正态性检验
安装 scikit-image
检测角点
检测边界
安装 Pandas
使用 Pandas 估计股票收益的相关性
从 Statsmodels 中将数据作为 pandas 对象加载
重采样时间序列数据

简介

Scikits 是小型的独立项目，以某种方式与 SciPy 相关，但不属于 SciPy。这些项目不是完全独立的，而是作为一个联合体在伞下运行的。在本章中，我们将讨论几个 Scikits 项目，例如：

scikit-learn，机器学习包
Statsmodels，统计数据包
scikit-image，图像处理包
Pandas，数据分析包

安装 scikit-learn

scikit-learn 项目旨在提供一种用于机器学习的 API。我最喜欢的是令人惊叹的文档。我们可以使用操作系统的包管理器安装 scikit-learn。根据操作系统的不同，此选项可能可用也可能不可用，但它应该是最方便的方法。

Windows 用户只需从项目网站下载安装程序即可。在 Debian 和 Ubuntu 上，该项目称为python-sklearn。在 MacPorts 上，这些端口称为py26-scikits-learn和py27-scikits-learn。我们也可以从源代码或使用easy_install安装。 PythonXY，Enthought 和 NetBSD 提供了第三方发行版。

准备

您需要安装 SciPy 和 NumPy。返回第 1 章，“使用 IPython”，以获取必要的说明。

操作步骤

现在让我们看一下如何安装 scikit-learn 项目：

使用easy_install进行安装：在命令行中键入以下命令之一：
```
$ pip install -U scikit-learn
$ easy_install -U scikit-learn
```
由于权限问题，这个可能无法工作，因此您可能需要在命令前面编写sudo，或以管理员身份登录。
从源代码安装：下载源代码，解压缩并使用cd进入下载的文件夹。然后发出以下命令：
```
$ python setup.py install
```

加载示例数据集

scikit-learn 项目附带了许多我们可以尝试的数据集和样例图像。在本秘籍中，我们将加载 scikit-learn 分发中包含的示例数据集。数据集将数据保存为 NumPy 二维数组，并将元数据链接到该数据。

操作步骤

我们将加载波士顿房价样本数据集。这是一个很小的数据集，因此，如果您要在波士顿寻找房子，请不要太兴奋！其他数据集在这个页面中进行了描述。

我们将查看原始数据的形状及其最大值和最小值。形状是一个元组，表示 NumPy 数组的大小。我们将对目标数组执行相同的操作，其中包含作为学习目标（确定房价）的值。以下sample_data.py 中的代码实现了我们的目标：

from __future__ import print_function
from sklearn import datasets

boston_prices = datasets.load_boston()
print("Data shape", boston_prices.data.shape)
print("Data max=%s min=%s" % (boston_prices.data.max(), boston_prices.data.min()))
print("Target shape", boston_prices.target.shape)
print("Target max=%s min=%s" % (boston_prices.target.max(), boston_prices.target.min()))

我们计划的结果如下：

Data shape (506, 13)
Data max=711.0 min=0.0
Target shape (506,)
Target max=50.0 min=5.0

利用 scikits-learn 对道琼斯股票进行聚类

聚类是一种机器学习算法，旨在基于相似度对项目进行分组。在此示例中，我们将使用道琼斯工业平均指数（DJI 或 DJIA）进行聚类。本秘籍的大多数步骤已通过前面各章的审查。

操作步骤

首先，我们将从 Yahoo 金融下载这些股票的 EOD 价格数据。然后，我们将计算平方亲和矩阵。最后，我们将股票与AffinityPropagation类聚类：

使用 DJI 指数的股票代码下载 2011 年价格数据。在此示例中，我们仅对收盘价感兴趣：

# 2011 to 2012
start = datetime.datetime(2011, 01, 01)
end = datetime.datetime(2012, 01, 01)

#Dow Jones symbols
symbols = ["AA", "AXP", "BA", "BAC", "CAT",
   "CSCO", "CVX", "DD", "DIS", "GE", "HD",
   "HPQ", "IBM", "INTC", "JNJ", "JPM", 
   "KO", "MCD", "MMM", "MRK", "MSFT", "PFE",
   "PG", "T", "TRV", "UTX", "VZ", "WMT", "XOM"]

quotes = []

for symbol in symbols:
    try:
        quotes.append(finance.quotes_historical_yahoo(symbol, start, end, asobject=True))
    except urllib2.HTTPError as e:
        print(symbol, e)

close = np.array([q.close for q in quotes]).astype(np.float)
print(close.shape)

使用对数收益率作为指标来计算不同股票之间的相似度。我们正在要做的是计算数据点的欧几里德距离：

logreturns = np.diff(np.log(close))
print(logreturns.shape)

logreturns_norms = np.sum(logreturns ** 2, axis=1)
S = - logreturns_norms[:, np.newaxis] - logreturns_norms[np.newaxis, :] + 2 * np.dot(logreturns, logreturns.T)

给AffinityPropagation类上一步的结果。此类为数据点或在我们的情况下的股票标记了适当的群集编号：

aff_pro = sklearn.cluster.AffinityPropagation().fit(S)
labels = aff_pro.labels_

for symbol, label in zip(symbols, labels):
    print('%s in Cluster %d' % (symbol, label))

完整的群集程序如下：

from __future__ import print_function
import datetime
import numpy as np
import sklearn.cluster
from matplotlib import finance
import urllib2

#1\. Download price data

# 2011 to 2012
start = datetime.datetime(2011, 01, 01)
end = datetime.datetime(2012, 01, 01)

#Dow Jones symbols
symbols = ["AA", "AXP", "BA", "BAC", "CAT",
   "CSCO", "CVX", "DD", "DIS", "GE", "HD",
   "HPQ", "IBM", "INTC", "JNJ", "JPM", 
   "KO", "MCD", "MMM", "MRK", "MSFT", "PFE",
   "PG", "T", "TRV", "UTX", "VZ", "WMT", "XOM"]

quotes = []

for symbol in symbols:
    try:
        quotes.append(finance.quotes_historical_yahoo(symbol, start, end, asobject=True))
    except urllib2.HTTPError as e:
        print(symbol, e)

close = np.array([q.close for q in quotes]).astype(np.float)
print(close.shape)

#2\. Calculate affinity matrix
logreturns = np.diff(np.log(close))
print(logreturns.shape)

logreturns_norms = np.sum(logreturns ** 2, axis=1)
S = - logreturns_norms[:, np.newaxis] - logreturns_norms[np.newaxis, :] + 2 * np.dot(logreturns, logreturns.T)

#3\. Cluster using affinity propagation
aff_pro = sklearn.cluster.AffinityPropagation().fit(S)
labels = aff_pro.labels_

for symbol, label in zip(symbols, labels):
    print('%s in Cluster %d' % (symbol, label))

具有每个股票的簇号的输出如下：

(29, 252)
(29, 251)
AA in Cluster 0
AXP in Cluster 6
BA in Cluster 6
BAC in Cluster 1
CAT in Cluster 6
CSCO in Cluster 2
CVX in Cluster 7
DD in Cluster 6
DIS in Cluster 6
GE in Cluster 6
HD in Cluster 5
HPQ in Cluster 3
IBM in Cluster 5
INTC in Cluster 6
JNJ in Cluster 5
JPM in Cluster 4
KO in Cluster 5
MCD in Cluster 5
MMM in Cluster 6
MRK in Cluster 5
MSFT in Cluster 5
PFE in Cluster 7
PG in Cluster 5
T in Cluster 5
TRV in Cluster 5
UTX in Cluster 6
VZ in Cluster 5
WMT in Cluster 5
XOM in Cluster 7

工作原理

下表是在秘籍中使用的函数的概述：

函数	描述
`sklearn.cluster.AffinityPropagation()`	创建一个`AffinityPropagation`对象。
`sklearn.cluster.AffinityPropagation.fit()`	从欧几里得距离计算亲和度矩阵，并应用亲和度传播聚类。
`diff()`	计算 NumPy 数组中数字的差。如果未指定，则计算一阶差。
`log()`	计算 NumPy 数组中元素的自然对数。
`sum()`	对 NumPy 数组的元素求和。
`dot()`	这对二维数组执行矩阵乘法。它还计算一维数组的内积。

另见

相关文档

安装 Statsmodels

statsmodels包专注于统计建模。我们可以将其与 NumPy 和 pandas 集成（在本章稍后的内容中将有更多关于 pandas 的信息）。

操作步骤

可以从这里下载源码和二进制文件。如果要从源代码安装，请运行以下命令：

$ python setup.py install

如果使用setuptools，则命令如下：

$ easy_install statsmodels

使用 Statsmodels 执行正态性检验

statsmodels包具有许多统计检验。我们将看到这样一个检验的示例 -- Anderson-Darling 检验用于正态性。

操作步骤

我们将像以前的秘籍一样下载价格数据，但这一次是单只股票。再次，我们将计算该股票收盘价的对数收益，并将其用作正态性检验函数的输入。

此函数返回一个包含第二个元素的元组，即 p 值，介于 0 和 1 之间。本教程的完整代码如下：

from __future__ import print_function
import datetime
import numpy as np
from matplotlib import finance
from statsmodels.stats.adnorm import normal_ad

#1\. Download price data

# 2011 to 2012
start = datetime.datetime(2011, 01, 01)
end = datetime.datetime(2012, 01, 01)

quotes = finance.quotes_historical_yahoo('AAPL', start, end, asobject=True)

close = np.array(quotes.close).astype(np.float)
print(close.shape)

print(normal_ad(np.diff(np.log(close))))

#Retrieving data for AAPL
#(252,)
#(0.57103805516803163, 0.13725944999430437)

下面显示了0.13为p-value的脚本的输出：

Retrieving data for AAPL
(252,)
(0.57103805516803163, 0.13725944999430437)

工作原理

如statsmodels中所示，此秘籍证明了 Anderson-Darling 统计检验的正态性。我们使用没有正态分布的股票价格数据作为输入。对于数据，我们获得了0.13的 p 值。由于概率在 0 到 1 之间，这证实了我们的假设。

安装 scikit-image

scikit-image 是用于图像处理的工具包，它依赖 PIL，SciPy，Cython 和 NumPy。 Windows 安装程序也可用。该工具包是 Enthought Python 发行版以及 PythonXY 发行版的一部分。

操作步骤

与往常一样，使用以下两个命令之一安装 scikit-image：

$ pip install -U scikit-image
$ easy_install -U scikit-image

同样，您可能需要以超级用户身份运行这些命令。

另一种选择是通过克隆 Git 存储库或从 Github 下载该存储库作为源归档来获取最新的开发版本。然后运行以下命令：

$ python setup.py install

检测角点

角点检测是计算机视觉中的标准技术。 scikit-image 提供了 Harris 角点检测器，它很棒，因为角检测非常复杂。显然，我们可以从头开始做，但是这违反了不重新发明轮子的基本原则。

准备

您可能需要在系统上安装jpeglib，才能加载 scikit-learn 图像（是 JPEG 文件）。如果您使用的是 Windows，请使用安装程序。否则，下载发行版，解压缩它，并使用以下命令从顶部文件夹中进行构建：

$ ./configure
$  make
$  sudo make install

操作步骤

我们将从 scikit-learn 加载示例图像。对于此示例，这不是绝对必要的；您可以改用其他任何图像：

scikit-learn 当前在数据集结构中有两个样例 JPEG 图像。只看第一张图片：
```
dataset = load_sample_images()
img = dataset.images[0]
```
自本书第一版以来，API 发生了变化。例如，对于 scikit-image 0.11.2，我们需要首先将彩色图像的值转换为灰度值。图像灰度如下：
```
gray_img = rgb2gray(img)
```

调用corner_harris()函数获取角点的坐标：

harris_coords = corner_peaks(corner_harris(gray_img))
y, x = np.transpose(harris_coords)

拐角检测的代码如下：

from sklearn.datasets import load_sample_images
import matplotlib.pyplot as plt
import numpy as np
from skimage.feature import corner_harris
from skimage.feature import corner_peaks
from skimage.color import rgb2gray

dataset = load_sample_images()
img = dataset.images[0]
gray_img = rgb2gray(img)
harris_coords = corner_peaks(corner_harris(gray_img))
y, x = np.transpose(harris_coords)
plt.axis('off')
plt.imshow(img)
plt.plot(x, y, 'ro')
plt.show()

我们得到一个带有点的图像，脚本在其中检测角点，如以下屏幕截图所示：

工作原理

我们对 scikit-image 的样例图像应用了 Harris 角点检测。如您所见，结果非常好。我们只能使用 NumPy 做到这一点，因为它只是一个简单的线性代数类型的计算。仍然，可能会变得凌乱。 scikit-image 工具包具有更多类似的功能，因此，如果需要图像处理例程，请查看 scikit-image 文档。另外请记住，API 可能会发生快速变化。

另见

相关的 scikit-image 文档

检测边界

边界检测是另一种流行的图像处理技术。 scikit-image 具有基于高斯分布的标准差的 Canny 过滤器实现，可以开箱即用地执行边界检测。除了将图像数据作为 2D 数组外，此过滤器还接受以下参数：

高斯分布的标准差
下限阈值
上限阈值

操作步骤

我们将使用与先前秘籍相同的图像。代码几乎相同（请参见edge_detection.py）。请特别注意我们称为 Canny 筛选器函数的行：

from sklearn.datasets import load_sample_images
import matplotlib.pyplot as plt
import skimage.feature

dataset = load_sample_images()
img = dataset.images[0] 
edges = skimage.feature.canny(img[..., 0])
plt.axis('off')
plt.imshow(edges)
plt.show()

该代码生成原始图像中边界的图像，如以下屏幕截图所示：

另见

相关文档

安装 Pandas

Pandas 是用于数据分析的 Python 库。它与 R 编程语言有一些相似之处，并非偶然。 R 是一种受数据科学家欢迎的专业编程语言。例如，R 启发了 Pandas 的核心DataFrame对象。

操作步骤

在 PyPi 上，该项目称为pandas。因此，您可以运行以下命令之一：

$ sudo easy_install -U pandas
$ pip install pandas

如果使用 Linux 包管理器，则需要安装python-pandas项目。在 Ubuntu 上，执行以下操作：

$ sudo apt-get install python-pandas

您也可以从源代码安装（除非下载源代码存档，否则需要 Git）：

$ git clone git://github.com/pydata/pandas.git
$ cd pandas
$ python setup.py install

另见

相关文档

使用 Pandas 估计股票收益的相关性

Pandas DataFrame是类似矩阵和字典的数据结构，类似于 R 中提供的功能。实际上，它是 Pandas 的中心数据结构，您可以应用各种操作。例如，查看投资组合的相关矩阵是很常见的，所以让我们开始吧。

操作步骤

首先，我们将为每个符号的每日对数回报创建带有 Pandas 的DataFrame。然后，我们将在约会中加入这些。最后，将打印相关性，并显示一个图：

要创建数据框，请创建一个包含股票代码作为键的字典，并将相应的日志作为值返回。数据框本身以日期作为索引，将股票代码作为列标签：

data = {}

for i, symbol in enumerate(symbols):
   data[symbol] = np.diff(np.log(close[i]))

# Convention: import pandas as pd
df = pd.DataFrame(data, index=dates[0][:-1], columns=symbols)

现在，我们可以执行诸如计算相关矩阵或在数据帧上绘制等操作：

print(df.corr())
df.plot()

完整的源代码（也可以下载价格数据）如下：

from __future__ import print_function
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime
from matplotlib import finance
import numpy as np

# 2011 to 2012
start = datetime(2011, 01, 01)
end = datetime(2012, 01, 01)

symbols = ["AA", "AXP", "BA", "BAC", "CAT"]

quotes = [finance.quotes_historical_yahoo(symbol, start, end, asobject=True)
          for symbol in symbols]

close = np.array([q.close for q in quotes]).astype(np.float)
dates = np.array([q.date for q in quotes])

data = {}

for i, symbol in enumerate(symbols):
   data[symbol] = np.diff(np.log(close[i]))

df = pd.DataFrame(data, index=dates[0][:-1], columns=symbols)

print(df.corr())
df.plot()
plt.legend(symbols)
plt.show()

#           AA       AXP        BA       BAC       CAT
#AA   1.000000  0.768484  0.758264  0.737625  0.837643
#AXP  0.768484  1.000000  0.746898  0.760043  0.736337
#BA   0.758264  0.746898  1.000000  0.657075  0.770696
#BAC  0.737625  0.760043  0.657075  1.000000  0.657113
#CAT  0.837643  0.736337  0.770696  0.657113  1.000000

这是相关矩阵的输出：

 AA       AXP        BA       BAC       CAT
AA   1.000000  0.768484  0.758264  0.737625  0.837643
AXP  0.768484  1.000000  0.746898  0.760043  0.736337
BA   0.758264  0.746898  1.000000  0.657075  0.770696
BAC  0.737625  0.760043  0.657075  1.000000  0.657113
CAT  0.837643  0.736337  0.770696  0.657113  1.000000

显示了五只股票的对数收益图：

工作原理

我们使用了以下DataFrame方法：

函数	描述
`pandas.DataFrame()`	此函数使用指定的数据，索引（行）和列标签构造`DataFrame`。
`pandas.DataFrame.corr()`	该函数计算列的成对相关，而忽略缺失值。默认情况下，使用 Pearson 相关。
`pandas.DataFrame.plot()`	此函数使用`matplotlib`绘制数据帧。

另见

相关文档
第 4 章，“Pandas 入门书”，摘自 Ivan Idris 的书“Python 数据分析”， Packt Publishing

从 Statsmodels 中将数据作为 pandas 对象加载

statsmodels 在其发行版中有很多样本数据集。完整列表可以在这个页面中找到。

在本教程中，我们将专注于铜数据集，其中包含有关铜价，世界消费量和其他参数的信息。

准备

在开始之前，我们可能需要安装 patsy。 patsy 是描述统计模型的库。很容易看出这个库是否是必需的。只需运行代码。如果您收到与 patsy 相关的错误，请执行以下任一命令：

$ sudo easy_install patsy
$ pip install --upgrade patsy

操作步骤

在本节中，我们将从 statsmodels 中加载数据集作为 Pandas DataFrame或Series对象。

我们需要调用的函数是load_pandas()。如下加载数据：
```
data = statsmodels.api.datasets.copper.load_pandas()
```
这会将数据加载到包含 Pandas 对象的DataSet对象中。
DataSet对象具有名为exog的属性，当作为 Pandas 对象加载时，该属性将成为具有多个列的DataFrame对象。在我们的案例中，它还有一个endog属性，其中包含世界铜消费量的值。

通过创建OLS对象并调用其fit()方法来执行普通的最小二乘计算，如下所示：
```
x, y = data.exog, data.endog

fit = statsmodels.api.OLS(y, x).fit()
print("Fit params", fit.params)
```
这应该打印拟合过程的结果：
```
Fit params COPPERPRICE         14.222028
INCOMEINDEX       1693.166242
ALUMPRICE          -60.638117
INVENTORYINDEX    2515.374903
TIME               183.193035
```

OLS 拟合的结果可以通过summary()方法总结如下：

print(fit.summary())

这将为我们提供以下回归结果输出：

加载铜数据集所需的代码如下：

from __future__ import print_function
import statsmodels.api

# See https://github.com/statsmodels/statsmodels/tree/master/statsmodels/datasets
data = statsmodels.api.datasets.copper.load_pandas()

x, y = data.exog, data.endog

fit = statsmodels.api.OLS(y, x).fit()
print("Fit params", fit.params)
print()
print("Summary")
print()
print(fit.summary())

工作原理

Statsmodels 的Dataset类中的数据遵循特殊格式。其中，此类具有endog和exog属性。 Statsmodels 具有load()函数，该函数将数据作为 NumPy 数组加载。相反，我们使用了load_pandas()方法，该方法将数据加载为pandas对象。我们进行了 OLS 拟合，基本上为我们提供了铜价和消费量的统计模型。

另见

相关文档

重采样时间序列数据

在此教程中，您将学习如何使用 Pandas 对时间序列进行重新采样。

操作步骤

我们将下载AAPL的每日价格时间序列数据，然后通过计算平均值将其重新采样为每月数据。我们将通过创建 Pandas DataFrame并调用其resample() 方法来做到这一点：

在创建 Pandas DataFrame之前，我们需要创建一个DatetimeIndex对象传递给DataFrame构造器。根据下载的报价数据创建索引，如下所示：
```
dt_idx = pandas.DatetimeIndex(quotes.date)
```
获得日期时间索引后，我们将其与收盘价一起使用以创建数据框：
```
df = pandas.DataFrame (quotes.close, index=dt_idx, columns=[symbol])
```

通过计算平均值，将时间序列重新采样为每月频率：

resampled = df.resample('M', how=numpy.mean)
print(resampled)

如下行所示，重新采样的时间序列每个月都有一个值：

 AAPL
2011-01-31  336.932500
2011-02-28  349.680526
2011-03-31  346.005652
2011-04-30  338.960000
2011-05-31  340.324286
2011-06-30  329.664545
2011-07-31  370.647000
2011-08-31  375.151304
2011-09-30  390.816190
2011-10-31  395.532381
2011-11-30  383.170476
2011-12-31  391.251429

使用DataFrame plot()方法绘制数据：

df.plot()
resampled.plot()
plt.show()

原始时间序列的图如下：

重采样的数据具有较少的数据点，因此，生成的图更加混乱，如以下屏幕截图所示：

完整的重采样代码如下：

from __future__ import print_function
import pandas
import matplotlib.pyplot as plt
from datetime import datetime
from matplotlib import finance
import numpy as np

# Download AAPL data for 2011 to 2012
start = datetime(2011, 01, 01)
end = datetime(2012, 01, 01)

symbol = "AAPL"
quotes = finance.quotes_historical_yahoo(symbol, start, end, asobject=True)

# Create date time index
dt_idx = pandas.DatetimeIndex(quotes.date)

#Create data frame
df = pandas.DataFrame(quotes.close, index=dt_idx, columns=[symbol])

# Resample with monthly frequency
resampled = df.resample('M', how=np.mean)
print(resampled) 

# Plot
df.plot()
plt.title('AAPL prices')
plt.ylabel('Price')

resampled.plot()
plt.title('Monthly resampling')
plt.ylabel('Price')
plt.grid(True)
plt.show()

工作原理

我们根据日期和时间列表创建了日期时间索引。然后，该索引用于创建 Pandas DataFrame。然后，我们对时间序列数据进行了重新采样。单个字符给出重采样频率，如下所示：

每天D
每月M
每年A

resample()方法的how参数指示如何采样数据。默认为计算平均值。

另见

相关 Pandas 文档

十一、最新最强的 NumPy

在本章中，我们涵盖以下秘籍：

用at()方法用花式索引代替 ufuncs
通过使用partition()函数选择快速中位数进行部分排序
使用nanmean()，nanvar()和nanstd()函数跳过 NaN
使用full()和full_like()函数创建值初始化的数组
numpy.random.choice()随机抽样
使用datetime64类型和相关的 API

简介

自《NumPy 秘籍》第一版以来，NumPy 团队引入了新功能；我将在本章中对其进行描述。您可能不太可能阅读本书的第一版，而现在正在阅读第二版。我在 2012 年撰写了第一版，并使用了当时可用的功能。 NumPy 具有许多功能，因此您不能期望涵盖所有功能，但是我在本章中介绍的功能相对重要。

使用`at()`方法为 ufuncs 建立花式索引

at()方法已添加到 NumPy 1.8 的 NumPy 通用函数类中。此方法允许就地进行花式索引。花式索引是不涉及整数或切片的索引，这是正常的索引。 “就地”是指将更改输入数组的数据。

at()方法的签名为ufunc.at(a, indices[, b])。索引数组对应于要操作的元素。我们仅必须为具有两个操作数的通用函数指定b数组。

操作步骤

以下步骤演示了at()方法的工作方式：

创建一个具有种子44的7个从-4到4的随机整数的数组。

np.random.seed(44)
a = np.random.random_integers(-4, 4, 7)
print(a)

该数组如下所示：

[ 0 -1 -3 -1 -4  0 -1]

将sign通用函数的at()方法应用于第三和第五个数组元素：
```
np.sign.at(a, [2, 4])
print(a)
```
我们得到以下更改后的数组：
```
[ 0 -1 -1 -1 -1  0 -1]
```

另见

NumPy 通用函数文档

使用`partition()`函数通过快速中位数的选择进行部分排序

partition()子例程进行部分排序。这应该比正常的分类工作少。

注意

有关更多信息，请参见这里。有用的情况是选择组中的前五项（或其他一些数字）。部分排序不能在顶部元素集中保留正确的顺序。

子例程的第一个参数是要排序的输入数组。第二个参数是整数或与数组元素的索引相对应的整数列表。 partition()子例程正确地对那些索引处的项目进行排序。一个指定的索引给出两个分区。多个索自举致两个以上的分区。该算法保证分区中小于正确排序项目的项目位于该项目之前。否则，它们将放在该项目的后面。

操作步骤

让我们举例说明此解释：

创建一个具有随机数的数组以进行排序：

np.random.seed(20)
a = np.random.random_integers(0, 7, 9)
print(a)

该数组具有以下元素：

[3 2 7 7 4 2 1 4 3]

通过将数组划分为两个大致相等的部分，对数组进行部分排序：
```
print(np.partition(a, 4))
```
我们得到以下结果：
```
[2 3 1 2 3 7 7 4 4]
```

工作原理

我们对 9 个元素的数组进行了部分排序。该函数保证索引4,的中间只有一个元素在正确的位置。这对应于尝试选择数组的前五项而不关心前五组中的顺序。由于正确排序的项目位于中间，因此这也将返回数组的中位数。

另见

相关 NumPy 文档

使用`nanmean()`，`nanvar()`和`nanstd()`函数跳过 NaN

试图估计一组数据的算术平均值，方差和标准差是很常见的。

一种简单但有效的方法称为 Jackknife 重采样。 Jackknife 重采样的想法是通过每次都遗漏一个值来从原始数据创建数据集。本质上，我们试图估计如果至少一个值不正确会发生什么。对于每个新数据集，我们都会重新计算我们感兴趣的统计估计量。这有助于我们了解估计量的变化方式。

操作步骤

我们将折刀重采样应用于随机数据。通过将其设置为 NaN（非数字），我们将跳过每个数组元素一次。然后，可以使用nanmean()，nanvar()和nanstd()计算算术平均值，方差和标准差：

首先为估算值初始化一个30 x 3的数组，如下所示：
```
estimates = np.zeros((len(a), 3))
```
遍历数组并通过在循环的每次迭代中将一个值设置为 NaN 来创建新的数据集。对于每个新数据集，计算估计值：
```
for i in xrange(len(a)):
   b = a.copy()
   b[i] = np.nan
   estimates[i,] = [np.nanmean(b), np.nanvar(b), np.nanstd(b)]
```

打印每个估计量的方差：

print("Estimator variance", estimates.var(axis=0))

屏幕上显示以下输出：

Estimator variance [ 0.00079905  0.00090129  0.00034604]

工作原理

我们用折刀重采样估计了数据集的算术平均值，方差和标准差的方差。这表明算术平均值，方差和标准差有多少变化。该秘籍的代码在本书的代码包的jackknife.py 文件中：

from __future__ import print_function
import numpy as np

np.random.seed(46)
a = np.random.randn(30)
estimates = np.zeros((len(a), 3))

for i in xrange(len(a)):
   b = a.copy()
   b[i] = np.nan

   estimates[i,] = [np.nanmean(b), np.nanvar(b), np.nanstd(b)]

print("Estimator variance", estimates.var(axis=0))

另见

使用`full()`和`full_like()`函数创建初始化值的数组

full()和full_like()函数是 NumPy 的新增函数，旨在促进初始化。这是文档中关于它们的内容：

>>> help(np.full)
Return a new array of given shape and type, filled with `fill_value`.
>>> help(np.full_like)
Return a full array with the same shape and type as a given array.

操作步骤

让我们看一下full()和full_like()函数：

用full()创建一个1x2数组，并填充幸运数字7：
```
print(np.full((1, 2), 7))
```
因此，我们得到以下数组：
```
array([[ 7.,  7.]])
```
数组元素是浮点数。
指定整数数据类型，如下所示：
```
print(np.full((1, 2), 7, dtype=np.int))
```
输出相应地更改：
```
array([[7, 7]])
```
full_like()函数检查数组的元数据，并将其重新用于新数组。例如，使用linspace()创建一个数组，并将其用作full_like()函数的模板：
```
a = np.linspace(0, 1, 5)
print(a)
array([ 0\.  ,  0.25,  0.5 ,  0.75,  1\.  ])
print(np.full_like(a, 7))
array([ 7.,  7.,  7.,  7.,  7.])
```
同样，我们用幸运数字7填充数组。要将数据类型修改为整数，请在以下行中使用：
```
print(np.full_like(a, 7, dtype=np.int))
array([7, 7, 7, 7, 7])
```

工作原理

我们用full()和full_like()产生了数组。 full()函数用数字7填充数组。 full_like()函数重新使用了数组的元数据来创建新的数组。这两个函数都可以指定数组的数据类型。

使用`numpy.random.choice()`进行随机采样

自举的过程类似于粗加工。基本的自举方法包括以下步骤：

从大小为 N 的原始数据生成样本。将原始数据样本可视化为一碗数字。我们通过从碗中随机抽取数字来创建新样本。取一个数字后，我们将其放回碗中。
对于每个生成的样本，我们计算感兴趣的统计估计量（例如，算术平均值）。

操作步骤

我们将应用numpy.random.choice()进行自举：

按照二项式分布生成数据样本，该数据样本模拟五次抛掷公平硬币：
```
N = 400
np.random.seed(28)
data = np.random.binomial(5, .5, size=N)
```

生成 30 个样本并计算其平均值（更多样本将得到更好的结果）：

bootstrapped = np.random.choice(data, size=(N, 30))
means = bootstrapped.mean(axis=0)

使用 matplotlib 箱形图可视化算术平均值分布：

plt.title('Bootstrapping demo')
plt.grid()
plt.boxplot(means)
plt.plot(3 * [data.mean()], lw=3, label='Original mean')
plt.legend(loc='best')
plt.show()

有关最终结果，请参考以下带注释的图：

工作原理

我们模拟了一个实验，该实验涉及掷出五次公平硬币。我们通过创建样本并计算相应的方法来自举数据。然后，我们使用numpy.random.choice()进行自举。我们用matplotlib箱形图直观地表示了均值。如果您不熟悉箱形图，图中的注释将对您有所帮助。箱形图中的以下元素很重要：

中位数由框中的一条线表示。
上下四分位数显示为框的边界。
胡须指示异常值的边界。默认情况下，这些值从框的边界设置为1.5 * (Q3 - Q1)，也称为四分位间距。

另见

使用`datetime64`类型和相关的 API

datetime64类型表示相应的日期和时间。您需要 NumPy 1.7.0 或更高版本才能使用此数据类型。

操作步骤

要熟悉datetime64，请按照下列步骤操作：

从字符串创建一个datetime64，如下所示：
```
print(np.datetime64('2015-05-21'))
```
前一行输出以下输出：
```
numpy.datetime64('2015-05-21')
```
我们使用YYYY-MM-DD格式在 2015 年 5 月 21 日创建了datetime64类型，其中Y对应于年份，M对应于月份，D对应于每月的一天。 NumPy 符合 ISO 8601 标准 -- 一种表示日期和时间的国际标准。

ISO 8601 还定义了YYYY-MM-DD，YYYY-MM和YYYYMMDD格式。自己检查这些，如下所示：

print(np.datetime64('20150521'))
print(np.datetime64('2015-05'))

该代码显示以下行：

numpy.datetime64('20150521')
numpy.datetime64('2015-05')

默认情况下，ISO 8601 使用本地时区。可以使用T[hh:mm:ss]格式指定时间。例如，我们可以定义 1578 年 1 月 1 日和晚上 9:18。如下：
```
local = np.datetime64('1578-01-01T21:18')
print(local)
```
以下行显示了结果：
```
numpy.datetime64('1578-01-01T21:18Z')
```
格式为-[hh:mm]的字符串定义相对于 UTC 时区的偏移量。我们可以使用 8 个小时的偏移量创建一个datetime64类型，如下所示：
```
with_offset = np.datetime64('1578-01-01T21:18-0800')
print(with_offset)
```
然后，我们在屏幕上看到以下行：
```
numpy.datetime64('1578-01-02T05:18Z')
```
最后的Z代表 Zulu 时间，有时也称为 UTC。
相互减去两个datetime64对象：
```
print(local - with_offset)
```
结果显示如下：
```
numpy.timedelta64(-480,'m')
```
减法创建一个timedelta64 NumPy 对象，在这种情况下，它表示 480 分钟的增量。

工作原理

您了解了datetime64 NumPy 类型。这种数据类型使我们可以轻松地操纵日期和时间。它的功能包括简单的算术运算和使用常规 NumPy 函数创建数组。请参考本书代码包中的datetime_demo.py文件：

import numpy as np

print(np.datetime64('2015-05-21'))
#numpy.datetime64('2015-05-21')

print(np.datetime64('20150521'))
print(np.datetime64('2015-05'))
#numpy.datetime64('20150521')
#numpy.datetime64('2015-05')

local = np.datetime64('1578-01-01T21:18')
print(local)
#numpy.datetime64('1578-01-01T21:18Z')

with_offset = np.datetime64('1578-01-01T21:18-0800')
print(with_offset)
#numpy.datetime64('1578-01-02T05:18Z')

print(local - with_offset)

另见

十二、使用 NumPy 进行探索性和预测性数据分析

在本章中，我们涵盖以下秘籍：

探索气压
探索日常气压范围
研究年度气压平均值
分析最大可见度
用自回归模型预测气压
使用移动平均模型预测气压
研究年内平均气压
研究气压的极值

简介

数据分析是 NumPy 最重要的用例之一。根据我们的目标，我们可以区分数据分析的许多阶段和类型。在本章中，我们将讨论探索性和预测性数据分析。探索性数据分析可探查数据的线索。在此阶段，我们可能不熟悉数据集。预测分析试图使用模型来预测有关数据的某些信息。

数据来自荷兰气象局 KNMI。特别是 KNMI 总部位于 De Bilt 的气象站。在这些秘籍中，我们将检查气压和最大可见度。

我修改了文本数据并将其从 KNMI 转换为 NumPy 特定的.npy格式，并保存为40996 x 5数组。该数组包含五个变量的每日值：

YYYYMMDD格式的日期
每日平均气压
每日最高气压
每日最低气压
每日最大可见度

探索气压

在此秘籍中，我们将查看由 24 小时值计算出的每日平均海平面气压（0.1 hPa）。这包括打印描述性统计数据和可视化概率分布。在自然界中，我们经常处理正态分布，因此使用第 10 章，“使用 Scikits 进行正态性检验”会派上用场。

完整的代码在本书代码包的exploring.py文件中：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.adnorm import normal_ad

data = np.load('cbk12.npy')

# Multiply to get hPa values
meanp = .1 * data[:,1]

# Filter out 0 values
meanp = meanp[ meanp > 0]

# Get descriptive statistics
print("Max", meanp.max())
print("Min", meanp.min())
mean = meanp.mean()
print("Mean", mean)
print("Median", np.median(meanp))
std = meanp.std()
print("Std dev", std)

# Check for normality
print("Normality", normal_ad(meanp))

#histogram with Gaussian PDF
plt.subplot(211)
plt.title('Histogram of average atmospheric pressure')
_, bins, _ = plt.hist(meanp, np.sqrt(len(meanp)), normed=True)
plt.plot(bins, 1/(std * np.sqrt(2 * np.pi)) * np.exp(- (bins - mean)**2/(2 * std**2)), 'r-', label="Gaussian PDF")
plt.grid()
plt.legend(loc='best')
plt.xlabel('Average atmospheric pressure (hPa)')
plt.ylabel('Frequency')

# boxplot
plt.subplot(212)
plt.boxplot(meanp)
plt.title('Boxplot of average atmospheric pressure')
plt.ylabel('Average atmospheric pressure (hPa)')
plt.grid()

# Improves spacing of subplots
plt.tight_layout()
plt.show()

准备

如果尚未安装，请安装statsmodels以进行正态性测试（请参阅第 10 章，“Scikits 的乐趣”的“安装 scikits-statsmodels”秘籍）。

操作步骤

请按照以下步骤探索每日气压：

使用load()函数加载数据：
```
data = np.load('cbk12.npy')
```
通常，数据需要进行处理和清理。在这种情况下，将这些值相乘以获得hPa中的值，然后删除与缺失值相对应的0值：
```
# Multiply to get hPa values
meanp = .1 * data[:,1]

# Filter out 0 values
meanp = meanp[ meanp > 0]
```

获取描述性统计信息，包括最大值，最小值，算术平均值，中位数和标准差：

print("Max", meanp.max())
print("Min", meanp.min())
mean = meanp.mean()
print("Mean", mean)
print("Median", np.median(meanp))
std = meanp.std()
print("Std dev", std)

您应该看到以下值：

Max 1048.3
Min 962.1
Mean 1015.14058231
Median 1015.8
Std dev 9.85889134337

如下应用第 10 章“Scikits 的乐趣”的正态性检验：
```
print("Normality", normal_ad(meanp))
```
屏幕上显示以下值：
```
Normality (72.685781095773564, 0.0)
```
用直方图和箱形图可视化值的分布也很好。最终结果请参考以下图表：

另见

“使用 Statsmodels”秘籍执行正态性测试，该秘籍来自第 10 章，“Scikits 的乐趣”
有关箱形图的说明，请参见第 11 章“最新最强的 NumPy”
load()函数的文档

探索日常气压范围

每日气压范围是每日最高点和最低点之差。使用实际数据，有时我们会缺少一些值。在这里，我们可能会缺少给定一天的高压和/或低压值。可以通过智能算法来填补这些空白。但是，让我们保持简单，而忽略它们。在计算了范围之后，我们将进行与前面的秘籍类似的分析，但是我们将使用可以处理NaN值的函数。另外，我们将研究月份和范围之间的关系。

本书代码包中的day_range.py文件位于对应的代码中：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
import calendar as cal

data = np.load('cbk12.npy')

# Multiply to get hPa values
highs = .1 * data[:,2]
lows = .1 * data[:,3]

# Filter out 0 values
highs[highs == 0] = np.nan
lows[lows == 0] = np.nan

# Calculate range and stats
ranges = highs - lows
print("Minimum daily range", np.nanmin(ranges))
print("Maximum daily range", np.nanmax(ranges))

print("Average daily range", np.nanmean(ranges))
print("Standard deviation", np.nanstd(ranges))

# Get months
dates = data[:,0]
months = (dates % 10000)/100
months = months[~np.isnan(ranges)]

monthly = []
month_range = np.arange(1, 13)

for month in month_range:
   indices = np.where(month == months)
   monthly.append(np.nanmean(ranges[indices]))

plt.bar(month_range, monthly)
plt.title('Monthly average of daily pressure ranges')
plt.xticks(month_range, cal.month_abbr[1:13])
plt.ylabel('Monthly Average (hPa)')
plt.grid()
plt.show()

操作步骤

此秘籍的第一步与先前的秘籍几乎相同，因此我们将跳过它们。继续分析每日气压范围：

我们可以将缺失值保留为当前的0值。但是，通常将它们设置为NaN较为安全，以避免造成混淆。将缺失值设置为NaN，如下所示：
```
highs[highs == 0] = np.nan
lows[lows == 0] = np.nan
```

使用nanmin()，nanmax()，nanmean()和nanstd()函数计算范围，最小值，最大值，均值和标准差：

ranges = highs - lows
print("Minimum daily range", np.nanmin(ranges))
print("Maximum daily range", np.nanmax(ranges))

print("Average daily range", np.nanmean(ranges))
print("Standard deviation", np.nanstd(ranges))

结果出现在屏幕上：

Minimum daily range 0.4
Maximum daily range 41.7
Average daily range 6.11945360571
Standard deviation 4.42162136692

如前所述，日期以YYYYMMDD格式给出。有了一点算术，我们就可以轻松地得到几个月。此外，我们忽略与NaN范围值相对应的月份值：
```
dates = data[:,0]
months = (dates % 10000)/100
months = months[~np.isnan(ranges)]
```
按月对范围进行平均，如下所示：
```
monthly = []
month_range = np.arange(1, 13)

for month in month_range:
   indices = np.where(month == months)
   monthly.append(np.nanmean(ranges[indices]))
```
在最后一步中，我们绘制了matplotlib条形图，显示了每日气压范围的每月平均值。最终结果请参考以下图表：

工作原理

我们分析了气压的每日范围。此外，我们可视化了每日范围的每月平均值。似乎有一种模式导致夏季的每日气压范围变小。当然，需要做更多的工作来确定。

另见

“探索气压”秘籍

研究年度气压平均值

您可能听说过，全球变暖表明温度逐年稳定上升。由于气压是另一个热力学变量，因此我们可以预期气压也会跟随趋势。该秘籍的完整代码在本书代码包的annual.py文件中：

import numpy as np
import matplotlib.pyplot as plt

data = np.load('cbk12.npy')

# Multiply to get hPa values
avgs = .1 * data[:,1]
highs = .1 * data[:,2]
lows = .1 * data[:,3]

# Filter out 0 values
avgs = np.ma.array(avgs, mask = avgs == 0)
lows = np.ma.array(lows, mask = lows == 0)
highs = np.ma.array(highs, mask = highs == 0)

# Get years
years = data[:,0]/10000

# Initialize annual stats arrays
y_range = np.arange(1901, 2014)
nyears = len(y_range)
y_avgs = np.zeros(nyears)
y_highs = np.zeros(nyears)
y_lows = np.zeros(nyears)

# Compute stats
for year in y_range:
   indices = np.where(year == years)
   y_avgs[year - 1901] = np.mean(avgs[indices])
   y_highs[year - 1901] = np.max(highs[indices])
   y_lows[year - 1901] = np.min(lows[indices])

plt.title('Annual atmospheric pressure for De Bilt(NL)')
plt.ticklabel_format(useOffset=900, axis='y')

plt.plot(y_range, y_avgs, label='Averages')

# Plot ignoring NaNs
h_mask = np.isfinite(y_highs)
plt.plot(y_range[h_mask], y_highs[h_mask], '^', label='Highs')

l_mask = np.isfinite(y_lows)
plt.plot(y_range[l_mask], y_lows[l_mask], 'v', label='Lows')

plt.xlabel('Year')
plt.ylabel('Atmospheric pressure (hPa)')
plt.grid()
plt.legend(loc='best')
plt.show()

操作步骤

要检查趋势，让我们按以下步骤绘制平均，最大和最小年度气气压：

初始化年度统计数据数组：

y_range = np.arange(1901, 2014)
nyears = len(y_range)
y_avgs = np.zeros(nyears)
y_highs = np.zeros(nyears)
y_lows = np.zeros(nyears)

计算年度统计数据：

for year in y_range:
   indices = np.where(year == years)
   y_avgs[year - 1901] = np.mean(avgs[indices])
   y_highs[year - 1901] = np.max(highs[indices])
   y_lows[year - 1901] = np.min(lows[indices])

绘制，忽略NaN值，如下所示：

h_mask = np.isfinite(y_highs)
plt.plot(y_range[h_mask], y_highs[h_mask], '^', label='Highs')

l_mask = np.isfinite(y_lows)
plt.plot(y_range[l_mask], y_lows[l_mask], 'v', label='Lows')

最终结果请参考以下图表：

工作原理

年平均气压似乎持平或有所波动，但没有任何趋势。我们使用isfinite()函数忽略了最终图中的NaN值。此函数检查无限值和NaN值。

另见

“探索气压”秘籍
isfinite()函数的文档

分析最大可见度

如果到目前为止，您在本章的所有秘籍中均未使用，则可能需要突破气压。因此，让我们来看看能见度。数据文件具有一列，以提供最大的可视性，KNMI 描述如下：

Maximum visibility; 0: <100 m, 1:100-200 m, 2:200-300 m,..., 49:4900-5000 m, 50:5-6 km, 56:6-7 km, 57:7-8 km,..., 79:29-30 km, 80:30-35 km, 81:35-40 km,..., 89: >70 km)

这里的能见度是一个离散变量，因此取平均值可能没有意义。此外，似乎我们几乎每天都有 1901 年到 1950 年之间的0值。我不相信 De Bilt 在那个时期会更加迷雾重重。出于本秘籍的目的，我们将雾定义为 1 至 2km 之间的能见度，这对应于数据文件中的10和20值。让我们还将雾度定义为 2 至 5 公里之间的能见度。这又对应于我们数据文件中的20和50。

空气污染会降低能见度，尤其是在晴天。我们可以将晴天定义为能见度高于 30km 的晴天，或数据文件中79的值。理想情况下，我们应该使用空气污染数据，但不幸的是，我们没有这些数据。据我所知，这个特殊气象站周围的空气污染水平不是很高。知道每年的晴天数很有趣。分析的代码在本书的代码包的visibility.py文件中：

import numpy as np
import matplotlib.pyplot as plt

data = np.load('cbk12.npy')

# Get minimum visibility
visibility = data[:,4]

# doy
doy = data[:,0] % 10000

doy_range = np.unique(doy)

# Initialize arrays
ndoy = len(doy_range)
mist = np.zeros(ndoy)
haze = np.zeros(ndoy)

# Compute frequencies
for i, d in enumerate(doy_range):
   indices = np.where(d == doy)
   selection = visibility[indices]

   mist_truth = (10 < selection) & (selection < 20)
   mist[i] = len(selection[mist_truth])/(1\. * len(selection))

   haze_truth = (20 < selection) & (selection < 50)
   haze[i] = len(selection[haze_truth])/(1\. * len(selection))

# Get years
years = data[:,0]/10000

# Initialize annual stats arrays
y_range = np.arange(1901, 2014)
nyears = len(y_range)
y_counts = np.zeros(nyears)

# Get annual counts
for year in y_range:
   indices = np.where(year == years)
   selection = visibility[indices]
   y_counts[year - 1901] = len(selection[selection > 79])

plt.subplot(211)
plt.plot(np.arange(1, 367), mist, color='.25', label='mist')
plt.plot(np.arange(1, 367), haze, color='0.75', lw=2, label='haze')
plt.title('Probability of mist and haze')
plt.xlabel('Day of the year')
plt.ylabel('Probability')
plt.grid()
plt.legend(loc='best')

plt.subplot(212)
plt.plot(y_range, y_counts)
plt.xlabel('Year')
plt.ylabel('Number of clear days')
plt.title('Annual counts of clear days')
plt.grid()
plt.tight_layout()
plt.show()

操作步骤

请按照以下步骤来绘制晴天的年度计数，即一年中的某一天（1-366）对霾和薄雾的可能性：

使用以下代码块计算雾霾的可能性：

for i, d in enumerate(doy_range):
   indices = np.where(d == doy)
   selection = visibility[indices]

   mist_truth = (10 < selection) & (selection < 20)
   mist[i] = len(selection[mist_truth])/(1\. * len(selection))

   haze_truth = (20 < selection) & (selection < 50)
   haze[i] = len(selection[haze_truth])/(1\. * len(selection))

使用此代码段获取年度计数：

for year in y_range:
   indices = np.where(year == years)
   selection = visibility[indices]
   y_counts[year - 1901] = len(selection[selection > 79])

Refer to the following plot for the end result:

工作原理

如您所见，我们在 1950 年以后开始变得晴朗。这不是由于 1950 年之前的大雾天气，而是由于数据丢失或无效的现象。去年的下降也是由于数据不完整。 1980 年以后，晴天明显增加。这应该是全球变暖和气候变化也加剧的时期。不幸的是，我们没有与空气污染直接相关的数据，但是我们的探索性分析表明存在趋势。

薄雾似乎主要发生在一年的头两个月中。您可以得出有关雾霾的类似结论。显然，雾气比雾气更有可能，这可能是一件好事。您也可以绘制直方图以确保。但是，请记住，正如我之前提到的，您需要忽略0值。

另见

“探索气压”秘籍
“研究年度气压平均值”秘籍
相关维基百科页面

使用自回归模型预测气压

非常简单的预测模型会获取变量的当前值，并将其外推到下一个周期。为了推断，我们可以使用一个简单的数学函数。由于可以通过多项式近似泰勒级数中的多项式，因此低阶多项式可以解决问题。归结为将先前值回归到下一个值。因此，相应的模型称为自回归。

我们必须对此谨慎，以免过拟合。 交叉验证是一种常见的方法，将数据分为训练和测试集。我们使用训练集拟合数据，并使用测试集测试拟合度。这样可以减少偏差（请参见本书代码包中的autoregressive.py文件）：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

data = np.load('cbk12.npy')

# Load average pressure
meanp = .1 * data[:,1]

# Split point for test and train data
cutoff = 0.9 * len(meanp)

for degree, marker in zip(xrange(1, 4), ['o', 'x','.']):
   poly = np.polyfit(meanp[:cutoff - 1], meanp[1:cutoff], degree)
   print('Polynomial coefficients', poly)

   fit = np.polyval(poly, meanp[cutoff:-1])
   error = np.abs(meanp[cutoff + 1:] - fit)/fit
   plt.plot(error, marker, color=str(.25* degree), label='Degree ' + str(degree))
   plt.plot(np.full(len(error), error.mean()), lw=degree, label='Mean for degree ' + str(degree))

   print("Absolute mean relative error", error.mean(), 'for polynomial of degree', degree)
   print()

plt.title('Relative test errors for polynomial fits')
plt.ylabel('Relative error')
plt.grid()
plt.legend(loc='best')
plt.show()

操作步骤

通过执行以下，我们将使用不同程度的多项式拟合气压：

定义测试和训练集的截止：
```
cutoff = 0.9 * len(meanp)
```

使用polyfit()和polyval()函数拟合数据：

poly = np.polyfit(meanp[:cutoff - 1], meanp[1:cutoff], degree)
print('Polynomial coefficients', poly)

fit = np.polyval(poly, meanp[cutoff:-1])

计算相对误差：

error = np.abs(meanp[cutoff + 1:] - fit)/fit

此代码输出以下输出：

Polynomial coefficients [ 0.995542    4.50866543]
Absolute mean relative error 0.00442472512506 for polynomial of degree 1

Polynomial coefficients [ -1.79946321e-04   1.17995347e+00   2.77195814e+00]
Absolute mean relative error 0.00421276856088 for polynomial of degree 2

Polynomial coefficients [  3.17914507e-06  -6.62444552e-03   4.44558056e+00   2.76520065e+00]
Absolute mean relative error 0.0041906802632 for polynomial of degree 3

Refer to the following plot for the end result:

工作原理

这三个多项式的平均相对误差非常接近 -- 约 .004 -- 因此我们在图中看到一条线（很有趣的是知道气压的典型测量误差是多少），小于百分之一。我们看到了一些潜在的异常值，但是没有太多。大部分繁重的工作都是通过polyfit()和polyval()函数完成的，它们分别将数据拟合到多项式并求值该多项式。

另见

使用移动平均模型预测气压

模拟气压的简单方法是假设值围绕均值μ跳动。然后，在最简单的情况下，我们假设连续值与平均值的偏差ε遵循以下公式：

该关系是线性的，在最简单的情况下，我们只需要估计一个参数 - θ。为此，我们将需要 SciPy 功能。该秘籍的完整代码在本书代码包的moving_average.py文件中：

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime as dt
from scipy.optimize import leastsq

data = np.load('cbk12.npy')

# Load average pressure
meanp = .1 * data[:,1]

cutoff = 0.9 * len(meanp)

def model(p, ma1):
   return p * ma1

def error(p, t, ma1):
   return t - model(p, ma1)

p0 = [.9]
mu = meanp[:cutoff].mean()
params = leastsq(error, p0, args=(meanp[1:cutoff] - mu, meanp[:cutoff-1] - mu))[0]
print(params)

abs_error = np.abs(error(params, meanp[cutoff+1:] - mu, meanp[cutoff:-1] - mu))
plt.plot(abs_error, label='Absolute error')
plt.plot(np.full_like(abs_error, abs_error.mean()), lw=2, label='Absolute mean error')
plt.title('Absolute error for the moving average model')
plt.ylabel('Absolute error (hPa)')
plt.grid()
plt.legend(loc='best')
plt.show()

入门

如有必要，请按照第 2 章，“高级索引和数组概念”的“安装 SciPy”秘籍中的说明安装 SciPy。

操作步骤

以下步骤将移动平均模型应用于气压。

定义以下函数：

def model(p, ma1):
   return p * ma1

def error(p, t, ma1):
   return t - model(p, ma1)

使用上一步中的函数将移动平均模型与leastsq()函数拟合，并将0.9的初始猜测作为模型参数：

p0 = [.9]
mu = meanp[:cutoff].mean()
params = leastsq(error, p0, args=(meanp[1:cutoff] - mu, meanp[:cutoff-1] - mu))[0]

使用测试数据集拟合后计算绝对误差：
```
abs_error = np.abs(error(params, meanp[cutoff+1:] - mu, meanp[cutoff:-1] - mu))
```
请参考以下数据集中每个数据点的绝对误差图：

工作原理

leastsq()函数通过最小化误差来拟合模型。它需要一个函数来计算拟合误差和模型参数的初始猜测。

另见

研究年内平均气压

在一年内探索气压很有趣。特别地，检查与变异性相关的模式可能是有益的，因此，与可预测性相关。原因是几个月中的气气压会发生很大变化，从而降低了可预测性。在此秘籍中，我们将绘制每月的箱形图和每月的气气压方差。

秘籍代码在本书代码包的intrayear.py文件中。请特别注意突出显示的部分：

import numpy as np
import matplotlib.pyplot as plt
import calendar as cal

data = np.load('cbk12.npy')

# Multiply to get hPa values
meanp = .1 * data[:,1]

# Get months
dates = data[:,0]
months = (dates % 10000)/100

monthly = []
vars = np.zeros(12)
month_range = np.arange(1, 13)

for month in month_range:
 indices = np.where(month == months)
 selection = meanp[indices]

 # Filter out 0 values
 selection = selection[selection > 0]

 monthly.append(selection)
 vars[month - 1] = np.var(selection)

def plot():
    plt.xticks(month_range, cal.month_abbr[1:13])
    plt.grid()
    plt.xlabel('Month')

plt.subplot(211)
plot()
plt.title('Atmospheric pressure box plots')
plt.boxplot(monthly)
plt.ylabel('Atmospheric pressure (hPa)')

plt.subplot(212)
plot()

# Display error bars using standard deviation
plt.errorbar(month_range, vars, yerr=vars.std())
plt.plot(month_range, np.full_like(month_range, np.median(vars)), lw=3, label='Median')

# Shades the region above the median
plt.fill_between(month_range, vars, where=vars>np.median(vars), color='0.5')
plt.title('Variance of atmospheric pressure')
plt.ylabel('Variance')
plt.legend(loc='best')

plt.show()

操作步骤

在我们探索时，往往会重复这些步骤，并且此秘籍与本书中的其他秘籍之间存在重叠。以下是此秘籍中的新步骤：

使用标准差显示误差线：

plt.errorbar(month_range, vars, yerr=vars.std())

用高于中值的值对图的区域进行着色：
```
plt.fill_between(month_range, vars, where=vars>np.median(vars), color='0.5')
```
Refer to the following plot for the end result:

工作原理

我们将几个月与气气压的测量相匹配。我们使用匹配来绘制箱形图并可视化每月变化。这项研究表明，在 1 月，2 月，11 月和 12 月的最冷月份，气压变化高于中值。从图中可以看出，在温暖的夏季，气压范围狭窄。这与其他秘籍的结果一致。

另见

“探索气压”秘籍
“研究年度气压平均值”秘籍
var()的文档

研究气压的极值

异常值是一个问题，因为它们会影响我们对数据的理解。在此秘籍中，我们将异常值定义为与数据的第一或第三四分位数相差至少四分位数范围的 1.5 倍。四分位数间距是第一和第三四分位数之间的距离。让我们计算一年中每个月的异常值。完整的代码在本书代码包的extreme.py文件中：

import numpy as np
import matplotlib.pyplot as plt
import calendar as cal

data = np.load('cbk12.npy')

# Multiply to get hPa values
meanp = .1 * data[:,1]

# Filter out 0 values
meanp = np.ma.array(meanp, mask = meanp == 0)

# Calculate quartiles and irq
q1 = np.percentile(meanp, 25)
median = np.percentile(meanp, 50)
q3 = np.percentile(meanp, 75)

irq = q3 - q1

# Get months
dates = data[:,0]
months = (dates % 10000)/100

m_low = np.zeros(12)
m_high = np.zeros(12)
month_range = np.arange(1, 13)

for month in month_range:
   indices = np.where(month == months)
   selection = meanp[indices]
   m_low[month - 1] = len(selection[selection < (q1 - 1.5 * irq)])
   m_high[month - 1] = len(selection[selection > (q3 + 1.5 * irq)])

plt.xticks(month_range, cal.month_abbr[1:13])
plt.bar(month_range, m_low, label='Low outliers', color='.25')
plt.bar(month_range, m_high, label='High outliers', color='0.5')
plt.title('Atmospheric pressure outliers')
plt.xlabel('Month')
plt.ylabel('# of outliers')
plt.grid()
plt.legend(loc='best')
plt.show()

操作步骤

要绘制一年中每个月的异常数，请执行以下步骤：

使用percentile()函数计算四分位数和四分位数间距：

q1 = np.percentile(meanp, 25)
median = np.percentile(meanp, 50)
q3 = np.percentile(meanp, 75)

irq = q3 - q1

计算异常值的数量，如下所示：

for month in month_range:
   indices = np.where(month == months)
   selection = meanp[indices]
   m_low[month - 1] = len(selection[selection < (q1 - 1.5 * irq)])
   m_high[month - 1] = len(selection[selection > (q3 + 1.5 * irq)])

Refer to the following plot for the end result:

工作原理

看起来异常值大多在下侧，并且夏天不太可能出现。较高端的异常值似乎仅在某些月份内发生。我们发现四分位数具有percentile()函数，因为四分之一对应 25%。

另见

“探索气压”秘籍
percentile()函数的文档

posted @ 2025-10-23 15:10 绝不原创的飞龙阅读(18) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈

NumPy-秘籍中文第二版-全-

NumPy 秘籍中文第二版（全）

零、前言

本书涵盖的内容

这本书需要什么

这本书适合谁？

约定

注意

提示

一、使用 IPython

简介

安装 IPython

操作步骤

工作原理

另见

使用 IPython 作为 Shell

操作步骤

工作原理

另见

阅读手册页

操作步骤

工作原理

安装 matplotlib

操作步骤

另见

运行 IPython 笔记本

准备

操作步骤

工作原理

另见

导出 IPython 笔记本

操作步骤

注意

导入网络笔记本

操作步骤

配置笔记本服务器

操作步骤

工作原理

另见

探索 SymPy 配置文件

准备

操作步骤

另见

二、高级索引和数组概念

简介

注意

安装 SciPy

准备

操作步骤

工作原理

另见

安装 PIL

操作步骤

另见

调整图像大小

注意

准备

操作步骤

工作原理

另见

创建视图和副本

准备

操作步骤

工作原理

另见

翻转 Lena

操作步骤

另见

花式索引

操作步骤

工作原理

另见

位置列表索引

操作步骤

另见

布尔值索引

操作步骤

工作原理