Python-数据分析教程-二-

Python 数据分析教程（二）

原文：Python Data Analysis

协议：CC BY-NC-SA 4.0

四、统计学和线性代数

统计学和线性代数为探索性数据分析奠定了基础。两种主要的统计方法，描述性的和推断性的，都有助于从原始数据中获得见解和做出推断。例如，我们可以计算出一个变量的数据有一定的算术平均值和标准偏差。从这些数字中，我们可以推断出这个变量的范围和期望值。然后，我们可以运行统计测试来检查我们得出正确结论的可能性有多大。

线性代数关注的是线性方程组。使用linalg包，用 NumPy 和 SciPy 很容易解决这些问题。例如，线性代数对于将数据拟合到模型中是有用的。我们将在本章中介绍用于随机数生成和屏蔽数组的其他 NumPy 和 SciPy 包。

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

• NumPy 的基本描述性统计
• 带有 NumPy 的 Linera 代数
• 用 NumPy 求特征值和特征向量
• 随机数
• 创建数字掩码数组

用 NumPy 进行基本描述性统计

在本书中，我们将尝试使用尽可能多的不同数据集。这取决于数据的可用性。不幸的是，这意味着数据的主题可能与你的兴趣不完全匹配。每个数据集都有自己的怪癖，但是你在本书中获得的一般技能应该转移到你自己的领域。在本章中，我们将把statsmodels库中的数据集加载到 NumPy 数组中，以便对数据进行分析。

莫纳罗亚公司 ₂ 测量是我们将从statsmodels数据集包中使用的第一个数据集。以下代码加载数据集并打印基本的描述性统计值:

import numpy as np
from scipy.stats import scoreatpercentile
import pandas as pd

data = pd.read_csv("co2.csv", index_col=0, parse_dates=True) 
co2 = np.array(data.co2) 

print("The statistical values for amounts of co2 in atmosphere: \n") 
print("Max method", co2.max()) 
print("Max function", np.max(co2)) 

print("Min method", co2.min()) 
print("Min function", np.min(co2)) 

print("Mean method", co2.mean()) 
print("Mean function", np.mean(co2)) 

print("Std method", co2.std()) 
print("Std function", np.std(co2)) 

print("Median", np.median(co2)) 
print("Score at percentile 50", scoreatpercentile(co2, 50))

前面的代码计算 NumPy 数组的平均值、中值、最大值、最小值和标准偏差。

注

如果这些术语听起来不熟悉，请花些时间从维基百科或任何其他来源了解它们。如前言所述，我们将假设您熟悉基本的数学和统计概念。

数据来自statsmodels包，包含美国夏威夷莫纳罗亚天文台的大气co2值

现在，让我们浏览一下代码:

1. 首先，加载 Python 包模块的常用导入语句如下:

           import numpy as np 
           from scipy.stats import scoreatpercentile 
           import pandas as pd 
   
   

2. Next, we will load the data with the following lines of code:

           data = pd.read_csv("co2.csv", index_col=0, parse_dates=True) 
           co2 = np.array(data.co2) 
   
   

   前面代码中的数据被复制到 NumPy 数组中，co2，包含大气中co2的值。

3. The maximum of an array can be obtained via a method of the ndarray and NumPy functions. The same goes for the minimum, mean, and standard deviations. The following code snippet prints the various statistics:

           print("Max method", co2.max())
           print("Max function", np.max(co2))
   
           print("Min method", co2.min())
           print("Min function", np.min(co2))
   
           print("Mean method", co2.mean())
           print("Mean function", np.mean(co2))
   
           print("Std method", co2.std())
           print("Std function", np.std(co2))
   
   

   输出如下:

              Max method 366.84 
              Max function 366.84 
              Min method 313.18 
              Min function 313.18 
              Mean method 337.053525641 
              Mean function 337.053525641
              Std method 14.9502216262
              Std function 14.9502216262
   

4. The median can be retrieved with a NumPy or SciPy function, which can estimate the 50th percentile of the data with the following lines:

           print("Median", np.median(co2))
           print("Score at percentile 50",scoreatpercentile(co2, 50)) 
   
   

   将打印以下内容:

              Median 335.17 
              Score at percentile 50 335.17
   

带 NumPy 的线性代数

线性代数是数学的一个重要分支。例如，我们可以使用线性代数来执行线性回归。numpy.linalg子包保存线性代数例程。有了这个子包，您可以反转矩阵、计算特征值、求解线性方程、寻找行列式等等。NumPy 中的矩阵由ndarray的子类表示。

用 NumPy 求矩阵的逆

线性代数中平方可逆矩阵A的逆是矩阵A<sup>-1</sup>，当与原矩阵相乘时，等于单位矩阵I。这可以写成下面的数学方程:

A A-1 = I 

numpy.linalg子包中的inv()功能可以为我们做到这一点。让我们反转一个示例矩阵。要反转矩阵，请执行以下步骤:

1. Create the demonstration matrix with the mat() function:

           A = np.mat("2 4 6;4 2 6;10 -4 18")
            print("A\n", A)
   

   A矩阵打印如下:

             A[[ 2  4  6] [ 4  2  6] [10 -4 18]]]
   

2. Invert the matrix.

   现在我们可以查看正在运行的inv()子程序:

             inverse = np.linalg.inv(A)
              print("inverse of A\n", inverse)
   

   逆矩阵显示如下:

             inverse of A[
             [-0.41666667 0.66666667 -0.08333333] 
             [ 0.08333333 0.16666667 -0.08333333] 
             [ 0.25 -0.33333333 0.08333333]]
   

   注

   如果矩阵是奇异的，或者不是正方形的，则会出现一个LinAlgError。如果您愿意，您可以手动检查解决方案。这是留给你的一个练习，在这个练习之外。pinv() NumPy 函数执行伪求逆，可以应用于任何矩阵，包括非正方形的矩阵。

3. Check by multiplication.

   让我们检查一下当我们将原始矩阵乘以inv()函数的结果时，我们得到了什么:

             print("Check\n", A * inverse) 
   
   

   结果是身份矩阵，正如预期的那样(忽略小的差异):

             Check
             [[  1.00000000e+00   0.00000000e+00  -5.55111512e-17] 
             [ -2.22044605e-16   1.00000000e+00  -5.55111512e-17] 
             [ -8.88178420e-16   8.88178420e-16   1.00000000e+00]]
   

   通过从前面的结果中减去 3×3 单位矩阵，我们得到反演过程的误差:

   print("Error\n", A * inverse - np.eye(3)) 
   
   

   一般来说，误差应该可以忽略不计，但在某些情况下，小误差可能会传播并带来不良副作用:

   [[ -1.11022302e-16   0.00000000e+00  -5.55111512e-17] 
    [ -2.22044605e-16   4.44089210e-16  -5.55111512e-17] 
    [ -8.88178420e-16   8.88178420e-16  -1.11022302e-16]]
   

   在这种情况下，更高精度的数据类型可能会有所帮助或切换到更好的算法。我们使用numpy.linalg子包的inv()例程计算矩阵的逆矩阵。我们用矩阵乘法确定，这是否真的是逆矩阵:

   import numpy as np 
   
   A = np.mat("2 4 6;4 2 6;10 -4 18") 
   print "A\n", A 
   
   inverse = np.linalg.inv(A) 
   print("inverse of A\n", inverse)
   
   print("Check\n", A * inverse)
   print("Error\n", A * inverse - np.eye(3))
   
   

用 NumPy 求解线性系统

矩阵以线性方式将一个向量转换成另一个向量。这个操作在数值上对应于一个线性方程组。numpy.linalg的solve()子程序求解形式为Ax = b的线性方程组；这里，A是矩阵，b可以是一维或二维数组，x是未知量。我们将见证dot()子程序的运行。这个函数计算两个浮点数数组的点积。

让我们解决一个线性系统的例子。要求解线性系统，请执行以下步骤:

1. Create matrix A and array b.

   下面的代码将创建矩阵A和数组b:

             A = np.mat("1 -2 1;0 2 -8;-4 5 9")
              print("A\n", A)
              b = np.array([0, 8, -9])
              print("b\n", b)
   

   矩阵A和数组(向量)b定义如下:

   

2. 调用solve()功能。

3. Solve this linear system with the solve() function:

             x = np.linalg.solve(A, b)
              print("Solution", x)
   

   线性系统的解如下:

             Solution [ 29\.  16\.   3.]
   

4. Check with the dot() function.

   使用dot()功能检查解决方案是否正确:

             print("Check\n", np.dot(A , x)) 
   
   

   结果是意料之中的:

             Check[[ 0\.  8\. -9.]]
   

我们通过应用 NumPy 的linalg子包中的solve()函数并使用dot()函数检查结果来求解线性系统:

import numpy as np 

A = np.mat("1 -2 1;0 2 -8;-4 5 9") 
print("A\n", A)

b = np.array([0, 8, -9]) 
print("b\n", b) 

x = np.linalg.solve(A, b) 
print("Solution", x)

print("Check\n", np.dot(A , x)) 

用 NumPy 求特征值和特征向量

eigenvalues是方程Ax = ax的标量解，其中A是二维矩阵，x是一维向量。eigenvectors是对应于eigenvalues的向量。

注

eigenvalues和eigenvectors是数学中的基础，用于许多重要的算法，如主成分分析 ( 主成分分析)。主成分分析可用于简化大型数据集的分析。

numpy.linalg包中的eigvals()子程序计算eigenvalues。eig()函数返回一个包含eigenvalues和eigenvectors的元组。

我们将获得具有numpy.linalg子包的eigvals()和eig()功能的矩阵的eigenvalues和eigenvectors。我们将通过应用dot()功能来检查结果:

import numpy as np 

A = np.mat("3 -2;1 0") 
print("A\n", A) 

print("Eigenvalues", np.linalg.eigvals(A)) 

eigenvalues, eigenvectors = np.linalg.eig(A) 
print("First tuple of eig", eigenvalues) 
print("Second tuple of eig\n", eigenvectors) 

for i in range(len(eigenvalues)): 
   print("Left", np.dot(A, eigenvectors[:,i])) 
   print("Right", eigenvalues[i] * eigenvectors[:,i]) 

让我们计算一个矩阵的eigenvalues:

1. Create the matrix.

   以下代码将创建一个矩阵:

             A = np.mat("3 -2;1 0")
              print("A\n", A)
   

   我们创建的矩阵如下所示:

             A[[ 3 -2] [ 1  0]]
   

2. Calculate eigenvalues with the eig() function.

   应用eig()子程序:

             print("Eigenvalues", np.linalg.eigvals(A)) 
   
   

   矩阵的eigenvalues如下:

             Eigenvalues [ 2\.  1.]
   

3. Get eigenvalues and eigenvectors with eig().

   使用eig()功能获取eigenvalues和eigenvectors。该例程返回一个元组，其中第一个元素包含eigenvalues，第二个元素包含匹配的eigenvectors，按列设置:

             eigenvalues, eigenvectors = np.linalg.eig(A)
              print("First tuple of eig", eigenvalues)
              print("Second tuple of eig\n", eigenvectors)
   

   eigenvalues和eigenvectors值如下:

             First tuple of eig [ 2\. 1.]
             Second tuple of eig
             [[ 0.89442719 0.70710678] 
             [ 0.4472136 0.70710678]]
   

4. Check the result.

   通过计算eigenvalues方程Ax = ax的两边，用dot()函数检查答案:

             for i in range(len(eigenvalues)):
              print("Left", np.dot(A, eigenvectors[:,i]))
              print("Right", eigenvalues[i] * eigenvectors[:,i])
   

   输出如下:

             Left [[ 1.78885438] [ 0.89442719]]
             Right [[ 1.78885438] [ 0.89442719]]
             Left [[ 0.70710678] [ 0.70710678]]
             Right [[ 0.70710678] [ 0.70710678]]
   

NumPy 随机数

随机数用于蒙特卡罗方法、随机微积分等。实随机数很难产生，所以在实践中，我们使用伪随机数。伪随机数对于大多数意图和目的来说都是足够随机的，除了一些非常特殊的情况，例如非常精确的模拟。随机数关联例程可以位于 NumPy random子包中。

注

核心随机数发生器基于默森扭转算法(参考en.wikipedia.org/wiki/Mersen…)。

随机数可以由离散或连续分布产生。分布函数有一个可选的size参数，它通知 NumPy 要创建多少个数字。您可以指定整数或元组作为大小。这将导致用随机数填充适当形状的数组。离散分布包括几何分布、超几何分布和二项式分布。连续分布包括正态分布和对数正态分布。

二项分布的赌博

二项式分布对一个实验的整数次独立运行的成功次数进行建模，其中每个实验的成功几率是一个设定值。

设想一个 17 世纪的赌场，在那里你可以赌掷 8 块钱。在一个受欢迎的游戏中，九枚硬币被抛了出来。如果少于五个硬币是头像，那么你就输了一块八；否则，你赚一个。让我们模拟一下，从我们拥有的一千枚硬币开始。为此，我们将使用random模块中的binomial()功能:

注

如果你想跟着代码走，看看本书代码包里的ch-04.ipynb。

import numpy as np 
from matplotlib.pyplot import plot, show 

cash = np.zeros(10000) 
cash[0] = 1000 
outcome = np.random.binomial(9, 0.5, size=len(cash)) 

for i in range(1, len(cash)): 

   if outcome[i] < 5: 
      cash[i] = cash[i - 1] - 1 
   elif outcome[i] < 10: 
      cash[i] = cash[i - 1] + 1 
   else: 
      raise AssertionError("Unexpected outcome " + outcome) 

print(outcome.min(), outcome.max()) 

plot(np.arange(len(cash)), cash) 
show() 

为了理解binomial()功能，请看以下步骤:

1. Calling the binomial() function.

   将充当现金余额的数组初始化为零。调用10000大小的binomial()函数。这代表我们赌场的 10，000 次抛硬币:

             cash = np.zeros(10000) 
             cash[0] = 1000 
             outcome = np.random.binomial(9, 0.5, size=len(cash)) 
   
   

2. Updating the cash balance.

   查看投币结果，更新cash数组。显示outcome数组的最高和最低值，只是为了确保我们没有任何异常的异常值:

             for i in range(1, len(cash)):
              if outcome[i] < 5:
              cash[i] = cash[i - 1] - 1
              elif outcome[i] < 10:
              cash[i] = cash[i - 1] + 1
              else:
              raise AssertionError("Unexpected outcome " + outcome)
              print(outcome.min(), outcome.max())
   

   不出所料，数值介于0和9之间:

              0 9
   

3. 接下来，用 matplotlib 绘制现金数组:

           plot(np.arange(len(cash)), cash) 
           show() 
   
   

在下面的图中，您可以确定我们的现金余额执行随机游走(不遵循模式的随机移动):

当然，每次执行代码，我们都会有不同的随机游走。如果您希望总是收到相同的结果，您将希望从 NumPy random子包中向binomial()函数传递一个种子值。

正态分布采样

连续分布由概率密度函数建模。指定间隔的机会是通过对 PDF 的整合找到的。NumPy random模块有许多代表连续分布的功能，如beta、chisquare、exponential、f、gamma、gumbel、laplace、lognormal、logistic、multivariate_normal、noncentral_chisquare、noncentral_f、normal等。

我们将通过应用random NumPy 子包中的normal()函数来可视化正态分布。我们将通过绘制钟形曲线和随机生成值的直方图来实现这一点:

import numpy as np 
import matplotlib.pyplot as plt 

N=10000 

normal_values = np.random.normal(size=N) 
dummy, bins, dummy = plt.hist(normal_values, np.sqrt(N), normed=True, lw=1) 
sigma = 1 
mu = 0 
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu)**2 / (2 * sigma**2) ),lw=2) 
plt.show() 

随机数可以由正态分布产生，它们的分布可以用直方图显示。要绘制正态分布，请执行以下步骤:

1. Generate values.

   借助random NumPy 子包中的normal()功能，为特定样本量创建随机数:

             N=100.00 
             normal_values = np.random.normal(size=N) 
   
   

2. Draw the histogram and theoretical PDF.

   绘制直方图和理论 PDF，中心值为0，标准偏差为1。我们将为此使用 matplotlib:

             dummy, bins, dummy = plt.hist(normal_values, np.sqrt(N),         
             normed=True, lw=1)
             sigma = 1
             mu = 0
             plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - 
             (bins - mu)**2 / (2 * sigma**2) ),lw=2)
             plt.show()
   

   在下面的图中，我们看到了著名的钟形曲线:

   

用 SciPy 进行正态性测试

正态分布广泛应用于科学和统计中。根据中心极限定理，具有独立观测值的大的随机样本将向正态分布收敛。正态分布的性质是众所周知的，并且被认为使用方便。然而，需要满足许多要求，例如必须独立的足够多的数据点。检查数据是否符合正态分布是很好的做法。存在大量的正态性测试，其中一些已经在scipy.stats包中实现。我们将在本节中应用这些测试。作为样本数据，我们将使用来自www.google.org/flutrends/d…的流感趋势数据。原始文件已缩减为仅包含两列-阿根廷的日期和值。几行如下:

Date,Argentina
29/12/02,
05/01/03,
12/01/03,
19/01/03,
26/01/03,
02/02/03,136

数据可以在代码包的goog_flutrends.csv文件中找到。我们还将从正态分布中采样数据，就像我们在前面的教程中所做的那样。由此产生的数组将具有与流感趋势数组相同的大小，并将作为黄金标准，它应该以优异的成绩通过常态测试。以下代码参见代码包中的ch-04.ipynb:

import numpy as np 
from scipy.stats import shapiro 
from scipy.stats import anderson 
from scipy.stats import normaltest 

flutrends = np.loadtxt("goog_flutrends.csv", delimiter=',', usecols=(1,), skiprows=1, converters = {1: lambda s: float(s or 0)}, unpack=True) 
N = len(flutrends) 
normal_values = np.random.normal(size=N) 

print("Normal Values Shapiro", shapiro(normal_values)) 
print("Flu Shapiro", shapiro(flutrends)) 

print("Normal Values Anderson", anderson(normal_values)) 
print("Flu Anderson", anderson(flutrends)) 

print("Normal Values normaltest", normaltest(normal_values)) 
print("Flu normaltest", normaltest(flutrends)) 

作为一个反面例子，我们将使用一个与前面提到的两个用零填充的数组大小相同的数组。在现实生活中，如果我们正在处理一个罕见的事件(例如，疫情疫情)，我们可以得到这些类型的价值。

在数据文件中，有些单元格是空的。当然，这类问题经常发生，因此我们必须习惯于清理数据。我们将假设正确的值应该是0。转换器可以为我们填写这些零值，如下所示:

flutrends = np.loadtxt("goog_flutrends.csv", delimiter=',', usecols=(1,), skiprows=1, converters = {1: lambda s: float(s or 0)}, unpack=True) 

夏皮罗-维尔克检验可以检验是否正常。对应的 SciPy 函数返回一个元组，其中第一个数字是测试统计量，第二个数字是 p 值。应该注意的是，用零填充的数组会引起警告。事实上，这个例子中使用的三个函数都有这个数组的问题，并给出了警告。我们得到以下结果:

Normal Values Shapiro (0.9967482686042786, 0.2774980068206787)
Flu Shapiro (0.9351990818977356, 2.2945883254311397e-15)

我们得到的 p 值类似于本例后面第三次测试的结果。分析基本一致。

安德森-达林检验可以检验正态性，也可以检验其他分布，如指数分布、逻辑分布和甘贝尔分布。SciPy 函数将测试统计与包含 15、10、5、2.5 和 1%显著性水平的临界值的数组相关联。如果统计量在显著性水平上大于临界值，我们可以拒绝正态性。我们得到以下值:

Normal Values Anderson (0.31201465602225653, array([ 0.572,  0.652,  0.782,  0.912,  1.085]), array([ 15\. ,  10\. ,   5\. ,   2.5,   1\. ]))
Flu Anderson (8.258614154768793, array([ 0.572,  0.652,  0.782,  0.912,  1.085]), array([ 15\. ,  10\. ,   5\. ,   2.5,   1\. ]))

正如我们所预料的那样，我们不能拒绝我们的黄金标准数组的常态。然而，为流感趋势数据返回的统计数据大于所有相应的临界值。因此，我们可以自信地拒绝常态。在三个测试函数中，这个似乎是最容易使用的。

达戈斯蒂诺-皮尔森测试也作为normaltest()函数在 SciPy 中实现。这个函数返回一个带有统计和 p 值的元组，就像shapiro()函数一样。p 值是双边卡方概率。卡方是另一个众所周知的分布。测试本身基于偏斜度和峰度测试的 z 分数。偏斜度衡量分布的对称性。正态分布是对称的，偏斜度为零。峰度告诉我们一些关于分布的形状(高峰，肥尾)。正态分布的峰度为 3(过度峰度为零)。通过测试获得以下值:

Normal Values normaltest (3.102791866779639, 0.21195189649335339)
Flu normaltest (99.643733363569538, 2.3048264115368721e-22)

因为我们处理的是 p 值的概率，所以我们希望这个概率尽可能高，并且接近 1。此外，如果 p 值至少为 0.5，我们可以接受正态性。对于黄金标准数组，我们得到一个较低的值，这意味着我们可能需要有更多的观察。留给你来证实这一点。

创建 NumPy 屏蔽数组

数据通常很混乱，包含我们不常处理的空白或字符。屏蔽数组可用于忽略缺失或无效的数据点。来自numpy.ma子包的屏蔽数组是带有屏蔽的ndarray的子类。在本节中，我们将使用face照片作为数据源，并表现得好像其中一些数据已损坏。以下是本书代码包中ch-04.ipynb文件中屏蔽数组示例的完整代码:

import numpy 
import scipy 
import matplotlib.pyplot as plt 

face = scipy.misc.face() 

random_mask = numpy.random.randint(0, 2, size=face.shape) 

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

masked_array = numpy.ma.array(face, mask=random_mask) 

plt.subplot(222) 
plt.title("Masked") 
plt.imshow(masked_array) 
plt.axis('off') 

plt.subplot(223) 
plt.title("Log") 
plt.imshow(numpy.ma.log(face).astype("float32")) 
plt.axis('off') 

plt.subplot(224) 
plt.title("Log Masked") 
plt.imshow(numpy.ma.log(masked_array).astype("float32")) 
plt.axis('off') 

plt.show() 

最后，我们将显示原始图片、原始图像的对数值、屏蔽数组及其对数值:

1. Create a mask.

   为了产生屏蔽数组，我们必须规定一个屏蔽。让我们创建一个随机掩码。该掩码的值为0或1:

             random_mask = numpy.random.randint(0, 2, size=face.shape) 
   
   

2. Create a masked array.

   通过应用前一步中的掩码，创建一个掩码数组:

             masked_array = numpy.ma.array(face, mask=random_mask) 
   
   

   结果图片展示如下:

   

我们对 NumPy 数组应用了随机掩码。这导致忽略与掩码匹配的数据。在numpy.ma子包中可以发现整个范围的屏蔽数组过程。在本教程中，我们只介绍了如何生成屏蔽数组。

忽略负值和极值

当我们想要忽略负值时，屏蔽数组很有用；例如，当取数组值的对数时。屏蔽数组的第二个用例是拒绝异常值。这基于极值的上限和下限。在本教程中，我们将把这些技术应用到 MLB 球员的工资数据中。数据来源于www.exploredata.net/Downloads/B…。数据被编辑为包含两列-球员姓名和工资。这导致了MLB2008.csv，可以在代码包中找到。本教程的完整脚本在本书代码包的ch-04.ipynb文件中:

import numpy as np 
from datetime import date 
import sys 
import matplotlib.pyplot as plt 

salary = np.loadtxt("MLB2008.csv", delimiter=',', usecols=(1,), skiprows=1, unpack=True) 
triples = np.arange(0, len(salary), 3) 
print("Triples", triples[:10], "...") 

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

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

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

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

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

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

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

plt.subplots_adjust(hspace=.9) 

plt.show() 

以下是将帮助您执行上述命令的步骤:

1. Taking the logarithm of negative numbers.

   我们将取包含负数的数组的对数。首先，让我们创建一个可以被 3 整除的数组:

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

   接下来，我们将生成一个数组，其大小与薪资数据数组相同:

             signs = numpy.ones(len(salary)) 
             print("Signs", signs[:10], "...") 
   
   

   借助我们在第 2 章、 NumPy 数组中获得的索引技巧，我们将每个第三数组元素设置为负:

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

   总之，我们将取这个数组的对数:

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

   这应该打印以下工资数据(注意-代表数据中的 NaN 值):

             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 [-- 14.970818190308929 15.830413578506539 --           
             13.458835614025542 15.319587954740548 -- 15.648092021712584 
             13.864300722133706 --] ...
   

2. Ignoring extreme values.

   让我们将异常值指定为“比平均值低一个标准差”或“比平均值高一个标准差”(这不一定是这里给出的正确定义，但很容易计算)。该定义指导我们编写以下代码，这些代码将掩盖极端点:

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

   以下代码显示了前十个元素的输出:

             Inside [3750000.0 3175000.0 7500000.0 3000000.0 700000.0 
             4500000.0 3000000.0 6250000.0 1050000.0 4600000.0] ...
   

   让我们绘制原始工资数据，再次取对数和指数后的数据，最后是应用基于标准差的掩码后的数据。

它看起来像这样:

numpy.ma子包中的函数屏蔽数组元素，我们认为这是无效的。例如，log()和sqrt()功能不允许负值。屏蔽值就像关系数据库和编程中的空值。具有屏蔽值的所有操作都会传递一个屏蔽值。

总结

在本章中，您学习了很多关于 NumPy 和 SciPy 子包的知识。我们研究了线性代数、统计学、连续和离散分布、屏蔽数组和随机数。

下一章第五章、检索、处理和存储数据，将教会我们一些至关重要的技能，尽管有些人可能不认为它们是数据分析。我们将使用一个更广泛的定义，考虑任何可以想象的与数据分析相关的东西。通常，当我们分析数据时，我们没有一个完整的助手团队来帮助我们检索和存储数据。然而，由于这些任务对于平稳的数据分析流程很重要，我们将详细描述这些活动。

五、检索、处理和存储数据

数据随处可见，有各种形状和形式。我们可以通过网络、电子邮件和文件传输协议获得它，或者我们可以自己在实验室实验或市场调查中创建它。全面概述如何获取各种格式的数据将需要比我们现有的更多的页面。有时，我们需要在分析数据之前或完成分析之后存储数据。我们将在本章中讨论存储数据。第 8 章、使用数据库，给出了关于各种数据库(关系数据库和 NoSQL 数据库)和相关应用编程接口的信息。以下是我们将在本章中讨论的主题列表:

• 用 NumPy 和 Pandas 编写 CSV 文件
• 二进制.npy和泡菜格式
• 用 PyTables 存储数据
• 向 HDF5 商店读写 Pandas 数据帧
• 和 Pandas 一起读写 Excel
• 使用 REST 网络服务和 JSON
• 和 Pandas 一起读写 JSON
• 解析 RSS 和 Atom 提要
• 用美丽的汤解析 HTML

用 NumPy 和 Pandas 编写 CSV 文件

在前面的章节中，我们学习了阅读 CSV 文件。编写 CSV 文件同样简单，但是使用不同的功能和方法。让我们首先生成一些以 CSV 格式存储的数据。在下面的代码片段中植入随机生成器后，生成一个 3x4 NumPy 数组。

将其中一个数组值设置为nan:

np.random.seed(42) 

a = np.random.randn(3, 4) 
a[2][2] = np.nan 
print(a) 

该代码将按如下方式打印数组:

[[ 0.49671415 -0.1382643   0.64768854  1.52302986]
 [-0.23415337 -0.23413696  1.57921282  0.76743473]
 [-0.46947439  0.54256004         nan -0.46572975]]

NumPy savetxt()函数是 NumPy loadtxt()函数的对应物，可以以分隔文件格式保存数组，例如 CSV。使用以下函数调用保存我们创建的数组:

np.savetxt('np.csv', a, fmt='%.2f', delimiter=',', header=" #1,  #2,  #3,  #4") 

在前面的函数调用中，我们指定了要保存的文件名、数组、可选格式、分隔符(默认为空格)和可选标题。

注

格式参数记录在http://docs . python . org/2/library/string . html # format-specification-mini-language。

查看我们使用cat命令(cat np.csv)或编辑器(如 Windows 中的记事本)创建的np.csv文件。文件的内容应显示如下:

    #  #1,  #2,  #3,  #4
    0.50,-0.14,0.65,1.52
    -0.23,-0.23,1.58,0.77
    -0.47,0.54,nan,-0.47

从随机值数组创建 Pandas 数据帧:

df = pd.DataFrame(a) 
print(df) 

如您所见，Pandas 会自动为我们的数据提供列名:

              0         1         2         3
    0  0.496714 -0.138264  0.647689  1.523030
    1 -0.234153 -0.234137  1.579213  0.767435
    2 -0.469474  0.542560NaN -0.465730

使用 Pandasto_csv()方法将数据帧写入 CSV 文件，如下所示:

df.to_csv('pd.csv', float_format='%.2f', na_rep="NAN!") 

我们给了这个方法文件名，一个类似于 NumPy savetxt()函数的格式参数的可选格式字符串，以及一个表示NaN的可选字符串。查看pd.csv文件可以看到以下内容:

    ,0,1,2,3
    0,0.50,-0.14,0.65,1.52
    1,-0.23,-0.23,1.58,0.77
    2,-0.47,0.54,NAN!,-0.47

看看本书代码包中ch-05.ipynb文件中的代码:

import numpy as np 
import pandas as pd 

np.random.seed(42) 

a = np.random.randn(3, 4) 
a[2][2] = np.nan 
print(a) 
np.savetxt('np.csv', a, fmt='%.2f', delimiter=',', header=" #1,  #2,  #3,  #4") 
df = pd.DataFrame(a) 
print(df) 
df.to_csv('pd.csv', float_format='%.2f', na_rep="NAN!") 

二进制。npy 和泡菜格式

大多数情况下，以 CSV 格式保存数据是可以的。交换 CSV 文件很容易，因为大多数编程语言和应用程序都可以处理这种格式。但是，效率不是很高；CSV 和其他明文格式占用大量空间。已经发明了许多提供高级压缩的文件格式，例如.zip、.bzip和.gzip。

以下是本次存储对比练习的完整代码，也可以在本书的代码包ch-05.ipynb文件中找到:

import numpy as np 
import pandas as pd 
from tempfile import NamedTemporaryFile 
from os.path import getsize 

np.random.seed(42) 
a = np.random.randn(365, 4) 

tmpf = NamedTemporaryFile() 
np.savetxt(tmpf, a, delimiter=',') 
print("Size CSV file", getsize(tmpf.name)) 

tmpf = NamedTemporaryFile() 
np.save(tmpf, a) 
tmpf.seek(0) 
loaded = np.load(tmpf) 
print("Shape", loaded.shape) 
print("Size .npy file", getsize(tmpf.name)) 

df = pd.DataFrame(a) 
df.to_pickle(tmpf.name) 
print("Size pickled dataframe", getsize(tmpf.name)) 
print("DF from pickle\n", pd.read_pickle(tmpf.name)) 

NumPy 提供了一种特定于 NumPy 的格式.npy，可以用来存储 NumPy 数组。在演示这种格式之前，我们将生成一个充满随机值的 365x4 NumPy 数组。该数组模拟一年中四个变量的每日测量值(例如，带有测量温度、湿度、降水和大气压力的传感器的天气数据站)。我们将使用标准的 Python NamedTemporaryFile来存储数据。临时文件应该被自动删除。

将数组存储在 CSV 文件中，并按如下方式检查其大小:

tmpf = NamedTemporaryFile() 
np.savetxt(tmpf, a, delimiter=',') 
print("Size CSV file", getsize(tmpf.name)) 

CSV 文件大小打印如下:

 Size CSV file 36693

以 NumPy .npy格式保存数组，加载数组，检查其形状和.npy文件的大小:

tmpf = NamedTemporaryFile() 
np.save(tmpf, a) 
tmpf.seek(0) 
loaded = np.load(tmpf) 
print("Shape", loaded.shape) 
print("Size .npy file", getsize(tmpf.name)) 

需要调用seek()方法来模拟关闭和重新打开临时文件。形状应以文件大小打印:

Shape (365, 4)
Size .npy file 11760

.npy文件大致比 CSV 文件小三倍，果然不出所料。Python 允许我们存储实际上任意复杂的数据结构。我们也可以将 Pandas 数据框或系列存储为泡菜。

注

Python pickle 是一种用于将 Python 对象存储到磁盘或其他介质的格式。这叫腌制。我们可以从存储中重新创建 Python 对象。这个反向过程叫做拆线(参考docs.python.org/2/library/p…)。腌制已经发展了多年，因此，各种腌制协议存在。不是所有 Python 对象都可以腌制；然而，也有替代的实现，比如 dill、，允许腌制更多类型的 Python 对象。如果可能，使用 cPickle(包含在标准 Python 发行版中)，因为它是用 C 实现的，因此速度更快。

从生成的 NumPy 数组中创建一个 DataFrame，用to_pickle()方法将其写入 pickle，用read_pickle()函数从 pickle 中检索:

df = pd.DataFrame(a) 
df.to_pickle(tmpf.name) 
print("Size pickled dataframe", getsize(tmpf.name)) 
print("DF from pickle\n", pd.read_pickle(tmpf.name)) 

数据框的泡菜略大于.npy文件，如以下打印输出所示:

Size pickled dataframe 12244
DF from pickle
           0         1         2         3
0   0.496714 -0.138264  0.647689  1.523030
[TRUNCATED]
59 -2.025143  0.186454 -0.661786  0.852433
         ...       ...       ...       ...

[365 rows x 4 columns]

用 PyTables 存储数据

分层数据格式 ( HDF )是存储大数值数据的规范和技术。HDF 是在超级计算社区创建的，现在是一个开放标准。HDF 的最新版本是 HDF5 ，也是我们将要使用的版本。HDF5 以组和数据集的形式组织数据。数据集是多维齐次数组。组可以包含其他组或数据集。组就像分层文件系统中的目录。

两个主要的 HDF5 Python 库如下:

• h5y
• iptables(iptables)

在这个例子中，我们将使用 PyTables。PyTables 有许多依赖项:

• 我们在第 1 章、Python 库中安装的 NumPy 包
• numexpr 包声称它计算多运算符数组表达式的速度比 NumPy 快很多倍
• HDF5

注

HDF5 的并行版本也需要 MPI。可以通过从www.hdfgroup.org/HDF5/releas…获得一个发行版并运行以下命令(可能需要几分钟)来安装 HDF5:

$ gunzip < hdf5-X.Y.Z.tar.gz | tar xf -
$ cd hdf5-X.Y.Z
$ make
$ make install

你最喜欢的包管理器很可能有一个 HDF5 的发行版。选择最新的稳定版本。在写这本书的时候，安装的版本是 1.8.12。

第二个依赖项 numexpr 声称能够比 NumPy 更快地执行某些操作。它支持多线程，并且有自己的用 c 实现的虚拟机，在 PyPi 上有 Numexpr 和 PyTables，所以我们可以用pip安装这些，如下所示:

$ pip3 install numexpr tables

同样，我们将生成随机值并用这些随机值填充 NumPy 数组。创建一个 HDF5 文件，并用以下代码将 NumPy 数组附加到根节点:

filename = "pytable_demo.h5" 
h5file = tables.openFile(filename, mode='w', ) 
root = h5file.root 
h5file.createArray(root, "array", a) 
h5file.close() 

读取 HDF5 文件并打印其文件大小:

h5file = tables.openFile(filename, "r") 
print(getsize(filename)) 

我们得到的文件大小值是13824。一旦我们读取了一个 HDF5 文件并获得了它的句柄，我们通常会遍历它来找到我们需要的数据。由于我们只有一个数据集，遍历非常简单。调用iterNodes()和read()方法取回 NumPy 数组:

for node in h5file.iterNodes(h5file.root): 
   b = node.read() 
   print(type(b), b.shape) 

数据集的类型和形状符合我们的预期:

    <class 'numpy.ndarray'> (365, 4)

以下代码可以在本书代码包的ch-05.ipynb文件中找到:

import numpy as np 
import tables 
from os.path import getsize 

np.random.seed(42) 
a = np.random.randn(365, 4) 

filename = "pytable_demo.h5" 
h5file = tables.open_file( filename, mode='w', ) 
root = h5file.root 
h5file.create_array(root, "array", a) 
h5file.close() 

h5file = tables.open_file(filename, "r") 
print(getsize(filename)) 

for node in h5file.root: 
   b = node.read() 
   print(type(b), b.shape) 
 h5file.close() 

读写 Pandas 数据帧到 HDF5 商店

HDFStore类是负责处理 HDF5 数据的 Pandas 抽象。使用随机数据，我们将演示这个功能。

给HDFStore构造函数一个演示文件的路径，并创建一个存储:

filename = "pytable_df_demo.h5"  
store = pd.io.pytables.HDFStore(filename) 
print(store) 

前面的代码片段将打印存储及其内容的文件路径，该路径目前为空:

    <class 'pandas.io.pytables.HDFStore'>
    File path: pytable_df_demo.h5
    Empty

HDFStore有一个类似字典的接口，这意味着我们可以存储值，比如 Pandas 数据帧和相应的查找键。将包含随机数据的数据帧存储在HDFStore中，如下所示:

store['df'] = df 
print(store) 

现在，该存储包含数据，如以下输出所示:

    <class 'pandas.io.pytables.HDFStore'>
    File path: pytable_df_demo.h5
                frame        (shape->[365,4])

我们可以通过三种方式访问数据框:使用get()方法、类似字典的查找或点状访问。让我们试试这个:

print("Get", store.get('df').shape) 
print("Lookup", store['df'].shape) 
print("Dotted", store.df.shape) 

数据框的形状对于所有三种访问方法都是相同的:

Get (365, 4)
Lookup (365, 4)
Dotted (365, 4)

我们可以通过调用remove()方法或使用del运算符来删除商店中的商品。显然，我们只能移除一个项目一次。从存储中删除数据帧:

del store['df'] 
print("After del\n", store) 

商店现在又空了:

After del
<class 'pandas.io.pytables.HDFStore'>
File path: pytable_df_demo.h5
Empty

is_open属性表示店铺是否开门。商店可以用close()方法关闭。关闭商店并检查其是否已关闭:

print("Before close", store.is_open) 
store.close() 
print("After close", store.is_open) 

一旦关闭，商店将不再营业，如下所述:

Before close True
After close False

Pandas 还提供了一个数据框架to_hdf()方法和一个顶层read_hdf()函数来读写 HDF 数据。调用to_hdf()方法并读取数据:

df.to_hdf(tmpf.name, 'data', format='table') 
print(pd.read_hdf(tmpf.name, 'data', where=['index>363'])) 

读写应用编程接口的参数是文件路径、存储中组的标识符和可选的格式字符串。格式可以是固定的，也可以是表格的。固定格式更快，但不能追加或搜索。表格格式对应于 PyTables Table结构，允许搜索和选择。我们在数据框上获得以下查询值:

                0         1         2         3
    364  0.753342  0.381158  1.289753  0.673181
    [1 rows x 4 columns]

本书代码包中的ch-05.ipynb文件包含以下代码:

import numpy as np 
import pandas as pd 

np.random.seed(42) 
a = np.random.randn(365, 4) 

filename = "pytable_df_demo.h5" 
store = pd.io.pytables.HDFStore(filename) 
print(store) 

df = pd.DataFrame(a) 
store['df'] = df 
print(store) 

print("Get", store.get('df').shape) 
print("Lookup", store['df'].shape) 
print( "Dotted", store.df.shape) 

del store['df'] 
print("After del\n", store) 

print("Before close", store.is_open) 
store.close() 
print("After close", store.is_open) 

df.to_hdf('test.h5', 'data', format='table') 
print(pd.read_hdf('test.h5', 'data', where=['index>363']))  

和 Pandas 一起读写 Excel

Excel 文件包含大量重要数据。当然，我们可以以其他更便携的格式导出这些数据，例如 CSV。但是用 Python 读写 Excel 文件更方便。正如在 Python 世界中常见的那样，目前有不止一个项目致力于提供 Excel 输入/输出功能的目标。我们将需要安装的模块，以获得 Excel 输入/输出工作与 Pandas 有点模糊的文件。原因是 Pandas 所依赖的项目是独立的，发展迅速。Pandas 包对它接受为 Excel 文件的文件很挑剔。这些文件必须有.xls或.xlsx后缀，否则，我们会得到以下错误:

ValueError: No engine for filetype: '' 

这很容易解决。例如，如果我们创建一个临时文件，我们只需给它一个适当的后缀。如果您没有安装任何东西，您将收到以下错误消息:

ImportError: No module named openpyxl.workbook

以下命令通过安装openpyxl消除错误:

$ pip3 install openpyxl xlsxwriter xlrd

openpyxl模块是 PHPExcel 的一个端口，支持.xlsx文件的读写。

型式

如果由于某种原因pip install方法对您不起作用，您可以在pythonhosted.org/openpyxl/找到替代安装说明。

读取.xlsx文件也需要xlsxwriter模块。xlrd模块能够从.xls和.xlsx文件中提取数据。

让我们生成随机值来填充 Pandas 数据框，从数据框创建一个 Excel 文件，从 Excel 文件重新创建数据框，并对其应用mean()方法。对于 Excel 文件的工作表，我们可以指定从零开始的索引或名称。

参考书籍代码包中的ch-05.ipynb文件，其中将包含以下代码:

import numpy as np 
import pandas as pd 

np.random.seed(42) 
a = np.random.randn(365, 4) 

filename="excel_demo.xlsx" 
df = pd.DataFrame(a) 
print(filename) 
df.to_excel(filename, sheet_name='Random Data') 
print("Means\n", pd.read_excel(filename, 'Random Data').mean())) 

使用to_excel()方法创建一个 Excel 文件:

df.to_excel(tmpf.name, sheet_name='Random Data') 

使用顶级read_excel()功能重新创建数据框:

print("Means\n", pd.read_excel(tmpf.name, 'Random Data').mean()) 

打印方式如下:

/var/folders/k_/xx_xz6xj0hx627654s3vld440000gn/T/tmpeBEfnO.xlsx
Means
0    0.037860
1    0.024483
2    0.059836
3    0.058417
dtype: float64

使用 REST web 服务和 JSON

具象 状态 转移 ( REST ) web 服务使用 REST 架构风格(更多信息请参考http://en . Wikipedia . org/wiki/具象 _state_transfer )。在 HTTP(S)协议的通常上下文中，我们有 GET 、 POST 、 PUT 和 DELETE 方法。这些方法可以与数据上的常见操作保持一致，以创建、请求、更新或删除数据项。

在 RESTful API 中，数据项由 URIs 标识，如http://example.com/resources或http://example.com/resources/item42。REST 不是一个官方标准，但它如此广泛，我们需要了解它。Web 服务经常使用 JavaScript 对象符号 ( JSON )(更多信息参考en.wikipedia.org/wiki/JSON)来交换数据。在这种格式中，数据是使用 JavaScript 符号编写的。该符号类似于 Python 列表和字典的语法。在 JSON 中，我们可以定义由列表和字典的组合组成的任意复杂的数据。为了说明这一点，我们将使用一个非常简单的 JSON 字符串，它对应于一个字典，给出特定 IP 地址的地理信息:

{"country":"Netherlands","dma_code":"0","timezone":"Europe\/Amsterdam","area_code":"0","ip":"46.19.37.108","asn":"AS196752","continent_code":"EU","isp":"Tilaa V.O.F.","longitude":5.75, "latitude":52.5,"country_code":"NL","country_code3":"NLD"}

以下是来自ch-05.ipynb文件的代码:

import json 

json_str = '{"country":"Netherlands","dma_code":"0","timezone":"Europe\/Amsterdam","area_code":"0","ip":"46.19.37.108","asn":"AS196752","continent_code":"EU","isp":"Tilaa V.O.F.","longitude":5.75,"latitude":52.5,"country_code":"NL","country_code3":"NLD"}' 

data = json.loads(json_str) 
print("Country", data["country"]) 
data["country"] = "Brazil" 
print(json.dumps(data)) 

Python 有一个标准的 JSON API，非常容易使用。使用loads()函数解析 JSON 字符串:

data = json.loads(json_str) 

使用以下代码访问country值:

print "Country", data["country"] 

前一行应打印以下内容:

Country Netherlands

覆盖country值，并根据新的 JSON 数据创建一个字符串:

data["country"] = "Brazil" 
printjson.dumps(data) 

结果是带有新country值的 JSON。这种顺序不会像字典中通常出现的那样被保留:

{"longitude": 5.75, "ip": "46.19.37.108", "isp": "Tilaa V.O.F.", "area_code": "0", "dma_code": "0", "country_code3": "NLD", "continent_code": "EU", "country_code": "NL", "country": "Brazil", "latitude": 52.5, "timezone": "Europe/Amsterdam", "asn": "AS196752"}

与 Pandas 一起读写 JSON

我们可以很容易地从前面例子中的 JSON 字符串创建一个 PandasSeries。Pandasread_json()功能可以创建 Pandas 系列或 Pandas 数据框。

以下示例代码可以在本书的代码包ch-05.ipynb中找到:

import pandas as pd 

json_str = '{"country":"Netherlands","dma_code":"0","timezone":"Europe\/Amsterdam","area_code":"0","ip":"46.19.37.108","asn":"AS196752","continent_code":"EU","isp":"Tilaa V.O.F.","longitude":5.75,"latitude":52.5,"country_code":"NL","country_code3":"NLD"}' 

data = pd.read_json(json_str, typ='series') 
print("Series\n", data) 

data["country"] = "Brazil" 
print("New Series\n", data.to_json()) 

我们可以指定一个 JSON 字符串或者 JSON 文件的路径。调用read_json()函数，从前面例子中的 JSON 字符串创建 PandasSeries:

data = pd.read_json(json_str, typ='series') 
print("Series\n", data) 

在生成的Series中，按键按字母顺序排列:

    Series
    area_code                        0
    asn                       AS196752
    continent_code                  EU
    country                Netherlands
    country_code                    NL
    country_code3                  NLD
    dma_code                         0
    ip                    46.19.37.108
    ispTilaa V.O.F.
    latitude                      52.5
    longitude                     5.75
    timezone          Europe/Amsterdam
    dtype: object

再次更改country值，并使用to_json()方法将 PandasSeries转换为 JSON 字符串:

data["country"] = "Brazil" 
print("New Series\n", data.to_json()) 

在新的 JSON 字符串中，键顺序被保留，但是我们还有一个不同的country值:

 New Series
    {"area_code":"0","asn":"AS196752","continent_code":"EU","country":"Brazil","country_code":"NL","country_code3":"NLD","dma_code":"0","ip":"46.19.37.108","isp":"Tilaa V.O.F.","latitude":52.5,"longitude":5.75,"timezone":"Europe\/Amsterdam"}

解析 RSS 和 Atom 提要

真正简单的联合 ( RSS )和 Atom feeds(参考en.wikipedia.org/wiki/RSS)经常用于博客和新闻。这些类型的订阅源遵循发布/订阅模式。例如，帕克特出版公司有一个文章和书籍公告的 RSS 源。我们可以订阅订阅源以获得及时的更新。Python feedparser模块允许我们轻松解析 RSS 和 Atom 提要，而无需处理很多技术细节。feedparser模块可以安装pip如下:

$ pip3 install feedparser

解析完一个 RSS 文件后，我们可以使用虚线符号访问底层数据。解析打包发布 RSS 源并打印条目数:

import feedparser as fp 

rss = fp.parse("http://www.packtpub.com/rss.xml") 

print("# Entries", len(rss.entries)) 

打印条目的数量(数量可能因每次程序运行而异):

# Entries 10

如果条目包含带有以下代码的单词Python，则打印条目标题和摘要:

for i, entry in enumerate(rss.entries): 
   if "Python" in entry.summary: 
      print(i, entry.title) 
      print(entry.summary) 

在此特定运行中，打印了以下内容(如果您自己尝试，则可能会得到其他内容，如果过滤器限制太多，则可能什么也没有):

42 Create interactive plots with matplotlib using Pack't new book and eBook
About the author: Alexandre Devert is a scientist. He is an enthusiastic Python coder as well and never gets enough of it! He used to teach data mining, software engineering, and research in numerical optimization.
Matplotlib is part of the Scientific Python modules collection. It provides a large library of customizable plots and a comprehensive set of backends. It tries to make easy things easy and make hard things possible. It can help users generate plots, add dimensions to plots, and also make plots interactive with just a few lines of code. Also, matplotlib integrates well with all common GUI modules.

以下代码可以在本书代码包的ch-05.ipynb文件中找到:

import feedparser as fp 

rss = fp.parse("http://www.packtpub.com/rss.xml") 

print("# Entries", len(rss.entries)) 

for i, entry in enumerate(rss.entries): 
   if "Java" in entry.summary: 
      print(i, entry.title) 
      print(entry.summary) 

用靓汤解析 HTML

超文本标记语言 ( HTML )是用来创建网页的基础技术。HTML 是由 HTML 元素组成的，这些元素由所谓的标记组成，这些标记封装在倾斜的括号中(例如，<html>)。通常，标签与一个开始和结束标签成对出现在一个分层的树状结构中。伯纳斯-李于 1991 年首次发布了一个与 HTML 相关的规范草案。最初，只有 18 个 HTML 元素。正式的 HTML 定义由互联网工程任务组 ( IETF )于 1993 年发布。IETF 在 1995 年完成了 HTML 2.0 标准。大约在 2013 年，指定了最新的 HTML 版本 HTML5。与 XHTML 和 XML 相比，HTML 并不是一个非常严格的标准。

现代浏览器容忍大量违反标准的行为，使得网页成为一种非结构化数据。例如，我们可以将 HTML 视为一个大字符串，并用正则表达式对其执行字符串操作。这种方法只适用于简单的项目。

我曾在专业的环境中从事过网页抓取项目，所以从个人经验来看，我可以告诉你，我们需要更复杂的方法。在现实场景中，可能需要以编程方式提交 HTML 表单，例如，登录、浏览页面和健壮地管理 cookies。从网络上抓取数据的问题是，如果我们不能完全控制我们正在抓取的网页，我们可能不得不经常更改代码。网站所有者也可能主动阻止程序访问，甚至可能是非法的。出于这些原因，您应该总是先尝试使用其他替代方法，例如 REST API。

如果您必须通过抓取来检索数据，建议您使用 Python 美丽汤应用编程接口。这个应用编程接口可以从 HTML 和 XML 文件中提取数据。新项目应该用美人汤 4，因为美人汤 3 已经不开发了。我们可以用以下命令安装美人汤 4(类似easy_install):

$ pip3 install beautifulsoup4 lxml

如果这不起作用，您可以简单地将美丽的汤和您自己的代码打包在一起。为了演示解析 HTML，我用来自 loripsum.net/ T2 的生成器生成了本书代码包中的 T0 文件。然后，我稍微编辑了一下文件。该文件的内容是西塞罗用拉丁语创作的公元前一世纪的文本，这是创建网站模型的传统方式。请参考以下网页顶部的截图:

在这个例子中，我们将使用美丽的汤 4 和标准的 Python 正则表达式库。

使用以下行导入这些库:

from bs4 import BeautifulSoup 
import re 

打开 HTML 文件，创建一个BeautifulSoup对象，如下行:

soup = BeautifulSoup(open('loremIpsum.html')) 

使用虚线符号，我们可以访问第一个<div>元素。<div> HTML 元素用于组织元素和设置元素的样式。访问第一个div元素，如下所示:

print("First div\n", soup.div) 

结果输出是一个 HTML 片段，带有第一个<div>标签和它包含的所有标签:

    First div
    <div class="tile">
    <h4>Development</h4>
         0.10.1 - July 2014<br/>
    </div>

注

这个特殊的div元素有一个值为tile的类属性。类属性属于应用于这个div元素的 CSS 样式。层叠样式表 ( CSS )是一种用于设计网页元素样式的语言。CSS 是一个广泛使用的规范，它通过 CSS 类来处理网页的外观和感觉。CSS 通过定义颜色、字体和元素布局来帮助分离内容和表示。这种分离导致更简单和更干净的设计。

标签的属性可以以类似字典的方式访问。如下打印<div>标签的类属性值:

print("First div class", soup.div['class']) 
First div class ['tile'] 

虚线符号允许我们访问任意深度的元素。例如，打印第一个<dfn>标签的文本如下:

print("First dfn text", soup.dl.dt.dfn.text) 

打印一行拉丁文文本(Solisten, I pray):

    First dfn text Quareattende, quaeso

有时候，我们只对 HTML 文档的超链接感兴趣。例如，我们可能只想知道哪个文档链接到哪个其他文档。在 HTML 中，链接用<a>标记指定。该标签的href属性保存链接指向的网址。BeautifulSoup班有个很好用的find_all()方法，我们会经常用到。使用find_all()方法定位所有超链接:

for link in soup.find_all('a'): 
      print("Link text", link.string, "URL", link.get('href')) 

文档中有三个链接具有相同的网址，但具有三种不同的文本:

Link text loripsum.net URL http://loripsum.net/ 
Link text Potera tautem inpune; URL http://loripsum.net/ 
Link text Is es profecto tu. URL http://loripsum.net/ 

我们可以省略find_all()方法作为捷径。访问所有<div>标签的内容，如下所示:

for i, div in enumerate(soup('div')): 
   print(i, div.contents) 

contents属性保存一个包含 HTML 元素的列表:

    0 [u'\n', <h4>Development</h4>, u'\n     0.10.1 - July 2014', <br/>, u'\n']
    1 [u'\n', <h4>Official Release</h4>, u'\n     0.10.0 June 2014', <br/>, u'\n']
    2 [u'\n', <h4>Previous Release</h4>, u'\n     0.09.1 June 2013', <br/>, u'\n']

带有唯一标识的标签很容易找到。选择official标识的<div>元素，打印第三个元素:

official_div = soup.find_all("div", id="official") 
print("Official Version", official_div[0].contents[2].strip()) 

许多网页是根据访问者的输入或外部数据动态创建的。网上购物网站的大部分内容都是这样提供的。如果我们在处理一个动态网站，我们必须记住，任何标签属性值都可能在一瞬间发生变化。通常，在大型网站中，会自动生成标识，从而产生长的字母数字字符串。最好不要寻找完全匹配，而是使用正则表达式。稍后我们将看到一个基于模式的匹配示例。前面的代码片段打印了一个版本号和月份，您可以在软件产品的网站上找到:

    Official Version 0.10.0 June 2014

众所周知，class是 Python 的一个关键词。要查询标签中的class属性，我们将其与class_进行匹配。获取已定义类别属性的<div>标签数量:

print("# elements with class",len(soup.find_all(class_=True))) 

正如预期的那样，我们发现了三个标签:

    # elements with class 3

用类别"tile"找到<div>标签的数量:

tile_class = soup.find_all("div", class_="tile") 
print("# Tile classes", len(tile_class)) 

有两个<div>类标签tile，一个<div>类标签notile:

  # Tile classes 2

定义一个匹配所有<div>标签的正则表达式:

print("# Divs with class containing tile", len(soup.find_all("div", class_=re.compile("tile")))) 

再次发现三种情况:

   # Divs with class containing tile 3

在 CSS 中，我们可以定义模式来匹配元素。这些模式被称为 CSS 选择器，并记录在www.w3.org/TR/selector…中。我们也可以用 CSS 选择器从BeautifulSoup类中选择元素。使用select()方法将<div>元素与notile类匹配:

print("Using CSS selector\n", soup.select('div.notile')) 

屏幕上会显示以下内容:

Using CSS selector
[<div class="notile">
<h4>Previous Release</h4>
 0.09.1 June 2013<br/>
</div>]

一个 HTML 排序的列表看起来像一个编号的项目符号列表。有序列表由每个列表项的一个<ol>标签和几个<li>标签组成。select()方法的结果可以被切片为任何 Python 列表。请参考下面的订购列表截图:

选择有序列表中的前两个列表项:

print("Selecting ordered list list items\n", soup.select("ol > li")[:2]) 

显示了以下两个列表项:

Selecting ordered list list items
[<li>Cur id non ita fit?</li>, <li>In qua si nihil est praeter rationem, sit in una virtute finis bonorum;</li>]

在 CSS 选择器迷你语言中，我们从1开始计数。选择第二个列表项，如下所示:

print("Second list item in ordered list", soup.select("ol>li:nth-of-type(2)")) 

第二个清单项可以翻译成英文其中，如果没有什么违背理智的事情，就让他成为美好事物终结的力量于一身:

Second list item in ordered list [<li>In qua si nihil est praeter rationem, sit in una virtute finis bonorum;</li>]

如果我们在浏览器中查看网页，我们可能会决定检索与某个正则表达式匹配的文本节点。查找包含带有text属性的字符串2014的所有文本节点:

print("Searching for text string", soup.find_all(text=re.compile("2014"))) 

这会打印以下文本节点:

Searching for text string [u'\n     0.10.1 - July 2014', u'\n     0.10.0 June 2014']

这只是BeautifulSoup课能为我们做什么的一个简单概述。美丽汤也可以用来修改 HTML 或 XML 文档。它有实用工具来排除故障，漂亮的打印，并处理不同的字符集。代码请参考ch-05.ipynb:

from bs4 import BeautifulSoup 
import re 

soup = BeautifulSoup(open('loremIpsum.html'),"lxml") 

print("First div\n", soup.div) 
print("First div class", soup.div['class']) 

print("First dfn text", soup.dl.dt.dfn.text) 

for link in soup.find_all('a'): 
   print("Link text", link.string, "URL", link.get('href')) 

# Omitting find_all 
for i, div in enumerate(soup('div')): 
   print(i, div.contents) 

#Div with id=official 
official_div = soup.find_all("div", id="official") 
print("Official Version", official_div[0].contents[2].strip()) 

print("# elements with class", len(soup.find_all(class_=True))) 

tile_class = soup.find_all("div", class_="tile") 
print("# Tile classes", len(tile_class)) 

print("# Divs with class containing tile", len(soup.find_all("div", class_=re.compile("tile")))) 

print("Using CSS selector\n", soup.select('div.notile')) 
print("Selecting ordered list list items\n", soup.select("ol > li")[:2]) 
print("Second list item in ordered list", soup.select("ol > li:nth-of-type(2)")) 

print("Searching for text string", soup.find_all(text=re.compile("2014"))) 

鼓励读者阅读参考部分提到的书籍，了解更多关于美丽汤功能的详细信息，例如搜索子节点的返回节点，获取返回节点的第 n 个父节点，获取返回节点的第 n 个同级节点，以及其他高级功能。

总结

在本章中，我们学习了以不同格式检索、处理和存储数据。这些是 CSV、NumPy .npy格式、Python pickle、JSON、RSS 和 HTML 格式。我们使用了 NumPy Pandas、JSON、feedparser 和美丽的汤库。

下一章第六章数据可视化，讲的是用 Python 可视化数据的重要话题。可视化是我们开始分析数据时经常做的事情。它有助于显示数据中变量之间的关系。通过可视化数据，我们还可以了解其统计特性。

参考

动词 （verb 的缩写）g .奈尔，*美丽汤入门，*帕克特出版，2014。

六、数据可视化

数据分析的第一步是可视化。即使是在查看数值表时，我们也可以在脑海中形成一个图表化数据的样子。数据可视化涉及信息的可视化表示的概念和分析，表示以某种形式模式抽象的数据，包括数据测量单位的属性或数量。数据可视化与科学可视化和统计图形密切相关。Python matplotlib库是一个著名的基于 NumPy 的绘图库，我们将在本章中使用它。它有一个面向对象的和类似于 Matlab 的程序性应用编程接口，可以并行使用。在matplotlib.org/gallery.htm…可以找到一个包含 matplotlib 示例的图库，这是一个基于云的数据可视化服务。我们将在本章末尾简要介绍这项服务。以下是本章将涵盖的主题列表:

• matplotlib子包
• 基本matplotlib地块
• 对数图
• 散点图
• 图例和注释
• 三维图
• Pandas 中的阴谋
• 滞后图
• 自相关图
• plot.ly

matplotlib 子包

如果我们更改ch-01.ipynb最后一节的代码来列出matplotlib子包，我们会得到以下结果:

matplotlib version 1.3.1
matplotlib.axes
matplotlib.backends
matplotlib.compat
matplotlib.delaunay DESCRIPTION :Author: Robert Kern   
    <robert.kern@gmail.com> :Copyright: Copyright 2005  
    Robert Kern. 
    :License: BSD-style license. See LICENSE.tx
matplotlib.projections
matplotlib.sphinxext
matplotlib.style
matplotlib.testing
matplotlib.tests
matplotlib.tri

子包的名称很容易解释。后端是指输出最终结果的方式。该输出可以是几种文件格式之一，也可以显示在图形用户界面的屏幕上。

基本 matplotlib 打印

我们在第 1 章、Python 库中安装了 matplotlib 和 IPython。如果你需要刷新记忆，你可以回到那一章。许多人认为类似于 Matlab 的程序性 matplotlib 应用编程接口比面向对象的应用编程接口更容易使用，所以我们将首先演示这个程序性应用编程接口。为了在 matplotlib 中创建一个非常基本的图，我们需要调用matplotlib.pyplot子包中的plot()函数。该函数为具有已知x和y 坐标的单个点列表或多个点列表生成二维图。

或者，我们可以传递一个格式参数，例如，指定一个虚线样式。plot()功能的格式选项和参数列表很长，但是使用以下命令很容易查找(在导入matplotlib.pyplot库之后):

In [1]: help(plot)

在本例中，我们将绘制两条线-一条是实线样式(默认)，另一条是虚线样式。

以下演示代码位于本书代码包的ch-06.ipynb文件中:

import matplotlib.pyplot as plt 
import numpy as np 

x = np.linspace(0, 20) 

plt.plot(x,  .5 + x) 
plt.plot(x, 1 + 2 * x, '--') 
plt.show() 

请按照以下步骤绘制上述线条:

1. 首先，我们用 NumPy linspace()函数指定x坐标。指定0的开始值和20的结束值:

           x = np.linspace(0, 20) 
   
   

2. 如下绘制线条:

           plt.plot(x,  .5 + x) 
           plt.plot(x, 1 + 2 * x, '--') 
   
   

3. At this juncture, we can either save the plot to a file with the savefig() function or show the plot on the screen with the show() function. Show the plot on the screen as follows:

           plt.show() 
   
   

   最终结果见下图:

   

对数图

对数图(或对数图)是使用对数标度的图。对数刻度显示变量的值，它使用与数量级匹配的间隔，而不是规则的线性刻度。对数图有两种类型。对数-对数图在两个轴上都采用对数标度，在 matplotlib 中由matplotlib.pyplot.loglog()函数表示。半对数图在一个轴上使用线性标度，在另一个轴上使用对数标度。这些图在 matplotlib 应用编程接口中由semilogx()和semilogy()函数表示。在对数-对数图上，幂律显示为直线。在半对数图上，直线代表指数规律。

摩尔定律就是这样的定律。这不是物理上的，更多的是经验上的观察。戈登·摩尔发现了集成电路中晶体管数量每两年翻一番的趋势。在http://en . Wikipedia . org/wiki/Transistor _ count #微处理器上，你可以看到一张表格，上面有各种微处理器的晶体管计数和相应的推出年份。

从这个表中，我准备了一个 CSV 文件，transcount.csv，只包含晶体管计数和年份。我们仍然需要平均每年的晶体管数量。平均和加载可以用 Pandas 来完成。如有需要，可参考第三章、Pandas 底火、获取提示。一旦我们在表中有了每年的平均晶体管计数，我们就可以试着在计数与年份的对数上画一条直线。NumPy polyfit()函数允许我们将数据拟合为多项式。

以下代码请参考本书代码包中的ch-06.ipynb文件:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 
years = df.index.values 
counts = df['trans_count'].values 
poly = np.polyfit(years, np.log(counts), deg=1) 
print("Poly", poly) 
plt.semilogy(years, counts, 'o') 
plt.semilogy(years, np.exp(np.polyval(poly, years))) 
plt.show() 

以下步骤将解释前面的代码:

1. 将数据拟合如下:

           poly = np.polyfit(years, np.log(counts), deg=1) 
           print("Poly", poly) 
   
   

2. 拟合的结果是一个多项式对象(见[。这个对象的字符串表示给出了多项式系数，系数的阶数是递减的，所以最高的系数排在第一位。对于我们的数据，我们获得以下多项式系数:

           Poly [  3.61559210e-01  -7.05783195e+02]
   ```](http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.polynomial.Polynomial.html#numpy.polynomial.polynomial.Polynomial) 
   

3. The NumPy polyval() function enables us to evaluate the polynomial we just obtained. Plot the data and fit it with the semilogy() function:

           plt.semilogy(years, counts, 'o') 
           plt.semilogy(years, np.exp(np.polyval(poly, years))) 
   
   

   趋势线绘制为实线，数据点绘制为实心圆。最终结果见下图:

   

散点图

一个散点图显示了笛卡尔坐标系中两个变量之间的关系。每个数据点的位置由这两个变量的值决定。散点图可以提示所研究的变量之间的任何相关性。上升趋势表明正相关。一张气泡图是散点图的延伸。在气泡图中，第三个变量的值由数据点周围气泡的相对大小来表示，因此得名。

在en.wikipedia.org/wiki/Transi…有一张表格，上面有图形处理器单元 ( 图形处理器)的晶体管计数。

图形处理器是用于高效显示图形的专用电路。由于现代显示硬件的工作方式，图形处理器可以用高度并行的操作来处理数据。图形处理器是计算领域的新发展。在本书代码包的gpu_transcount.csv文件中，你会注意到我们没有太多的数据点。处理缺失数据是一个反复出现的气泡图问题。我们将为缺失值定义默认气泡大小。同样，我们将每年加载和平均数据。然后，我们将把年指数上的中央处理器和图形处理器数据帧的晶体管计数与外部连接合并。NaN值将被设置为0(本例适用，但有时将NaN值设置为0可能不是一个好主意)。前文中描述的所有功能都包含在第 3 章、**Pandas 初级套装中，因此如果您需要，请参考本章。matplotlib 应用编程接口为散点图和气泡图提供scatter()功能。我们可以使用以下命令查看该函数的文档:

$ ipython3
ln [1]: import matplotlib as mpl
In [2]: help(mpl.scatter)

在本例中，我们将指定s参数，该参数与气泡的大小有关。c参数指定颜色。不幸的是，你在这本书里看不到颜色，所以你必须自己运行例子来看不同的颜色。alpha参数决定了图中气泡的透明度。该值在0(完全透明)和1(不透明)之间变化。按如下方式创建气泡图:

plt.scatter(years, cnt_log, c= 200 * years, s=20 + 200 * gpu_counts/gpu_counts.max(), alpha=0.5) 

本例的以下代码也可以在本书代码包的ch-06.ipynb文件中找到:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 

gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 
print(df) 
years = df.index.values 
counts = df['trans_count'].values 
gpu_counts = df['gpu_trans_count'].values 
cnt_log = np.log(counts) 
plt.scatter(years, cnt_log, c= 200 * years, s=20 + 200 * gpu_counts/gpu_counts.max(), alpha=0.5) 
plt.show() 

最终结果见下图:

传说和注释

图例和注释是一目了然地显示理解情节所需信息的有效工具。典型图将包含以下附加信息元素:

• 描述剧情中各种数据系列的图例。这是通过调用 matplotlib legend()函数并为每个数据系列提供标签来实现的。
• 情节中重要点的注释。matplotlib annotate()功能可用于此目的。matplotlib 注释由一个标签和一个箭头组成。这个函数有很多参数描述标签和箭头的样式和位置，所以可能需要调用help(annotate)进行详细描述。
• 水平轴和垂直轴上的标签。这些标签可以通过xlabel()和ylabel()功能绘制。我们需要给这些函数一个字符串形式的标签文本，以及可选的参数，比如标签的字体大小。
• 具有 matplotlib title()功能的图形的描述性标题。通常，我们只给这个函数一个表示标题的字符串。
• 有一个网格也很好，这样可以很容易地定位点。matplotlib grid()功能打开和关闭绘图网格。

我们将修改上一个示例中的气泡图代码，并添加本章第二个示例中的直线拟合。在此设置中，将标签添加到数据系列，如下所示:

plt.plot(years, np.polyval(poly, years), label='Fit') 
plt.scatter(years, cnt_log, c= 200 * years, s=20 + 200 * gpu_counts/gpu_counts.max(), alpha=0.5, label="Scatter Plot") 

让我们注释一下数据集中的第一个 GPU。为此，获取相关的点，定义注释的标签，指定箭头的样式(arrowprops参数)，并确保注释悬停在有问题的点上方:

gpu_start = gpu.index.values.min() 
y_ann = np.log(df.at[gpu_start, 'trans_count']) 
ann_str = "First GPU\n %d" % gpu_start 
plt.annotate(ann_str, xy=(gpu_start, y_ann), arrowprops=dict(arrowstyle="->"), xytext=(-30, +70), textcoords='offset points') 

完整的代码示例在本书代码包的ch-06.ipynb文件中:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 
gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 
years = df.index.values 
counts = df['trans_count'].values 
gpu_counts = df['gpu_trans_count'].values 

poly = np.polyfit(years, np.log(counts), deg=1) 
plt.plot(years, np.polyval(poly, years), label='Fit') 

gpu_start = gpu.index.values.min() 
y_ann = np.log(df.at[gpu_start, 'trans_count']) 
ann_str = "First GPU\n %d" % gpu_start 
plt.annotate(ann_str, xy=(gpu_start, y_ann), arrowprops=dict(arrowstyle="->"), xytext=(-30, +70), textcoords='offset points') 

cnt_log = np.log(counts) 
plt.scatter(years, cnt_log, c= 200 * years, s=20 + 200 * gpu_counts/gpu_counts.max(), alpha=0.5, label="Scatter Plot") 
plt.legend(loc='upper left') 
plt.grid() 
plt.xlabel("Year") 
plt.ylabel("Log Transistor Counts", fontsize=16) 
plt.title("Moore's Law & Transistor Counts") 
plt.show() 

最终结果见下图:

三维图

二维图是数据可视化的基础。然而，如果你想炫耀，没有什么比得上一个好的三维剧情。我曾经负责过一个可以画等高线图和三维图的软件包。该软件甚至可以绘制出当用特殊的眼镜观看时会突然出现在你面前的图表。

matplotlib API 具有用于三维绘图的Axes3D类。通过演示这个类是如何工作的，我们也将展示面向对象的 matplotlib API 是如何工作的。matplotlib Figure类是图表元素的顶级容器:

1. 创建一个figure对象，如下所示:

           fig = plt.figure() 
   
   

2. 从figure对象创建一个Axes3D对象:

           ax = Axes3D(fig) 
   
   

3. 年份和 CPU 晶体管计数将是我们的X和Y轴。我们有必要根据年份和中央处理器晶体管计数数组创建坐标矩阵。使用meshgrid()功能

           X, Y = np.meshgrid(X, Y) 
   
   

   创建坐标矩阵

4. 用Axes3D类的plot_surface()方法绘制数据:

           ax.plot_surface(X, Y, Z) 
   
   

5. 面向对象的 API 方法的命名约定是以set_开始，以过程对应的函数名结束，如下面的代码片段所示:

           ax.set_xlabel('Year') 
           ax.set_ylabel('Log CPU transistor counts') 
           ax.set_zlabel('Log GPU transistor counts') 
           ax.set_title("Moore's Law & Transistor Counts") 
   
   

您还可以查看本书代码包中ch-06.ipynb文件中的以下代码:

from mpl_toolkits.mplot3d.axes3d import Axes3D 
import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 
gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 

fig = plt.figure() 
ax = Axes3D(fig) 
X = df.index.values 
Y = np.where(df['trans_count'].values>0,np.ma.log(df['trans_count'].values), 0) 
X, Y = np.meshgrid(X, Y) 
Z = np.where(df['gpu_trans_count'].values>0,np.ma.log(df['gpu_trans_count'].values), 0) 
ax.plot_surface(X, Y, Z) 
ax.set_xlabel('Year') 
ax.set_ylabel('Log CPU transistor counts') 
ax.set_zlabel('Log GPU transistor counts') 
ax.set_title("Moore's Law & Transistor Counts") 
plt.show() 

最终结果见下图:

在 Pandas 中绘图

PandasSeries和DataFrame类中的plot()方法包装了相关的 matplotlib 函数。在最基本的形式中，没有任何参数，plot()方法显示了我们在本章中一直使用的数据集的以下图:

要创建半对数图，添加logy参数:

df.plot(logy=True) 

这导致我们的数据如下图所示:

要创建散点图，请将kind参数指定为scatter。我们还需要指定两列。将loglog参数设置为True以生成对数图(对于此代码，我们至少需要 Pandas v0.13.0):

df[df['gpu_trans_count'] > 0].plot(kind='scatter', x='trans_count', y='gpu_trans_count', loglog=True) 

最终结果见下图:

以下程序在本书代码包的ch-06.ipynb文件中:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 

gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 
df.plot() 
df.plot(logy=True) 
df[df['gpu_trans_count'] > 0].plot(kind='scatter', x='trans_count', y='gpu_trans_count', loglog=True) 
plt.show() 

滞后图

一个滞后图是一个时间序列和相同数据滞后的散点图。例如，通过这样的图，我们可以检查今年的中央处理器晶体管计数和前一年的计数之间是否存在可能的相关性。pandas.tools.plotting中的lag_plot()Pandas 功能可以画出滞后图。用1的默认滞后绘制一个 CPU 晶体管计数的滞后图，如下所示:

lag_plot(np.log(df['trans_count'])) 

最终结果见下图:

滞后图示例的以下代码也可以在本书代码包的ch-06.ipynb文件中找到:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 
from pandas.tools.plotting import lag_plot 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 

gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 
lag_plot(np.log(df['trans_count'])) 
plt.show() 

自相关图

自相关图绘制不同时滞的时间序列数据的自相关图。通俗地说，自相关就是时间 n 的值和时间 n 的值+l 的相关性，其中 l 是时滞。通常，这些图用于检查时间序列在其进展中是否具有随机性。在随机时间序列的情况下，自相关在所有时滞分离中都接近于零，在非随机时间序列的一些或所有时滞分离中具有非零的显著性值。我们在第 7 章、信号处理和时间序列中进一步解释自相关。

pandas.tools.plotting中的autocorrelation_plot()Pandas 函数可以画出自相关图。以下是本书代码包中ch-06.ipynb文件的代码:

import matplotlib.pyplot as plt 
import numpy as np 
import pandas as pd 
from pandas.tools.plotting import autocorrelation_plot 

df = pd.read_csv('transcount.csv') 
df = df.groupby('year').aggregate(np.mean) 

gpu = pd.read_csv('gpu_transcount.csv') 
gpu = gpu.groupby('year').aggregate(np.mean) 

df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True) 
df = df.replace(np.nan, 0) 
autocorrelation_plot(np.log(df['trans_count'])) 
plt.show() 

绘制中央处理器晶体管计数的自相关图，如下所示:

autocorrelation_plot(np.log(df['trans_count'])) 

最终结果见下图。正如我们在下图中看到的，较新的值(较小的滞后)与当前值的相关性比较旧的值(较大的滞后)更强，并且当滞后非常大时，相关性衰减到0:

Plot.ly

Plot.ly 是在线数据可视化工具的软件即服务 ( SaaS )云服务。Plot.ly 提供了一个相关的 python 库，可以在用户的机器上与 Python 一起使用。我们可以通过网络界面导入和分析数据，或者完全在本地环境中工作，并在 Plot.ly 网站上发布最终结果。情节可以很容易地在一个团队的网站上分享，允许合作，这确实是网站的第一点。在本节中，我们将给出一个如何用 Python API 绘制方框图的例子。

箱线图是使用四分位数可视化数据集的一种特殊方式。如果我们将排序后的数据集分成四个相等的部分，第一个四分位数将是具有最小数字的部分的最大值。第二个四分位数将是数据集中间的值，也称为中位数。第三个四分位数是中间值和最高值之间的值。箱线图的底部和顶部由第一个和第三个四分位数组成。穿过盒子的线是中线。盒子两端的胡须通常是数据集的最小值和最大值。在这一节的最后，我们将看到一个带注释的方框图，它将澄清问题。使用以下命令安装绘图应用编程接口:

$ sudo pip3 install plotly

安装完应用编程接口后，注册获取应用编程接口密钥。以下代码片段在提供有效密钥后让您登录:

# Change the user and api_key to your own username and api_key
py.sign_in('username', 'api_key')

使用绘图应用编程接口创建方框图，如下所示:

data = Data([Box(y=counts), Box(y=gpu_counts)])
plot_url = py.plot(data, filename='moore-law-scatter')

请参考本书代码包中ch-06.ipynb文件的以下代码:

import plotly.plotly as py
from plotly.graph_objs import *
import numpy as np
import pandas as pd

df = pd.read_csv('transcount.csv')
df = df.groupby('year').aggregate(np.mean)

gpu = pd.read_csv('gpu_transcount.csv')
gpu = gpu.groupby('year').aggregate(np.mean)
df = pd.merge(df, gpu, how='outer', left_index=True, right_index=True)
df = df.replace(np.nan, 0)

# Change the user and api_key to your own username and api_key
py.sign_in('username', 'api_key')

counts = np.log(df['trans_count'].values)
gpu_counts = np.log(df['gpu_trans_count'].values)

data = Data([Box(y=counts), Box(y=gpu_counts)])
plot_url = py.plot(data, filename='moore-law-scatter')
print(plot_url)

最终结果见下图:

总结

在本章中，我们讨论了使用 Python 使用绘图来可视化数据。为此，我们使用了 matplotlib 和 Pandas。我们介绍了箱线图、散点图、气泡图、对数图、自相关图、滞后图、三维图、图例和注释。

对数图(或对数图)是使用对数刻度的图。半对数图在一个轴上使用线性标度，在另一个轴上使用对数标度。散点图将两个变量相对绘制。气泡图是一种特殊类型的散点图。在气泡图中，第三个变量的值相对由数据点周围气泡的大小来表示。自相关绘制不同滞后时间序列数据的自相关图。

我们了解了 plot.ly，一个基于云的在线数据可视化服务，并使用该服务构建了一个箱式图。箱线图基于数据的四分位数可视化数据。

下一章第七章信号处理与时间序列讲述的是一种特殊类型的数据——时间序列。时间序列由已打上时间戳的有序数据点组成。我们测量的很多物理世界数据都是以时间序列的形式出现的，可以认为是一种信号，比如声音、光或电信号。你将在下一章学习如何过滤信号和建模时间序列。*

七、信号处理和时间序列

信号处理是工程和应用数学的一个领域，包括分析随时间变化的变量，这种数据也被称为模拟和数字信号。信号处理技术的一个类别是时间序列分析。一个时间序列是一个有序的数据点列表，从最早的测量值开始。数据点通常是等距的，例如，每小时、每天、每周、每月或每年取样。在时间序列分析中，值的顺序很重要。在同一时间序列中，试图导出一个值与另一个数据点或数据点组合(过去的固定周期数)之间的关系是很常见的。

本章中的时间序列例子使用了年太阳黑子周期数据。该数据由statsmodels包(一个开源 Python 项目)提供。示例使用了 NumPy/SciPy、Pandas 以及statsmodels。

我们将在本章中讨论以下主题:

• statsmodels 模块
• 移动平均线
• 窗口功能
• 定义协整
• 自相关
• 自回归模型
• ARMA 模型
• 产生周期性信号
• 傅里叶分析
• 光谱分析
• 过滤

stats models 模块

要安装statsmodels，执行以下命令:

$ pip3 install statsmodels

在所附的ch-07.ipynb文件中，我们列出了statsmodels模块，得到如下结果:

statmodels version 0.6.1
statsmodels.base 
statsmodels.compat 
statsmodels.datasets 
statsmodels.discrete 
statsmodels.distributions 
statsmodels.duration 
statsmodels.emplike 
statsmodels.formula 
statsmodels.genmod 
statsmodels.graphics 
statsmodels.interface 
statsmodels.iolib 
statsmodels.miscmodels 
statsmodels.nonparametric 
DESCRIPTION For an overview of this module, see docs/source/nonparametric.rst PACKAGE CONTENTS _kernel_base _smoothers_lowess api bandwidths statsmodels.regression 
statsmodels.resampling 
statsmodels.robust 
statsmodels.sandbox 
statsmodels.stats 
statsmodels.tests 
statsmodels.tools 
statsmodels.tsa 

移动平均线

移动平均线常用于分析时间序列。移动平均线指定了以前看到的数据窗口，每次窗口向前滑动一个周期时，移动平均线就会取平均值:

不同类型的移动平均线在用于平均的权重上有本质的不同。例如，指数移动平均线的权重随时间呈指数递减:

这意味着旧值比新值的影响小，这有时是可取的。

本书代码包中ch-07.ipynb文件的以下代码绘制了 11 年和 22 年太阳黑子周期的简单移动平均线:

import matplotlib.pyplot as plt 
import statsmodels.api as sm 
from pandas.stats.moments import rolling_mean 

data_loader = sm.datasets.sunspots.load_pandas() 
df = data_loader.data 
year_range = df["YEAR"].values 
plt.plot(year_range, df["SUNACTIVITY"].values, label="Original") 
plt.plot(year_range, df.rolling(window=11).mean()["SUNACTIVITY"].values, label="SMA 11") 
plt.plot(year_range, df.rolling(window=22).mean()["SUNACTIVITY"].values, label="SMA 22") 
plt.legend() 
plt.show() 

我们可以表示指数移动平均线的指数递减权重策略，如下面的 NumPy 代码所示:

weights = np.exp(np.linspace(-1., 0., N)) 
weights /= weights.sum() 

简单的移动平均线使用相等的权重，在代码中看起来如下:

def sma(arr, n): 
   weights = np.ones(n) / n 

   return np.convolve(weights, arr)[n-1:-n+1]  

因为我们可以将数据加载到Pandas 数据框中，所以使用 PandasDataFrame.rolling().mean()功能更方便。使用statsmodels加载数据如下:

data_loader = sm.datasets.sunspots.load_pandas() 
df = data_loader.data 

最终结果见下图:

窗口功能

NumPy 有许多窗口例程，可以像我们在上一节中所做的那样，在滚动窗口中计算权重。

一个窗口函数是一个在一个区间(窗口)内定义的函数，或者是零值。我们可以使用窗口函数进行光谱分析和滤波器设计(更多背景信息，请参考en.wikipedia.org/wiki/Window…)。棚车车窗为矩形车窗，公式如下:

w(n) = 1 

三角窗形状像三角形，公式如下:

在上式中，L可以等于N、N+1或N-1。在最后一种情况下，窗口函数被称为巴特利特窗口。布莱克曼窗为钟形，定义如下:

汉宁窗也是钟形的，定义如下:

在 Pandas 应用编程接口中，DataFrame.rolling()函数提供相同的功能，不同的窗口函数对应不同的win_type字符串参数值。另一个参数是窗口的大小，对于太阳黑子数据的中间周期(根据研究，有 11 年、22 年和 100 年三个周期)，窗口将被设置为22。代码很简单，在本书代码包的ch-07.ipynb文件中给出(这里的数据仅限于过去 150 年，以便在图中进行比较):

import matplotlib.pyplot as plt 
import statsmodels.api as sm 
from pandas.stats.moments import rolling_window 
import pandas as pd 

data_loader = sm.datasets.sunspots.load_pandas() 
df = data_loader.data.tail(150) 
df = pd.DataFrame({'SUNACTIVITY':df['SUNACTIVITY'].values}, index=df['YEAR']) 
ax = df.plot() 

def plot_window(wintype): 
    df2 = df.rolling(window=22,win_type=wintype, 
center=False,axis=0).mean() 
    df2.columns = [wintype] 
    df2.plot(ax=ax) 

plot_window('boxcar') 
plot_window('triang') 
plot_window('blackman') 
plot_window('hanning') 
plot_window('bartlett') 
plt.show() 

最终结果见下图:

定义协整

协整类似于相关性，但被许多人视为定义两个时间序列相关性的优越指标。如果两个时间序列x(t)和y(t)的线性组合是平稳的，则它们是协整的。在这种情况下，下面的等式应该是固定的:

y(t) - a x(t) 

想想一个醉汉和他的狗出去散步。相关性告诉我们它们是否在同一个方向。协整告诉我们一些关于人和他的狗之间的时间距离。我们将使用随机生成的时间序列和真实数据来展示协整性。增广 Dickey-Fuller ( ADF )测试(见http://en . Wikipedia . org/wiki/增广 _Dickey%E2%80%93Fuller_test )测试时间序列中的单位根，可用于确定时间序列的协整性。

对于下面的代码，请看一下本书代码包中的ch-07.ipynb文件:

import statsmodels.api as sm 
from pandas.stats.moments import rolling_window 
import pandas as pd 
import statsmodels.tsa.stattools as ts 
import numpy as np 

def calc_adf(x, y): 
    result = sm.OLS(x, y).fit()     
    return ts.adfuller(result.resid) 

data_loader = sm.datasets.sunspots.load_pandas() 
data = data_loader.data.values 
N = len(data) 

t = np.linspace(-2 * np.pi, 2 * np.pi, N) 
sine = np.sin(np.sin(t)) 
print("Self ADF", calc_adf(sine, sine)) 

noise = np.random.normal(0, .01, N) 
print("ADF sine with noise", calc_adf(sine, sine + noise)) 

cosine = 100 * np.cos(t) + 10 
print("ADF sine vs cosine with noise", calc_adf(sine, cosine + noise)) 

print("Sine vs sunspots", calc_adf(sine, data)) 

让我们从协整演示开始:

1. 定义以下函数来计算 ADF 统计量:

           def calc_adf(x, y): 
               result = sm.OLS(x, y).fit()     
               return ts.adfuller(result.resid) 
   
   

2. 将太阳黑子数据加载到 NumPy 数组中:

           data_loader = sm.datasets.sunspots.load_pandas() 
           data = data_loader.data.values 
           N = len(data) 
   
   

3. 生成sine并计算sine与自身的协整:

           t = np.linspace(-2 * np.pi, 2 * np.pi, N) 
           sine = np.sin(np.sin(t)) 
           print("Self ADF", calc_adf(sine, sine)) 
   
   

4. The code should print the following:

           Self ADF (-5.0383000037165746e-16, 0.95853208606005591, 0, 308,  
           {'5%': -2.8709700936076912, '1%': -3.4517611601803702, '10%':  
           -2.5717944160060719}, -21533.113655477719)
   
   

   打印输出中的第一个值是 ADF 指标，第二个值是 p 值。如你所见，p 值非常高。以下值是滞后和样本大小。最后的字典给出了这个精确样本量的 t 分布值。

5. 现在，向sine添加噪声，演示噪声将如何影响信号:

           noise = np.random.normal(0, .01, N) 
           print("ADF sine with noise", calc_adf(sine, sine + noise)) 
   
   

6. With the noise, we get the following results:

          ADF sine with noise (-7.4535502402193075, 5.5885761455106898e-
          11, 3, 305, {'5%': -2.8710633193086648, '1%': 
          -3.4519735736206991, '10%': -2.5718441306100512}, 
          -1855.0243977703672)
   
   

   p 值已经大幅下降。这里的 ADF 度量-7.45低于字典中所有的临界值。所有这些都是拒绝协整的有力论据。

7. 让我们生成一个更大幅度和偏移的cosine。再次，让我们添加噪音:

        cosine = 100 * np.cos(t) + 10 
        print("ADF sine vs cosine with noise", calc_adf(sine, cosine + 
        noise)) 

将打印以下值:

ADF sine vs cosine with noise (-17.927224617871534, 2.8918612252729532e-30, 16, 292, {'5%': -2.8714895534256861, '1%': -3.4529449243622383, '10%': -2.5720714378870331}, -11017.837238220782)

同样，我们有充分的理由拒绝协整。检查sine和sunspots之间的协整给出以下输出:

Sine vs sunspots (-6.7242691810701016, 3.4210811915549028e-09, 16, 292, {'5%': -2.8714895534256861, '1%': -3.4529449243622383, '10%': -2.5720714378870331}, -1102.5867415291168)

这里使用的对的置信水平大致相同，因为它们取决于数据点的数量，而数据点的数量变化不大。下表总结了结果:

| **配对** | **统计** | **p 值** | **5%** | **1%** | **10%** | **拒绝** | | 自正弦 | -5.03E-16 | Zero point nine five | -2.87 | -3.45 | -2.57 | 不 | | 正弦与正弦噪声 | -7.45 | 5.58E-11 | -2.87 | -3.45 | -2.57 | 是 | | 正弦和余弦与噪声的关系 | -17.92 | 2.89E-30 | -2.87 | -3.45 | -2.57 | 是 | | 正弦与太阳黑子 | -6.72 | 3.42E-09 | -2.87 | -3.45 | -2.57 | 是 |

自相关

自相关是数据集内的相关性，可以指示趋势。

对于给定的时间序列，在已知均值和标准差的情况下，我们可以使用期望值运算符定义时间s和t的自相关，如下所示:

本质上，这是应用于时间序列和滞后时间序列的相关公式。

例如，如果我们有一个周期的滞后，我们可以检查前一个值是否影响当前值。要做到这一点，自相关值必须相当高。

在前一章第 6 章、数据可视化中，我们已经使用了绘制自相关的 Pandas 函数。在这个例子中，我们将使用 NumPy correlate()函数来计算太阳黑子周期的实际自相关值。最后，我们需要将我们收到的值规范化。按如下方式应用数字按钮correlate()功能:

y = data - np.mean(data) 
norm = np.sum(y ** 2) 
correlated = np.correlate(y, y, mode='full')/norm 

我们还对与最高相关性相对应的指数感兴趣。这些索引可以通过 NumPy argsort()函数找到，该函数返回对数组进行排序的索引:

print np.argsort(res)[-5:] 

这些是最大自相关的指数:

[ 9 11 10  1  0]  

根据定义，最大的自相关是零滞后，即信号与其自身的相关性。第二大值是滞后 1 年和 10 年。查看本书代码包中的ch-07.ipynb文件:

import numpy as np 
import pandas as pd 
import statsmodels.api as sm 
import matplotlib.pyplot as plt 
from pandas.tools.plotting import autocorrelation_plot 

data_loader = sm.datasets.sunspots.load_pandas() 
data = data_loader.data["SUNACTIVITY"].values 
y = data - np.mean(data) 
norm = np.sum(y ** 2) 
correlated = np.correlate(y, y, mode='full')/norm 
res = correlated[len(correlated)/2:] 

print(np.argsort(res)[-5:]) 
plt.plot(res) 
plt.grid(True) 
plt.xlabel("Lag") 
plt.ylabel("Autocorrelation") 
plt.show() 
autocorrelation_plot(data) 
plt.show() 

最终结果见下图:

将之前的剧情与 Pandas 制作的剧情进行对比:

自回归模型

一个自回归模型可以用来表示一个时间序列，目标是预测未来值。在这样的模型中，假设一个变量依赖于它以前的值。该关系也被假设为线性的，并且我们需要拟合数据以便找到数据的参数。自回归模型的数学公式如下:

在前面的公式中，c是常数，最后一项是随机分量，也称为白噪声。

这给我们带来了非常常见的线性回归问题。出于实际原因，保持模型简单并且只包含必要的滞后组件是很重要的。在机器学习术语中，这些被称为特征。对于回归问题，Python 机器学习 scikit-learn 库即使不是最好的选择，也是一个不错的选择。我们将在第 10 章、预测性分析和机器学习中使用该应用编程接口。

在回归设置中，我们经常会遇到过度拟合的问题——当我们对一个样本有完美的拟合时，这个问题就会出现，当我们引入新的数据点时，这个样本的性能会很差。标准的解决方案是应用交叉验证(或者使用避免过度拟合的算法)。在这种方法中，我们估计部分样本的模型参数。其余的数据用于测试和评估模型。这其实是一个简化的解释。还有更复杂的交叉验证方案，其中很多都得到 scikit-learn 的支持。为了评估模型，我们可以计算适当的评估指标。可以想象，有许多度量标准，由于实践者的不断调整，这些度量标准可以有不同的定义。我们可以在书上或维基百科上找到这些定义。重要的是要记住，对预测或拟合的评估不是一门精确的科学。有这么多指标的事实只能证实这一点。

我们将使用前一节中找到的前两个滞后分量，使用scipy.optimize.leastsq()函数建立模型。我们可以选择线性代数函数来代替。然而，leastsq()功能更加灵活，让我们可以指定几乎任何类型的模型。按照以下步骤设置模型:

def model(p, x1, x10): 
   p1, p10 = p 
   return p1 * x1 + p10 * x10 

def error(p, data, x1, x10): 
   return data - model(p, x1, x10) 

要拟合模型，初始化参数列表并将其传递给leastsq()函数，如下所示:

def fit(data): 
   p0 = [.5, 0.5] 
   params = leastsq(error, p0, args=(data[10:], data[9:-1], data[:-10]))[0] 
   return params 

根据部分数据训练模型:

cutoff = .9 * len(sunspots) 
params = fit(sunspots[:cutoff]) 
print "Params", params 

以下是我们得到的参数:

Params [ 0.67172672  0.33626295]

有了这些参数，我们将绘制预测值并计算各种指标。以下是我们获得的指标值:

Root mean square error 22.8148122613
Mean absolute error 17.6515446503
Mean absolute percentage error 60.7817800736
Symmetric Mean absolute percentage error 34.9843386176
Coefficient of determination 0.799940292779

有关最终结果，请参考下图:

似乎我们有许多几乎是准确无误的预测，但也有一堆相当遥远的预测。总的来说，我们没有完美的契合；然而，这并不是一场彻底的灾难。它在中间的某个地方。

以下代码在本书代码包的ch-07.ipynb文件中:

from scipy.optimize import leastsq 
import statsmodels.api as sm 
import matplotlib.pyplot as plt 
import numpy as np 

def model(p, x1, x10): 
   p1, p10 = p 
   return p1 * x1 + p10 * x10 

def error(p, data, x1, x10): 
   return data - model(p, x1, x10) 

def fit(data): 
   p0 = [.5, 0.5] 
   params = leastsq(error, p0, args=(data[10:], data[9:-1], data[:-10]))[0] 
   return params 

data_loader = sm.datasets.sunspots.load_pandas() 
sunspots = data_loader.data["SUNACTIVITY"].values 

cutoff = .9 * len(sunspots) 
params = fit(sunspots[:cutoff]) 
print("Params", params) 

pred = params[0] * sunspots[cutoff-1:-1] + params[1] * sunspots[cutoff-10:-10] 
actual = sunspots[cutoff:] 
print("Root mean square error", np.sqrt(np.mean((actual - pred) ** 2))) 
print("Mean absolute error", np.mean(np.abs(actual - pred))) 
print("Mean absolute percentage error", 100 * np.mean(np.abs(actual - pred)/actual)) 
mid = (actual + pred)/2 
print("Symmetric Mean absolute percentage error", 100 * np.mean(np.abs(actual - pred)/mid)) 
print("Coefficient of determination", 1 - ((actual - pred) ** 2).sum()/ ((actual - actual.mean()) ** 2).sum()) 
year_range = data_loader.data["YEAR"].values[cutoff:] 
plt.plot(year_range, actual, 'o', label="Sunspots") 
plt.plot(year_range, pred, 'x', label="Prediction") 
plt.grid(True) 
plt.xlabel("YEAR") 
plt.ylabel("SUNACTIVITY") 
plt.legend() 
plt.show() 

ARMA 模型

ARMA 模型常用于预测一个时间序列。这些模型结合了自回归和移动平均模型。在移动平均模型中，我们假设一个变量是时间序列平均值和噪声成分线性组合的和。

注

自回归和移动平均模型可以有不同的顺序。一般来说，我们可以用p自回归项和q移动平均项定义一个 ARMA 模型，如下所示:

在前面的公式中，就像在自回归模型公式中一样，我们有一个常数和一个白噪声分量；然而，我们也试图拟合滞后噪声分量。

幸运的是，可以使用statsmodelssm.tsa.ARMA()例程进行此分析。将数据拟合到ARMA(10,1)模型，如下所示:

model = sm.tsa.ARMA(df, (10,1)).fit() 

执行预测(statsmodels 大量使用字符串):

prediction = model.predict('1975', str(years[-1]), dynamic=True) 

最终结果见下图:

拟合度很差，因为坦率地说，我们过度拟合了数据。上一节中更简单的模型效果更好。示例代码可以在本书代码包的ch-07.ipynb文件中找到:

import pandas as pd 
import matplotlib.pyplot as plt 
import statsmodels.api as sm 
import datetime 

data_loader = sm.datasets.sunspots.load_pandas() 
df = data_loader.data 
years = df["YEAR"].values.astype(int) 
df.index = pd.Index(sm.tsa.datetools.dates_from_range(str(years[0]), str(years[-1]))) 
del df["YEAR"] 

model = sm.tsa.ARMA(df, (10,1)).fit() 
prediction = model.predict('1975', str(years[-1]), dynamic=True) 

df['1975':].plot() 
prediction.plot(style='--', label='Prediction') 
plt.legend() 
plt.show() 

产生周期信号

许多自然现象就像一个精确的时钟一样有规律，值得信赖。有些现象呈现出似乎有规律的模式。一组科学家通过希尔伯特-黄变换发现了太阳黑子活动的三个周期。周期的持续时间大约为 11 年、22 年和 100 年。通常，我们会使用三角函数(如正弦函数)来模拟周期信号。你可能还记得高中时的一些三角学。这就是我们这个例子所需要的。因为我们有三个周期，所以创建一个模型似乎是合理的，它是三个正弦函数的线性组合。这仅仅需要对自回归模型的代码进行微小的调整。以下代码请参考本书代码包中的periodic.py文件:

from scipy.optimize import leastsq 
import statsmodels.api as sm 
import matplotlib.pyplot as plt 
import numpy as np 
def model(p, t): 
   C, p1, f1, phi1 , p2, f2, phi2, p3, f3, phi3 = p 
   return C + p1 * np.sin(f1 * t + phi1) + p2 * np.sin(f2 * t + phi2) +p3 * np.sin(f3 * t + phi3) 

def error(p, y, t): 
   return y - model(p, t) 

def fit(y, t): 
   p0 = [y.mean(), 0, 2 * np.pi/11, 0, 0, 2 * np.pi/22, 0, 0, 2 * np.pi/100, 0] 
   params = leastsq(error, p0, args=(y, t))[0] 
   return params 

data_loader = sm.datasets.sunspots.load_pandas() 
sunspots = data_loader.data["SUNACTIVITY"].values 
years = data_loader.data["YEAR"].values 

cutoff = .9 * len(sunspots) 
params = fit(sunspots[:cutoff], years[:cutoff]) 
print("Params", params) 

pred = model(params, years[cutoff:]) 
actual = sunspots[cutoff:] 
print("Root mean square error", np.sqrt(np.mean((actual - pred) ** 2))) 
print("Mean absolute error", np.mean(np.abs(actual - pred))) 
print("Mean absolute percentage error", 100 * np.mean(np.abs(actual - pred)/actual)) 
mid = (actual + pred)/2 
print("Symmetric Mean absolute percentage error", 100 * np.mean(np.abs(actual - pred)/mid)) 
print("Coefficient of determination", 1 - ((actual - pred) ** 2).sum()/ ((actual - actual.mean()) ** 2).sum()) 
year_range = data_loader.data["YEAR"].values[cutoff:] 
plt.plot(year_range, actual, 'o', label="Sunspots") 
plt.plot(year_range, pred, 'x', label="Prediction") 
plt.grid(True) 
plt.xlabel("YEAR") 
plt.ylabel("SUNACTIVITY") 
plt.legend() 
plt.show() 

我们得到以下输出:

Params [ 47.18800285  28.89947419   0.56827284   6.51168446   4.55214999
0.29372077 -14.30926648 -18.16524041   0.06574835  -4.37789602]
Root mean square error 59.5619175499
Mean absolute error 44.5814573306
Mean absolute percentage error 65.1639657495
Symmetric Mean absolute percentage error 78.4477263927
Coefficient of determination -0.363525210982

第一行显示了我们尝试的模型的系数。我们的平均绝对误差为 44，这意味着我们平均在任一方向上偏离了这个量。我们还希望确定系数尽可能接近 1，以获得良好的拟合。相反，我们得到一个负值，这是不可取的。有关最终结果，请参考下图:

傅里叶分析

傅里叶分析是基于以数学家约瑟夫·傅里叶命名的傅里叶级数而产生的。傅立叶级数是一种数学方法，用于将函数表示为正弦和余弦项的无穷级数。所讨论的函数可以是实数或复数:

傅里叶分析最有效的算法是快速傅里叶变换 ( 快速傅里叶变换)。该算法在 SciPy 和 NumPy 中实现。当应用于时间序列数据时，傅立叶分析将映射转换到频域，产生频谱。频谱将谐波显示为特定频率的明显尖峰。例如，音乐是由不同的频率组成的，音符 A 在 440 赫兹。音符 A 可以由音叉产生。我们可以用钢琴之类的乐器来演奏这个和其他的音符。白噪声是由许多频率组成的信号，这些频率被相等地表示。白光是光的所有可见频率的混合，也同样表现出来。

在下面的例子中，我们将导入两个函数(参考ch-07.ipynb):

from scipy.fftpack import rfft 
from scipy.fftpack import fftshift 

rfft()功能对实值数据进行快速傅立叶变换。我们也可以使用fft()函数，但是它在这个特定的数据集上给出了一个警告。fftshift()功能将零频率分量(数据的平均值)移动到频谱的中间，以便更好地可视化。我们还将看一下正弦波，因为这很容易理解。创建正弦波并对其应用快速傅立叶变换:

t = np.linspace(-2 * np.pi, 2 * np.pi, len(sunspots)) 
mid = np.ptp(sunspots)/2 
sine = mid + mid * np.sin(np.sin(t)) 

sine_fft = np.abs(fftshift(rfft(sine))) 
print "Index of max sine FFT", np.argsort(sine_fft)[-5:] 

以下输出显示了对应于最大振幅的指数:

Index of max sine FFT [160 157 166 158 154]

对太阳黑子数据进行快速傅立叶变换:

transformed = np.abs(fftshift(rfft(sunspots))) 
print "Indices of max sunspots FFT", np.argsort(transformed)[-5:] 

光谱中的五个最大峰值可以在以下指数中找到:

Indices of max sunspots FFT [205 212 215 209 154]

最大的山峰也位于 154。最终结果见下图:

完整的代码位于本书代码包的ch-07.ipynb文件中:

import numpy as np 
import statsmodels.api as sm 
import matplotlib.pyplot as plt 
from scipy.fftpack import rfft 
from scipy.fftpack import fftshift 

data_loader = sm.datasets.sunspots.load_pandas() 
sunspots = data_loader.data["SUNACTIVITY"].values 

t = np.linspace(-2 * np.pi, 2 * np.pi, len(sunspots)) 
mid = np.ptp(sunspots)/2 
sine = mid + mid * np.sin(np.sin(t)) 

sine_fft = np.abs(fftshift(rfft(sine))) 
print("Index of max sine FFT", np.argsort(sine_fft)[-5:]) 

transformed = np.abs(fftshift(rfft(sunspots))) 
print("Indices of max sunspots FFT", np.argsort(transformed)[-5:]) 

plt.subplot(311) 
plt.plot(sunspots, label="Sunspots") 
plt.plot(sine, lw=2, label="Sine") 
plt.grid(True) 
plt.legend() 
plt.subplot(312) 
plt.plot(transformed, label="Transformed Sunspots") 
plt.grid(True) 
plt.legend() 
plt.subplot(313) 
plt.plot(sine_fft, lw=2, label="Transformed Sine") 
plt.grid(True) 
plt.legend() 
plt.show() 

光谱分析

在前一节中，我们绘制了数据集的振幅谱。物理信号的功率谱可视化了信号的能量分布。我们可以很容易地修改代码来绘制功率谱，只需将值平方如下:

plt.plot(transformed ** 2, label="Power Spectrum") 

相位谱将相位(正弦函数的初始角度)可视化，可以绘制如下:

plt.plot(np.angle(transformed), label="Phase Spectrum") 

有关最终结果，请参考下图:

完整代码请参考本书代码包中的ch-07.ipynb文件。

过滤

滤波是一种信号处理，包括去除或抑制信号的一部分。在应用快速傅立叶变换后，我们可以过滤高频或低频，或者我们可以尝试去除白噪声。白噪声是具有恒定功率谱的随机信号，因此不包含任何有用的信息。scipy.signal包有许多用于过滤的实用程序。在本例中，我们将演示这些例程的一个小示例:

• 中值滤波器计算滚动窗口中的中值(参见en.wikipedia.org/wiki/Median…)。它由medfilt()函数实现，该函数有一个可选的窗口大小参数。
• 维纳滤波器使用统计数据去除噪声(参见en.wikipedia.org/wiki/Wiener…)。对于滤波器g(t)和信号s(t)，通过卷积(g * [s + n])(t)计算输出。通过wiener()功能实现。该函数还有一个可选的窗口大小参数。
• 去趋势过滤器去除趋势。这可以是线性或恒定趋势。通过detrend()功能实现。

以下代码请参考本书代码包中的ch-07.ipynb文件:

import statsmodels.api as sm 
import matplotlib.pyplot as plt 
from scipy.signal import medfilt 
from scipy.signal import wiener 
from scipy.signal import detrend 

data_loader = sm.datasets.sunspots.load_pandas() 
sunspots = data_loader.data["SUNACTIVITY"].values 
years = data_loader.data["YEAR"].values 

plt.plot(years, sunspots, label="SUNACTIVITY") 
plt.plot(years, medfilt(sunspots, 11), lw=2, label="Median") 
plt.plot(years, wiener(sunspots, 11), '--', lw=2, label="Wiener") 
plt.plot(years, detrend(sunspots), lw=3, label="Detrend") 
plt.xlabel("YEAR") 
plt.grid(True) 
plt.legend() 
plt.show() 

有关最终结果，请参考下图:

总结

在本章中，时间序列示例使用了年太阳黑子周期数据。

您了解到，在过去的固定时间段内，在同一时间序列中，试图推导一个值与另一个数据点或数据点组合之间的关系是很常见的。

移动平均值指定了以前看到的数据的窗口，每次窗口向前滑动一个周期时，移动平均值就会被平均。在 Pandas 应用编程接口中，DataFrame.rolling()函数为窗口函数功能提供了与不同窗口函数对应的不同的win_type字符串参数值。

协整类似于相关性，是定义两个时间序列相关性的度量。在回归设置中，我们经常遇到过度拟合的问题。当我们完全适合一个样本时，这个问题就出现了，当我们引入新的数据点时，这个样本表现不佳。为了评估模型，我们可以计算适当的评估指标。

数据库是数据分析的重要工具。关系数据库从 20 世纪 70 年代就已经出现了。最近，NoSQL 数据库已经成为一个可行的选择。下一章，第 8 章、使用数据库，包含关于各种数据库(关系数据库和 NoSQL 数据库)和相关应用编程接口的信息。