精通-NumPy-数值分析-全-
精通 NumPy 数值分析（全）
零、前言
如果您想磨练数据科学领域的技能，那么有很多书籍和课程的难度会有所不同。 通常会发生的情况是，您开始研究入门资源，然后继续使用更深入的技术资源，以尝试新的领域或技术。 如果您一段时间以来一直遵循这种学习方式，那么您一定已经意识到这将变得非常耗时。 作为终身学习者，我们需要书本能更紧凑地表示知识和经验，这需要在理论和实践之间取得适当的平衡。 本书旨在将初学者，中级和高级概念结合在一起，这是我们从头开始积累您的知识的谦虚努力。
本书假定以前没有科学计算的背景，并且将通过实际例子介绍各种主题。 有时可能感觉像是将各个主题随机组合在一起，而本书的流程并没有遵循一条一致的道路。 这是我们故意做出的决定，目的是让您对几种不同的主题和应用有所了解。
我们希望您能读这本书，以更广泛地了解科学计算，并掌握 NumPy 的精髓和其他支持 Python 的科学库，例如 SciPy 和 Scikit-Learn。
这本书是给谁的
本书适用于希望在数据科学领域获得更多知识的每个人。 如果您是 Python 程序员，数据分析师，数据工程师或数据科学爱好者，并且希望掌握 NumPy 的复杂性并为您的数字和科学构建解决方案，那么《使用 NumPy 掌握数字计算问题》非常适合您。 希望您熟悉数学以充分利用本书。
充分利用这本书
- 基本的 Python 编程知识肯定会有所帮助，尽管并非绝对必要
- Python 3 的 Anaconda 发行版足以覆盖本书中使用的大多数示例
使用约定
本书中使用了许多文本约定。
CodeInText：表示文本中的词，数据库表名称，文件夹名称，文件名，文件扩展名，路径名，伪 URL，用户输入和 Twitter 句柄。 这是一个示例：“此函数中的另一个重要参数是learning_rate。”
代码块设置如下：
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)'])
任何命令行输入或输出的编写方式如下：
$ sudo apt-get update
$ sudo apt-get upgrade
粗体：表示您在屏幕上看到的新术语，重要单词或顺序。 例如，菜单或对话框中的单词会出现在这样的文本中。 这是一个示例：“因变量是我们要预测的变量。”
警告或重要提示显示如下，提示和技巧显示如下。
一、使用 NumPy 数组
科学计算是一个多学科领域，其应用跨越数值分析，计算金融和生物信息学等学科。
让我们考虑一下金融市场的情况。 当您考虑金融市场时，会有巨大的相互联系的互动网络。 政府，银行，投资基金，保险公司，养老金，个人投资者和其他人都参与了这种金融工具的交换。 您不能简单地模拟市场参与者之间的所有互动，因为参与金融交易的每个人都有不同的动机和不同的风险/回报目标。 还有其他因素会影响金融资产的价格。 即使为一个资产价格建模也需要您做大量工作，并且不能保证成功。 用数学术语来说，这没有封闭形式的解决方案，这为利用科学计算提供了一个很好的案例，您可以在其中使用高级计算技术来解决此类问题。
通过编写计算机程序，您将有能力更好地了解正在使用的系统。 通常，您将要编写的计算机程序将是某种模拟，例如蒙特卡洛模拟。 通过使用诸如蒙特卡洛的模拟，您可以对期权合约的价格进行建模。 仅仅由于金融市场的复杂性，对金融资产进行定价是进行模拟的良好材料。 在进行计算时，所有这些数学计算都需要一个功能强大，可扩展且方便的数据结构（大多数为矩阵形式）。 换句话说，您需要比列表更紧凑的结构，以简化您的任务。 NumPy 是高效向量/矩阵运算的理想选择，其广泛的数学运算库使数值计算变得简单而高效。
在本章中，我们将涵盖以下主题：
- NumPy 的重要性
- 关于向量和矩阵的理论和实践信息
- NumPy 数组操作及其在多维数组中的用法
问题是，我们应该从哪里开始练习编码技能？ 在本书中，由于 Python 在科学界的广泛应用，您将在使用它，并且您将主要使用名为 NumPy 的特定库（它代表数字 Python）。
技术要求
在本书中，我们将使用 Jupyter 笔记本。 我们将通过网络浏览器编辑和运行 Python 代码。 这是一个开源平台，您可以按照此链接中的说明进行安装。
本书将使用 Python 3.x，因此在打开新笔记本时，应选择 Python 3 内核。 另外，也可以使用 Anaconda（Python 版本 3.6）安装 Jupyter 笔记本。 您可以按照此链接中的说明进行安装。
为什么我们需要 NumPy？
Python 最近成为一种摇滚明星的编程语言，不仅因为它具有友好的语法和可读性，而且因为它可以用于多种用途。 各种库的 Python 生态系统使程序员相对容易进行各种计算。 堆栈溢出是程序员最受欢迎的网站之一。 用户可以通过标记与之相关的编程语言来提问。 下图通过计算这些标签显示了主要编程语言的增长，并绘制了多年来主流编程语言的流行程度。 通过 Stack Overflow 进行的研究可以通过其官方博客链接进行进一步分析：
主要编程语言的发展
NumPy 是 Python 中科学计算的最基本软件包，也是许多其他软件包的基础。 由于 Python 最初并不是为数字计算而设计的，因此在 90 年代后期开始出现这种需求时，Python 开始在需要更快的向量运算的工程师和程序员中流行。 从下图可以看到，许多流行的机器学习和计算包都使用了 NumPy 的某些功能，最重要的是它们在其方法中大量使用了 NumPy 数组，这使 NumPy 成为科学项目的基本库。
该图显示了一些使用 NumPy 功能的知名库：
NumPy 栈
对于数值计算，您主要处理向量和矩阵。 您可以通过使用一系列数学函数以不同的方式来操作它们。 NumPy 非常适合此类情况，因为它允许用户有效地完成其计算。 尽管 Python 列表很容易创建和操作，但它们不支持向量化操作。 Python 在列表中没有固定的类型元素，例如，for循环效率不高，因为在每次迭代中都需要检查数据类型。 但是，在 NumPy 数组中，数据类型是固定的，并且还支持向量化运算。 与 Python 列表相比，NumPy 在多维数组操作中不仅效率更高。 它还提供了许多数学方法，您可以在导入后立即应用它们。 NumPy 是用于科学 Python 数据科学堆栈的核心库。
SciPy 与 NumPy 有着密切的关系，因为它使用 NumPy 多维数组作为其线性代数，优化，插值，积分，FFT，信号，图像处理和其它的科学函数的基础数据结构。 SciPy 建立在 NumPy 数组框架之上，并凭借其先进的数学函数提升了科学编程能力。 因此，NumPy API 的某些部分已移至 SciPy。 在许多情况下，与 NumPy 的这种关系使 SciPy 更便于进行高级科学计算。
综上所述，我们可以总结出 NumPy 的优势如下：
- 它是开源的，零成本
- 这是一种具有用户友好语法的高级编程语言
- 比 Python 列表更有效
- 它具有更高级的内置函数，并与其他库很好地集成在一起
谁使用 NumPy？
在学术界和商业界，您都会听到人们谈论他们在工作中使用的工具和技术。 根据环境和条件，您可能需要使用特定技术。 例如，如果您的公司已经投资了 SAS，则需要在适合您问题的 SAS 开发环境中进行项目。
但是，NumPy 的优点之一是它是开源的，在您的项目中使用它无需花费任何成本。 如果您已经使用 Python 编写过代码，那么它非常容易学习。 如果您关心性能，则可以轻松嵌入 C 或 Fortran 代码。 此外，它将为您介绍其他完整的库集，例如 SciPy 和 Scikit-learn，您可以使用它们来解决几乎所有问题。
由于数据挖掘和预测分析在最近变得非常重要，因此数据科学家和数据分析师等角色在 21 世纪最热门的工作中被提及，例如《福布斯》，彭博社， 等等。 需要处理数据并进行分析，建模或预测的人员应该熟悉 NumPy 的用法及其功能，因为它将帮助您快速创建原型并测试您的想法。 如果您是专业工作人员，那么您的公司很可能希望使用数据分析方法，以使其领先于竞争对手。 如果他们能够更好地理解他们拥有的数据，那么他们就可以更好地理解业务，这将使他们做出更好的决策。 NumPy 在这里起着至关重要的作用，因为它能够执行各种操作，并使您的项目及时有效。
向量和矩阵简介
矩阵是一组数字或元素，它们以矩形数组的形式排列。 矩阵的行和列通常按字母索引。 对于n x m矩阵，n表示行数，m表示列数。 如果我们有一个假设的n×m矩阵，则其结构如下：
如果n = m，则称为称为方阵：
向量实际上是具有多于一个元素的一行或一列的矩阵。 也可以将其定义为1 x m或n x 1矩阵。 您可以将向量解释为m维空间中的箭头或方向。 通常，大写字母表示矩阵，例如X；小写字母带有下标，例如x[11]，表示矩阵X的元素。
此外，还有一些重要的特殊矩阵：零矩阵（空矩阵）和恒等矩阵。0表示零矩阵，它是所有0的矩阵（MacDufee 1943 p.27）。 在0矩阵中，添加下标是可选的：
由I表示的单位矩阵及其对角线元素为1，其他元素为0：
将矩阵X与单位矩阵相乘时，结果将等于X：
单位矩阵对于计算矩阵的逆非常有用。 当您将任何给定矩阵与其逆矩阵相乘时，结果将是一个单位矩阵：
让我们简要地看一下 NumPy 数组上的矩阵代数。 矩阵的加减法运算与具有普通单数的数学方程式相似。 举个例子：
标量乘法也非常简单。 例如，如果将矩阵X乘以4，则唯一要做的就是将每个元素乘以4值，如下所示：
开始时矩阵处理看似复杂的部分是矩阵乘法。
假设您有两个矩阵，分别为X和Y，其中X是 
矩阵，Y是 
矩阵：
这两个矩阵的乘积将如下所示：
因此，乘积矩阵的每个元素的计算方式如下：
如果您不了解该符号，请不要担心。 以下示例将使事情变得更清楚。 您有矩阵X和Y，目标是获得这些矩阵的矩阵乘积：
基本思想是X的i行与Y的j列的乘积将成为结果矩阵的第i, j个元素。 乘法将从X的第一行和Y的第一列开始，因此它们的乘积将为Z[1, 1]：
您可以使用以下四行代码轻松地交叉检查结果：
In [1]: import numpy as np
x = np.array([[1,0,4],[3,3,1]])
y = np.array([[2,5],[1,1],[3,2]])
x.dot(y)
Out[1]: array([[14, 13],[12, 20]])
先前的代码块只是演示如何使用 NumPy 计算两个矩阵的点积。 在后面的章节中，我们将更深入地研究矩阵运算和线性代数。
NumPy 数组对象的基础
如上一节所述，使 NumPy 与众不同的是使用称为ndarrays的多维数组。 所有ndarray项目都是同类的，并且在内存中使用相同的大小。 让我们首先导入 NumPy，然后通过创建数组来分析 NumPy 数组对象的结构。 您可以通过在控制台中键入以下语句来轻松导入该库。 您可以使用任何命名约定代替np，但是在本书中，将使用np作为标准约定。 让我们创建一个简单的数组，并说明 Python 在幕后所具有的属性作为所创建数组的元数据，即所谓的属性：
In [2]: import numpy as np
x = np.array([[1,2,3],[4,5,6]])
x
Out[2]: array([[1, 2, 3],[4, 5, 6]])
In [3]: print("We just create a ", type(x))
Out[3]: We just create a <class 'numpy.ndarray'>
In [4]: print("Our template has shape as" ,x.shape)
Out[4]: Our template has shape as (2, 3)
In [5]: print("Total size is",x.size)
Out[5]: Total size is 6
In [6]: print("The dimension of our array is " ,x.ndim)
Out[6]: The dimension of our array is 2
In [7]: print("Data type of elements are",x.dtype)
Out[7]: Data type of elements are int32
In [8]: print("It consumes",x.nbytes,"bytes")
Out[8]: It consumes 24 bytes
如您所见，我们对象的类型是 NumPy 数组。x.shape返回一个元组，该元组为您提供数组的维度作为输出，例如(n, m)。 您可以使用x.size获得数组中元素的总数。 在我们的示例中，我们总共有六个元素。 知道形状和大小等属性非常重要。 您了解的越多，您对计算的适应就越多。 如果您不知道数组的大小和大小，那么开始使用数组进行计算就不明智了。 在 NumPy 中，您可以使用x.ndim来检查数组的维数。 还有其他属性，例如dtype和nbytes，这些属性在检查内存使用情况并验证应该在数组中使用哪种数据类型时非常有用。 在我们的示例中，每个元素的数据类型为int32，总共消耗 24 个字节。 您还可以在创建数组时强制使用其中某些属性，例如dtype。 以前，数据类型是整数。 让我们将其切换为float，complex或uint（无符号整数）。 为了查看数据类型更改的作用，让我们分析一下字节消耗量，如下所示：
In [9]: x = np.array([[1,2,3],[4,5,6]], dtype = np.float)
print(x)
Out[9]: print(x.nbytes)
[[ 1\. 2\. 3.]
[ 4\. 5\. 6.]]
48
In [10]: x = np.array([[1,2,3],[4,5,6]], dtype = np.complex)
print(x)
print(x.nbytes)
Out[10]: [[ 1.+0.j 2.+0.j 3.+0.j]
[ 4.+0.j 5.+0.j 6.+0.j]]
96
In [11]: x = np.array([[1,2,3],[4,-5,6]], dtype = np.uint32)
print(x)
print(x.nbytes)
Out[11]: [[ 1 2 3]
[ 4 4294967291 6]]
24
如您所见，每种类型消耗不同数量的字节。 假设您有一个如下矩阵，并且您正在使用int64或int32作为数据类型：
In [12]: x = np.array([[1,2,3],[4,5,6]], dtype = np.int64)
print("int64 consumes",x.nbytes, "bytes")
x = np.array([[1,2,3],[4,5,6]], dtype = np.int32)
print("int32 consumes",x.nbytes, "bytes")
Out[12]: int64 consumes 48 bytes
int32 consumes 24 bytes
如果使用int64，则内存需求将增加一倍。 向自己问这个问题； 哪种数据类型就足够了？ 在您的数字大于 2,147,483,648 且小于 -2,147,483,647 之前，使用int32就足够了。 假设您有一个大小超过 100MB 的巨大数组。 在这种情况下，这种转换对性能起着至关重要的作用。
正如您在上一个示例中可能已经注意到的那样，当您更改数据类型时，每次都在创建一个数组。 从技术上讲，创建数组后无法更改dtype。 但是，您可以做的是再次创建它，或者使用新的dtype和astype属性复制现有的文件。 让我们使用新的dtype创建数组的副本。 下面是一个示例，说明如何使用astype属性也可以更改dtype：
In [13]: x_copy = np.array(x, dtype = np.float)
x_copy
Out[13]: array([[ 1., 2., 3.],
[ 4., 5., 6.]])
In [14]: x_copy_int = x_copy.astype(np.int)
x_copy_int
Out[14]: array([[1, 2, 3],
[4, 5, 6]])
请记住，使用astype属性时，即使将其应用于x_copy，它也不会更改x_copy的dtype。它保留x_copy，但创建一个x_copy_int：
In [15]: x_copy
Out[15]: array([[ 1., 2., 3.],
[ 4., 5., 6.]])
假设有一个案例，您正在一个研究小组中工作，该研究小组试图确定和计算每个患癌症患者的风险。 您有 100,000 条记录（行），其中每一行代表一位患者，每位患者具有 100 个特征（某些测试的结果）。 结果，您有(100000, 100)的数组：
In [16]: Data_Cancer= np.random.rand(100000,100)
print(type(Data_Cancer))
print(Data_Cancer.dtype)
print(Data_Cancer.nbytes)
Data_Cancer_New = np.array(Data_Cancer, dtype = np.float32)
print(Data_Cancer_New.nbytes)
Out[16]: <class 'numpy.ndarray'>
float64
80000000
40000000
从前面的代码中可以看到，仅通过更改dtype即可将它们的大小从 80MB 减小到 40MB。 我们得到的回报是小数点后的精度较低。 除了精确到 16 位小数点外，您将只有 7 位小数。 在某些机器学习算法中，精度可以忽略不计。 在这种情况下，您应该随意调整dtype，以最大程度地减少内存使用量。
NumPy 数组操作
本节将指导您使用 NumPy 创建和处理数字数据。 让我们从列表中创建一个 NumPy 数组开始：
In [17]: my_list = [2, 14, 6, 8]
my_array = np.asarray(my_list)
type(my_array)
Out[17]: numpy.ndarray
让我们用标量值做一些加法，减法，乘法和除法：
In [18]: my_array + 2
Out[18]: array([ 4, 16, 8, 10])
In [19]: my_array - 1
Out[19]: array([ 1, 13, 5, 7])
In [20]: my_array * 2
Out[20]: array([ 4, 28, 12, 16, 8])
In [21]: my_array / 2
Out[21]: array([ 1\. , 7\. , 3\. , 4\. ])
在列表中执行相同的操作要困难得多，因为该列表不支持向量化操作，并且您需要迭代其元素。 创建 NumPy 数组的方法有很多，现在您将使用这些方法之一来创建一个充满零的数组。 稍后，您将执行一些算术运算，以查看 NumPy 在两个数组之间的按元素运算中的行为：
In [22]: second_array = np.zeros(4) + 3
second_array
Out[22]: array([ 3., 3., 3., 3.])
In [23]: my_array - second_array
Out[23]: array([ -1., 11., 3., 5.])
In [24]: second_array / my_array
Out[24]: array([ 1.5 , 0.21428571, 0.5 , 0.375 ])
就像我们在前面的代码中所做的那样，您可以创建一个充满np.ones的数组或一个充满np.identity的标识数组，并执行与之前相同的代数运算：
In [25]: second_array = np.ones(4) + 3
second_array
Out[25]: array([ 4., 4., 4., 4.])
In [26]: my_array - second_array
Out[26]: array([ -2., 10., 2., 4.])
In [27]: second_array / my_array
Out[27]: array([ 2\. , 0.28571429, 0.66666667, 0.5 ])
它可以通过np.ones方法按预期工作，但是当您使用单位矩阵时，计算将返回(4, 4)数组，如下所示：
In [28]: second_array = np.identity(4)
second_array
Out[28]: array([[ 1., 0., 0., 0.],
[ 0., 1., 0., 0.],
[ 0., 0., 1., 0.],
[ 0., 0., 0., 1.]])
In [29]: second_array = np.identity(4) + 3
second_array
Out[29]: array([[ 4., 3., 3., 3.],
[ 3., 4., 3., 3.],
[ 3., 3., 4., 3.],
[ 3., 3., 3., 4.]])
In [30]: my_array - second_array
Out[30]: array([[ -2., 11., 3., 5.],
[ -1., 10., 3., 5.],
[ -1., 11., 2., 5.],
[ -1., 11., 3., 4.]])
这是从second_array的第一列的所有元素和第二列的second_element的所有元素中减去my_array的第一个元素，依此类推。 同样的规则也适用于除法。 请记住，即使形状不完全相同，也可以成功执行数组操作。 在本章的后面，您将了解当两个数组由于形状不同而无法进行计算时的广播错误：
In [31]: second_array / my_array
Out[31]: array([[ 2\. , 0.21428571, 0.5 , 0.375 ],
[ 1.5 , 0.28571429, 0.5 , 0.375 ],
[ 1.5 , 0.21428571, 0.66666667, 0.375 ],
[ 1.5 , 0.21428571, 0.5 , 0.5 ]])
创建 NumPy 数组最有用的方法之一是arange。 这将在起始值和结束值之间的给定间隔内返回一个数组。 第一个参数是数组的起始值，第二个参数是终止值（在该位置它不再创建值），第三个参数是间隔。 您可以选择将dtype定义为第四个参数。 默认间隔值为 1：
In [32]: x = np.arange(3,7,0.5)
x
Out[32]: array([ 3\. , 3.5, 4\. , 4.5, 5\. , 5.5, 6\. , 6.5])
当您无法确定间隔应该是多少时，还有另一种方法可以在起点和终点之间创建具有固定间隔的数组，但是您应该知道数组应该有多少个拆分：
In [33]: x = np.linspace(1.2, 40.5, num=20)
x
Out[33]: array([ 1.2 , 3.26842105, 5.33684211, 7.40526316, 9.47368421,
11.54210526, 13.61052632, 15.67894737, 17.74736842, 19.81578947,
21.88421053, 23.95263158, 26.02105263, 28.08947368, 30.15789474,
32.22631579, 34.29473684, 36.36315789, 38.43157895, 40.5 ])
有两种用法相似的不同方法，但由于它们的基本尺度不同，它们返回不同的数字序列。 这意味着数字的分布也将不同。 第一个是geomspace，它以几何级数返回对数刻度上的数字：
In [34]: np.geomspace(1, 625, num=5)
Out[34]: array([ 1., 5., 25., 125., 625.])
另一个重要的方法是logspace，您可以在其中返回起点和终点的值，这些值在以下位置均匀缩放：
In [35]: np.logspace(3, 4, num=5)
Out[35]: array([ 1000\. , 1778.27941004, 3162.27766017, 5623.4132519 ,
10000\. ])
这些参数是什么？ 如果起始点为3，结束点为4，则这些函数将返回比初始范围大得多的数字。 实际上，您的起点默认设置为10，终点也设置为10。 因此，从技术上讲，在此示例中，起点为10**3，终点为10**4。 您可以避免这种情况，并使起点和终点与在方法中将其用作参数时相同。 诀窍是使用给定参数的以 10 为底的对数：
In [36]: np.logspace(np.log10(3) , np.log10(4) , num=5)
Out[36]: array([ 3\. , 3.2237098 , 3.46410162, 3.72241944, 4\. ])
到目前为止，您应该熟悉创建具有不同分布的数组的不同方法。 您还学习了如何对它们进行一些基本操作。 让我们继续您在日常工作中肯定会使用的其他有用函数。 大多数时候，您将不得不使用多个数组，并且需要非常快速地比较它们。 NumPy 对这个问题有很好的解决方案。 您可以像比较两个整数一样比较数组：
In [37]: x = np.array([1,2,3,4])
y = np.array([1,3,4,4])
x == y
Out[37]: array([ True, False, False, True], dtype=bool)
比较是逐个元素进行的，无论元素是否在两个不同的数组中匹配，它都会返回一个布尔向量。 此方法在小大小数组中效果很好，以及也为您提供了更多详细信息。 您可以从输出数组中看到以False表示的值，这些索引值在这两个数组中不匹配。 如果数组很大，则还可以选择对问题的单个答案，即元素是否在两个不同的数组中匹配：
In [38]: x = np.array([1,2,3,4])
y = np.array([1,3,4,4])
np.array_equal(x,y)
Out[38]: False
在这里，您只有一个布尔输出。 您只知道数组不相等，但是您不知道哪些元素完全不相等。 比较不仅限于检查两个数组是否相等。 您还可以在两个数组之间进行逐元素的上下比较：
In [39]: x = np.array([1,2,3,4])
y = np.array([1,3,4,4])
x < y
Out[39]: array([False, True, True, False], dtype=bool)
当需要进行逻辑比较（AND，OR和XOR）时，可以在数组中使用它们，如下所示：
In [40]: x = np.array([0, 1, 0, 0], dtype=bool)
y = np.array([1, 1, 0, 1], dtype=bool)
np.logical_or(x,y)
Out[40]: array([ True, True, False, True], dtype=bool)
In [41]: np.logical_and(x,y)
Out[41]: array([False, True, False, False], dtype=bool)
In [42]: x = np.array([12,16,57,11])
np.logical_or(x < 13, x > 50)
Out[42]: array([ True, False, True, True], dtype=bool)
到目前为止，已经介绍了诸如加法和乘法之类的代数运算。 我们如何将这些运算与指数函数，对数函数或三角函数之类的超越函数一起使用？
In [43]: x = np.array([1, 2, 3,4 ])
np.exp(x)
Out[43]: array([ 2.71828183, 7.3890561 , 20.08553692, 54.59815003])
In [44]: np.log(x)
Out[44]: array([ 0\. , 0.69314718, 1.09861229, 1.38629436])
In [45]: np.sin(x)
Out[45]: array([ 0.84147098, 0.90929743, 0.14112001, -0.7568025 ])
矩阵的转置呢？ 首先，您将reshape函数与arange一起使用以设置所需的矩阵形状：
In [46]: x = np.arange(9)
x
Out[46]: array([0, 1, 2, 3, 4, 5, 6, 7, 8])
In [47]: x = np.arange(9).reshape((3, 3))
x
Out[47]: array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [48]: x.T
Out[48]: array([[0, 3, 6],
[1, 4, 7],
[2, 5, 8]])
您转置了3 * 3数组，因此形状不会改变，因为两个大小均为3。 让我们看看没有平方数组时会发生什么：
In [49]: x = np.arange(6).reshape(2,3)
x
Out[49]: array([[0, 1, 2],
[3, 4, 5]])
In [50]: x.T
Out[50]: array([[0, 3],
[1, 4],
[2, 5]])
转置将按预期工作，并且大小也会切换。 您还可以从数组中获取摘要统计信息，例如均值，中位数和标准差。 让我们从 NumPy 提供的用于计算基本统计信息的方法开始：
| 方法 | 说明 |
|---|---|
np.sum |
返回数组所有值沿指定轴的总和 |
np.amin |
返回数组所有值或沿指定轴的最小值 |
np.amax |
返回数组所有值或沿指定轴的最大值 |
np.percentile |
返回数组所有值或沿指定轴的第 n 个的给定百分数 |
np.nanmin |
与np.amin相同，但忽略数组中的 NaN 值 |
np.nanmax |
与np.amax相同，但忽略数组中的 NaN 值 |
np.nanpercentile |
与np.percentile相同，但忽略数组中的 NaN 值 |
以下代码块给出了前面的 NumPy 统计方法的示例。 这些方法非常有用，因为您可以根据需要在整个数组中或轴向操作这些方法。 您应该注意，由于 SciPy 使用 NumPy 多维数组作为数据结构，因此可以在 SciPy 中找到这些方法的功能更完善的更好的实现：
In [51]: x = np.arange(9).reshape((3,3))
x
Out[51]: array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [52]: np.sum(x)
Out[52]: 36
In [53]: np.amin(x)
Out[53]: 0
In [54]: np.amax(x)
Out[54]: 8
In [55]: np.amin(x, axis=0)
Out[55]: array([0, 1, 2])
In [56]: np.amin(x, axis=1)
Out[56]: array([0, 3, 6])
In [57]: np.percentile(x, 80)
Out[57]: 6.4000000000000004
axis参数确定此函数将要作用的大小。 在此示例中，轴= 0代表第一个轴，即x轴，轴= 1代表第二个轴，即y。 当我们使用常规amin(x)时，我们返回单个值，因为它会计算所有数组中的最小值，但是当我们指定轴时，它将开始沿轴方向计算函数并返回一个数组，该数组显示每行或列的结果。 想象一下，您有很多人。 您可以使用amax找到最大值，但是如果您需要将此值的索引传递给另一个函数，将会发生什么？ 在这种情况下，argmin和argmax可以提供帮助，如以下代码片段所示：
In [58]: x = np.array([1,-21,3,-3])
np.argmax(x)
Out[58]: 2
In [59]: np.argmin(x)
Out[59]: 1
让我们继续更多的统计函数：
| 方法 | 说明 |
|---|---|
np.mean |
返回数组所有值或沿特定轴的平均值 |
np.median |
返回数组所有值或沿特定轴的中值 |
np.std |
返回数组所有值或沿特定轴的标准差 |
np.nanmean |
与np.mean相同，但忽略数组中的 NaN 值 |
np.nanmedian |
与np.nanmedian相同，但忽略数组中的 NaN 值 |
np.nonstd |
与np.nanstd相同，但忽略数组中的 NaN 值 |
以下代码提供了 NumPy 的先前统计方法的更多示例。 这些方法在数据发现阶段中大量使用，您可以在其中分析数据特征和分布：
In [60]: x = np.array([[2, 3, 5], [20, 12, 4]])
x
Out[60]: array([[ 2, 3, 5],
[20, 12, 4]])
In [61]: np.mean(x)
Out[61]: 7.666666666666667
In [62]: np.mean(x, axis=0)
Out[62]: array([ 11\. , 7.5, 4.5])
In [63]: np.mean(x, axis=1)
Out[63]: array([ 3.33333333, 12\. ])
In [64]: np.median(x)
Out[64]: 4.5
In [65]: np.std(x)
Out[65]: 6.3944420310836261
处理多维数组
本节将通过执行不同的矩阵运算使您对多维数组有一个简要的了解。
为了在 NumPy 中进行矩阵乘法，您必须使用dot()而不是*。 让我们看一些例子：
In [66]: c = np.ones((4, 4))
c*c
Out[66]: array([[ 1., 1., 1., 1.],
[ 1., 1., 1., 1.],
[ 1., 1., 1., 1.],
[ 1., 1., 1., 1.]])
In [67]: c.dot(c)
Out[67]: array([[ 4., 4., 4., 4.],
[ 4., 4., 4., 4.],
[ 4., 4., 4., 4.],
[ 4., 4., 4., 4.]])
使用多维数组时，最重要的主题是堆栈，即如何合并两个数组。hstack用于水平（列方式）堆叠数组，vstack用于垂直（行方式）堆叠数组。 您还可以使用hsplit和vsplit方法拆分列，方法与堆叠它们相同：
In [68]: y = np.arange(15).reshape(3,5)
x = np.arange(10).reshape(2,5)
new_array = np.vstack((y,x))
new_array
Out[68]: array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14],
[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9]])
In [69]: y = np.arange(15).reshape(5,3)
x = np.arange(10).reshape(5,2)
new_array = np.hstack((y,x))
new_array
Out[69]: array([[ 0, 1, 2, 0, 1],
[ 3, 4, 5, 2, 3],
[ 6, 7, 8, 4, 5],
[ 9, 10, 11, 6, 7],
[12, 13, 14, 8, 9]])
这些方法在机器学习应用中非常有用，尤其是在创建数据集时。 堆叠数组后，可以使用Scipy.stats检查其描述性统计信息。 假设有一个100记录，并且每个记录都具有10个特征，这意味着您有一个 2D 矩阵，其中包含100行和10列。 下面的示例说明如何轻松获取每个特征的一些描述性统计信息：
In [70]: from scipy import stats
x= np.random.rand(100,10)
n, min_max, mean, var, skew, kurt = stats.describe(x)
new_array = np.vstack((mean,var,skew,kurt,min_max[0],min_max[1]))
new_array.T
Out[70]: array([[ 5.46011575e-01, 8.30007104e-02, -9.72899085e-02,
-1.17492785e+00, 4.07031246e-04, 9.85652100e-01],
[ 4.79292653e-01, 8.13883169e-02, 1.00411352e-01,
-1.15988275e+00, 1.27241020e-02, 9.85985488e-01],
[ 4.81319367e-01, 8.34107619e-02, 5.55926602e-02,
-1.20006450e+00, 7.49534810e-03, 9.86671083e-01],
[ 5.26977277e-01, 9.33829059e-02, -1.12640661e-01,
-1.19955646e+00, 5.74237697e-03, 9.94980830e-01],
[ 5.42622228e-01, 8.92615897e-02, -1.79102183e-01,
-1.13744108e+00, 2.27821933e-03, 9.93861532e-01],
[ 4.84397369e-01, 9.18274523e-02, 2.33663872e-01,
-1.36827574e+00, 1.18986562e-02, 9.96563489e-01],
[ 4.41436165e-01, 9.54357485e-02, 3.48194314e-01,
-1.15588500e+00, 1.77608372e-03, 9.93865324e-01],
[ 5.34834409e-01, 7.61735119e-02, -2.10467450e-01,
-1.01442389e+00, 2.44706226e-02, 9.97784091e-01],
[ 4.90262346e-01, 9.28757119e-02, 1.02682367e-01,
-1.28987137e+00, 2.97705706e-03, 9.98205307e-01],
[ 4.42767478e-01, 7.32159267e-02, 1.74375646e-01,
-9.58660574e-01, 5.52410464e-04, 9.95383732e-01]])
NumPy 有一个很棒的模块numpy.ma，用于屏蔽数组元素。 当您要在计算时屏蔽（忽略）某些元素时，它非常有用。 当 NumPy 数组被屏蔽时，它将被视为无效，并且不考虑计算：
In [71]: import numpy.ma as ma
x = np.arange(6)
print(x.mean())
masked_array = ma.masked_array(x, mask=[1,0,0,0,0,0])
masked_array.mean()
2.5
Out[71]: 3.0
在前面的代码中，您有一个数组x = [0,1,2,3,4,5]。 您要做的是屏蔽数组的第一个元素，然后计算平均值。 当一个元素被屏蔽为1(True)时，数组中的关联索引值将被屏蔽。 在替换 NAN 值时，此方法也非常有用：
In [72]: x = np.arange(25, dtype = float).reshape(5,5)
x[x<5] = np.nan
x
Out[72]: array([[ nan, nan, nan, nan, nan],
[ 5., 6., 7., 8., 9.],
[ 10., 11., 12., 13., 14.],
[ 15., 16., 17., 18., 19.],
[ 20., 21., 22., 23., 24.]])
In [73]: np.where(np.isnan(x), ma.array(x, mask=np.isnan(x)).mean(axis=0), x)
Out[73]: array([[ 12.5, 13.5, 14.5, 15.5, 16.5],
[ 5\. , 6\. , 7\. , 8\. , 9\. ],
[ 10\. , 11\. , 12\. , 13\. , 14\. ],
[ 15\. , 16\. , 17\. , 18\. , 19\. ],
[ 20\. , 21\. , 22\. , 23\. , 24\. ]])
在前面的代码中，我们通过放置一个带有索引的条件将前五个元素的值更改为nan。x[x<5]指的是索引为 0、1、2、3 和 4 的元素。然后，我们用每列的平均值覆盖这些值（nan值除外）。 为了帮助您的代码更简洁，数组操作中还有许多其他有用的方法：
| 方法 | 说明 |
|---|---|
np.concatenate |
连接列表中的给定数组 |
np.repeat |
沿特定轴重复数组的元素 |
np.delete |
返回一个带有删除的子数组的新数组 |
np.insert |
在指定轴之前插入值 |
np.unique |
在数组中查找唯一值 |
np.tile |
通过以给定的重复次数重复给定的输入来创建数组 |
索引，切片，重塑，调整大小和广播
当您在机器学习项目中使用大量数组时，通常需要索引，切片，调整形状和调整大小。
索引是数学和计算机科学中使用的一个基本术语。 一般而言，索引可帮助您指定如何返回各种数据结构的所需元素。 下面的示例显示列表和元组的索引：
In [74]: x = ["USA","France", "Germany","England"]
x[2]
Out[74]: 'Germany'
In [75]: x = ('USA',3,"France",4)
x[2]
Out[75]: 'France'
在 NumPy 中，索引的主要用途是控制和操纵数组的元素。 这是一种创建通用查找值的方法。 索引包含三个子操作，分别是字段访问，基本切片和高级索引。 在字段访问中，只需指定数组中元素的索引即可返回给定索引的值。
NumPy 在索引和切片方面非常强大。 在许多情况下，您需要在数组中引用所需的元素，然后对该切片区域执行操作。 您可以对数组进行索引，类似于对元组或带有方括号的列表进行索引。 让我们从字段访问和使用一维数组的简单切片开始，然后继续使用更高级的技术：
In [76]: x = np.arange(10)
x
Out[76]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [77]: x[5]
Out[77]: 5
In [78]: x[-2]
Out[78]: 8
In [79]: x[2:8]
Out[79]: array([2, 3, 4, 5, 6, 7])
In [80]: x[:]
Out[80]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [81]: x[2:8:2]
Out[81]: array([2, 4, 6])
索引从0开始，因此在创建带有元素的数组时，第一个元素被索引为x[0]，与最后一个元素x[n-1]相同。 如您在前面的示例中看到的，x[5]指的是第六个元素。 您也可以在索引中使用负值。 NumPy 将这些值理解为倒数第n个。 像示例中一样，x[-2]指倒数第二个元素。 您还可以通过声明开始索引和结束索引来选择数组中的多个元素，也可以通过将增量级别作为第三个参数声明来创建顺序索引，如代码的最后一行所示。
到目前为止，我们已经看到了一维数组中的索引和切片。 逻辑不会改变，但是为了演示起见，我们也对多维数组进行一些练习。 当您拥有多维数组时，唯一发生变化的就是只有更多的轴。 您可以在以下代码中将 n 维数组切片为[x 轴切片，y 轴切片]：
In [82]: x = np.reshape(np.arange(16),(4,4))
x
Out[82]: array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
In [83]: x[1:3]
Out[83]: array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
In [84]: x[:,1:3]
Out[84]: array([[ 1, 2],
[ 5, 6],
[ 9, 10],
[13, 14]])
In [85]: x[1:3,1:3]
Out[85]: array([[ 5, 6],
[ 9, 10]])
您按行和列对数组进行了切片，但没有以更加不规则或更具动态性的方式对元素进行切片，这意味着您始终以矩形或正方形的方式对其进行切片。 想象一下我们要切片的4 * 4数组，如下所示：
为了获得前面的切片，我们执行以下代码：
In [86]: x = np.reshape(np.arange(16),(4,4))
x
Out[86]: array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
In [87]: x[[0,1,2],[0,1,3]]
Out[87]: array([ 0, 5, 11])
在高级索引中，第一部分指示要切片的行的索引，第二部分指示相应的列。 在前面的示例中，您首先切片了第一，第二和第三行（[0,1,2]），然后切片第一，第二和第四列。
整形和调整大小方法似乎很相似，但是这些操作的输出有所不同。 重新排列数组的形状时，只是输出会临时更改数组的形状，但不会更改数组本身。 调整数组大小时，它将永久更改数组的大小，如果新数组的大小大于旧数组的大小，则新数组元素将填充旧数组的重复副本。 相反，如果新数组较小，则新数组将从旧数组中获取元素，并按索引顺序填充新数组。 请注意，相同的数据可以由不同的ndarray共享。这意味着一个ndarray可以是另一个ndarray的视图。 在这种情况下，对一个数组进行的更改将对其他视图产生影响。
下面的代码给出了一个示例，说明当大小大于或小于原始数组时如何填充新数组元素：
In [88]: x = np.arange(16).reshape(4,4)
x
Out[88]: array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
In [89]: np.resize(x,(2,2))
Out[89]: array([[0, 1],
[2, 3]])
In [90]: np.resize(x,(6,6))
Out[90]: array([[ 0, 1, 2, 3, 4, 5],
[ 6, 7, 8, 9, 10, 11],
[12, 13, 14, 15, 0, 1],
[ 2, 3, 4, 5, 6, 7],
[ 8, 9, 10, 11, 12, 13],
[14, 15, 0, 1, 2, 3]])
本小节的最后一个重要术语是广播，它解释了当 NumPy 具有不同形状时在数组的算术运算中如何表现。 NumPy 有两个广播规则：数组的大小相等或其中之一为 1。如果不满足以下条件之一，则将得到以下两个错误之一：frames are not aligned或operands could not be broadcast together：
In [91]: x = np.arange(16).reshape(4,4)
y = np.arange(6).reshape(2,3)
x+y
--------------------------------------------------------------- ------------
ValueError Traceback (most recent call last)
<ipython-input-102-083fc792f8d9> in <module>()
1 x = np.arange(16).reshape(4,4)
2 y = np.arange(6).reshape(2,3)
----> 3 x+y
12
ValueError: operands could not be broadcast together with shapes (4,4) (2,3)
您可能已经看到可以将两个矩阵的形状分别为(4, 4)和(4,)或(2, 2)和(1, 2)。 第一种情况满足具有一维的条件，因此乘法成为向量*数组，这不会引起任何广播问题：
In [92]: x = np.ones(16).reshape(4,4)
y = np.arange(4)
x*y
Out[92]: array([[ 0., 1., 2., 3.],
[ 0., 1., 2., 3.],
[ 0., 1., 2., 3.],
[ 0., 1., 2., 3.]])
In [93]: x = np.arange(4).reshape(2,2)
x
Out[93]: array([[0, 1],
[2, 3]])
In [94]: y = np.arange(2).reshape(1,2)
y
Out[94]: array([[0, 1]])
In [95]: x*y
Out[95]: array([[0, 1],
[0, 3]])
前面的代码块提供了第二种情况的示例，其中在计算过程中，小数组迭代通过大数组，并且输出扩展到整个数组。 这就是为什么有(4, 4)和(2, 2)输出的原因：在乘法过程中，两个数组都以较大的大小广播。
总结
在本章中，您熟悉了 NumPy 数组操作的基础知识，并刷新了有关基本矩阵操作的知识。 NumPy 是用于 Python 科学堆栈的极其重要的库，它具有用于数组操作的广泛方法。 您已经学习了如何使用多维数组，并涵盖了重要的主题，例如索引，切片，整形，调整大小和广播。 本章的主要目的是让您简要了解有关数字数据集的 NumPy 的工作方式，这将对您的日常数据分析工作有所帮助。
在下一章中，您将学习线性代数的基础知识以及使用 NumPy 完整的实际示例。
二、NumPy 线性代数
数学的主要分类之一是代数，尤其是线性代数专注于线性方程和映射线性空间（即向量空间）。 当我们在向量空间之间创建线性映射时，实际上是在创建称为矩阵的数据结构。 线性代数的主要用途是求解联立线性方程，但也可用于非线性系统的近似。 想象一下您要理解的复杂模型或系统，将其视为非线性模型。 在这种情况下，您可以将问题的复杂的非线性特征简化为联立的线性方程，并且可以在线性代数的帮助下求解它们。
在计算机科学中， 线性代数在机器学习（ML）应用中大量使用。 在 ML 应用中，您处理的是高维数组，可以轻松地将其转换为线性方程式，您可以在其中分析给定空间中要素的相互作用。 想像一下您正在从事图像识别项目并且您的任务是从 MRI 图像中检测出大脑中的肿瘤的情况。 从技术上讲，您的算法应该像医生一样工作，在其中扫描给定的输入并检测大脑中的肿瘤。 医生的优势在于能够发现异常情况。 人类的大脑已经进化了数千年，以解释视觉输入。 无需付出太多努力，人类就可以直观地捕获异常。 但是，对于执行相似任务的算法，您应该尽可能详细地考虑此过程，以了解如何正式表达它，以便机器可以理解。
首先，您应该考虑 MRI 数据如何存储在仅处理 0 和 1 的计算机中。 计算机实际上将像素强度存储在称为矩阵的结构中。 换句话说，您将 MRI 转换为大小为N2的向量，其中每个元素均由像素值组成。 如果此 MRI 大小为512 x 512，则每个像素将是 262,144 像素中的一个点。 因此，您将在此矩阵中进行的任何计算操作都极有可能会使用线性代数原理。 如果该示例不足以证明 ML 中的线性代数的重要性，那么让我们看一下深度学习中的一个流行示例。 简而言之，深度学习是一种算法，该算法使用神经网络结构通过不断更新各层之间神经元连接的权重来学习所需的输出（标签）。 一个简单的深度学习算法的图形表示如下：
神经网络存储各层之间的权重和矩阵中的偏差值。 由于这些参数是您尝试调整到深度学习模型中以最小化损失函数的参数，因此您不断进行计算并更新它们。 通常，机器学习模型需要大量计算，并且需要针对大型数据集进行训练以提供有效的结果。 这就是为什么线性代数是 ML 的基本部分的原因。
在本章中，我们将使用 numpy 库，但请注意，大多数线性代数函数也是通过 scipy 导入的，它们更恰当地属于它。 理想情况下，在大多数情况下，您都可以导入这两个库并执行计算。 scipy 的一个重要特征是拥有功能齐全的线性代数模块。 在本章中，我们强烈建议您阅读 scipy 文档并使用与 scipy 相同的操作进行练习。 这是 scipy 的线性代数模块的链接。
在本章中，我们将介绍以下主题：
- 向量和矩阵数学
- 什么是特征值？我们如何计算它？
- 计算范数和行列式
- 求解线性方程
- 计算梯度
向量和矩阵运算
在上一章中，您练习了向量和矩阵的入门操作。 在本节中，您将练习在线性代数中大量使用的更高级的向量和矩阵运算。 让我们记住关于矩阵操作的点积透视图，以及当您具有 2D 数组时如何使用不同的方法来完成。 以下代码块显示了执行点积计算的替代方法：
In [1]: import numpy as np
a = np.arange(12).reshape(3,2)
b = np.arange(15).reshape(2,5)
print(a)
print(b)
Out[1]:
[[ 0 1]
[ 2 3]
[ 4 5]]
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
In [2]: np.dot(a,b)
Out[2]: array([[ 5, 6, 7, 8, 9],
[15, 20, 25, 30, 35],
[25, 34, 43, 52, 61]])
In [3]: np.matmul(a,b)
Out[3]: array([[ 5, 6, 7, 8, 9],
[15, 20, 25, 30, 35],
[25, 34, 43, 52, 61]])
In [4]: a@b
Out[4]: array([[ 5, 6, 7, 8, 9],
[15, 20, 25, 30, 35],
[25, 34, 43, 52, 61]])
内积和点积在 ML 算法（如监督学习）中非常重要。 让我们回到有关肿瘤检测的示例。 想象我们有三张图像（MRI）：第一张图像有肿瘤（A），第二张图像无肿瘤（B），第三张图像是您要标记为肿瘤或无肿瘤的未知 MRI。下图显示了向量a和b的点积的几何表示：
作为一个非常简单的例子，点积将向您显示这两个向量之间的相似性，因此，如果您未知的 MRI 方向接近向量 A，则您的算法将被归类为具有肿瘤的 MRI。 否则，它将分类为无肿瘤。 如果要在一个函数中乘以两个或多个数组，linalg.multi_dot()非常方便。 进行此操作的另一种方法是将np.dot()嵌套的数组相乘，但您需要知道哪个计算顺序最快，因为在linalg.multi_dot()中此优化是自动进行的。 如前所述，以下代码块使用不同的方法执行相同的点积运算：
In [5]: from numpy.linalg import multi_dot
a = np.arange(12).reshape(4,3)
b = np.arange(15).reshape(3,5)
c = np.arange(25).reshape(5,5)
multi_dot([a, b, c])
Out[5]: array([[ 1700, 1855, 2010, 2165, 2320],
[ 5300, 5770, 6240, 6710, 7180],
[ 8900, 9685, 10470, 11255, 12040],
[12500, 13600, 14700, 15800, 16900]])
In [6]: a.dot(b).dot(c)
Out[6]: array([[ 1700, 1855, 2010, 2165, 2320],
[ 5300, 5770, 6240, 6710, 7180],
[ 8900, 9685, 10470, 11255, 12040],
[12500, 13600, 14700, 15800, 16900]])
如下面的代码所示，即使使用三个小型矩阵，multi_dot()方法也将运行时间减少了 40%。 如果矩阵数量和大小增加，则该时间间隔会大大增加。 因此，您应该使用multi_dot()方法以确保最快的求值顺序。 以下代码将比较这两种方法的执行时间：
In [7]: import numpy as np
from numpy.linalg import multi_dot
import time
a = np.arange(120000).reshape(400,300)
b = np.arange(150000).reshape(300,500)
c = np.arange(200000).reshape(500,400)
start = time.time()
multi_dot([a,b,c])
ft = time.time()-start
print ('Multi_dot tooks', time.time()-start,'seconds.')
start_ft = time.time()
a.dot(b).dot(c)
print ('Chain dot tooks', time.time()-start_ft,'seconds.')
Out[7]:
Multi_dot tooks 0.14687418937683105 seconds.
Chain dot tooks 0.1572890281677246 seconds.
NumPy 的线性代数库中还有另外两种重要的方法，即outer()和inner()。outer方法计算两个向量的外积。inner方法的行为取决于所采用的参数。 如果您有两个向量作为参数，它将产生一个普通的点积，但是当您具有一个更高维的数组时，它会像tensordot()一样返回最后一个轴的积的和。您将在本章后面看到tensordot()。 现在，让我们首先关注outer()和inner()方法。 以下示例将帮助您了解这些函数的作用：
In [8]: a = np.arange(9).reshape(3,3)
b = np.arange(3)
print(a)
print(b)
Out[8]:
[[0 1 2]
[3 4 5]
[6 7 8]]
[0 1 2]
In [9]: np.inner(a,b)
Out[9]: array([ 5, 14, 23])
In [10]: np.outer(a,b)
Out[10]: array([[ 0, 0, 0],
[ 0, 1, 2],
[ 0, 2, 4],
[ 0, 3, 6],
[ 0, 4, 8],
[ 0, 5, 10],
[ 0, 6, 12],
[ 0, 7, 14],
[ 0, 8, 16]])
在上述inner()方法的示例中，数组的第i行与向量产生标量积，总和成为输出数组的第i个元素，因此输出数组的构造如下：
[0x0 + 1x1 + 2x2, 0x3 + 1x4 + 2x5, 0x6 + 1x7 + 2x8] = [5, 14, 23]
在下面的代码块中，我们对同一数组但具有一维执行outer()方法。 如您所见，结果与我们对二维数组所做的完全相同：
In [11]: a = np.arange(9)
np.ndim(a)
Out[11]: 1
In [12]: np.outer(a,b)
Out[12]: array([[ 0, 0, 0],
[ 0, 1, 2],
[ 0, 2, 4],
[ 0, 3, 6],
[ 0, 4, 8],
[ 0, 5, 10],
[ 0, 6, 12],
[ 0, 7, 14],
[ 0, 8, 16]])
outer()方法计算向量的外积，在我们的示例中为 2D 数组。 这不会改变方法的功能，而是将 2D 数组展平为向量并进行计算。
在进行分解之前，本小节介绍的最后一件事是tensordot()方法。 在数据科学项目中，我们通常处理需要进行发现和应用 ML 算法的 n 维数据。 以前，您了解了向量和矩阵。 张量是向量和矩阵的通用数学对象，可以将向量的关系保持在高维空间中。tensordot()方法用于两个张量的收缩； 换句话说，它通过将指定轴上两个张量的乘积相加来减小维数（张量级）。 以下代码块显示了两个数组的tensordot()操作示例：
In [13]: a = np.arange(12).reshape(2,3,2)
b = np.arange(48).reshape(3,2,8)
c = np.tensordot(a,b, axes =([1,0],[0,1]))
print(a)
print(b)
Out[13]:
[[[ 0 1]
[ 2 3]
[ 4 5]]
[[ 6 7]
[ 8 9]
[10 11]]]
[[[ 0 1 2 3 4 5 6 7]
[ 8 9 10 11 12 13 14 15]]
[[16 17 18 19 20 21 22 23]
[24 25 26 27 28 29 30 31]]
[[32 33 34 35 36 37 38 39]
[40 41 42 43 44 45 46 47]]]
In [14]: c
Out[14]: array([[ 800, 830, 860, 890, 920, 950, 980, 1010],
[ 920, 956, 992, 1028, 1064, 1100, 1136, 1172]])
什么是特征值，我们如何计算？
特征值是特征向量的系数。 根据定义，特征向量是一个非零向量，在应用线性变换时仅会按标量因子进行更改。 通常，将线性变换应用于向量时，其跨度（穿过其原点的线）会移动，但是某些特殊向量不受这些线性变换的影响，而是保留在自己的跨度上。 这些就是我们所说的特征向量。 线性变换仅在将向量与标量相乘时才通过拉伸或挤压它们来影响它们。 该标量的值称为特征值。 假设我们有一个 A 矩阵，该矩阵将用于线性变换。 我们可以用以下数学表达式表示特征值和特征向量：
在这里，
是特征向量，
是特征值。 在等式的左侧，向量通过矩阵进行变换，结果只是同一向量的标量形式。 另外，请注意，左侧实际上是矩阵向量乘法，而右侧是标量向量乘法。 为了使边的乘法类型相同，我们可以将 与恒等矩阵相乘，这不会改变右侧。 但是将 从标量更改为矩阵，并且两边都是矩阵向量乘法：
如果我们减去右侧并分解为 
，您将具有以下等式：
因此，我们的新矩阵将如下所示：
如您所知，您正在尝试计算0以外的特征向量的特征值； （如果 
= 0，则没有意义），因此特征向量（
）不能为0。 因此，您正在尝试求解以下等式：
前面的公式正在计算矩阵的行列式，因为只有一个条件，即非零矩阵等于0，即矩阵的行列式等于零。 给定矩阵的行列式为零意味着使用该矩阵进行的转换会将所有内容压缩为较小的维度。 在下一部分中，您将详细了解行列式。 此处的真正目的是找到使行列式为零并将空间压缩为较小维度的 。 找到 后，我们可以根据以下公式计算特征向量（
）：
在 ML 算法中，您需要处理大量维度。 主要问题不是规模很大，而是算法与它们的兼容性和性能。 例如，在 PCA（主成分分析）中，您尝试发现大小最有意义的线性组合。 PCA 的主要思想是减少数据集的规模，同时最大程度地减少信息丢失。 此处的信息丢失实际上是您特征的差异。 假设您具有五个特征，并且具有五个示例的类标签，如下所示：
| 特征 1 | 特征 2 | 特征 3 | 特征 4 | 特征 5 | 类别 |
|---|---|---|---|---|---|
| 1.45 | 42 | 54 | 1.001 | 1.05 | 狗 |
| 2 | 12 | 34 | 1.004 | 24.1 | 猫 |
| 4 | 54 | 10 | 1.004 | 13.4 | 狗 |
| 1.2 | 31 | 1 | 1.003 | 42.1 | 猫 |
| 5 | 4 | 41 | 1.003 | 41.4 | 狗 |
在上表中，在特征 4 中，类别标签是狗还是猫实际上并没有明显的区别。 在我的分析中，此特征将是多余的。 这里的主要目标是使各个类之间的特征保持很大的差异，因此特征值在我的决定中起着重要的作用。
在以下示例中，我们将使用一个名为scikit-learn的附加库，以便加载数据集并通过使用该库来执行 PCA。 Scikit-learn 是一个免费的 python 机器学习库，它基于 SciPy 构建。 它具有用于许多机器学习算法（例如分类，回归，聚类）的许多内置函数和模块。 SVM，DBSCAN，K 均值等。
Scikit-learn 与 numpy，scipy 和 pandas 具有很好的兼容性。 您可以将sklearn中的 numpy 数组用作数据结构。 此外，您还可以使用sklearn_pandas将 scikit-learn 管道输出转换为 pandas 数据帧。 Scikit-learn 还通过auto_ml和auto-sklearn库支持自动化机器学习算法。sklearn是在 python 中键入 scikit-learn 名称的一种方法。 这就是为什么在导入 skicit-learn 和使用函数时使用sklearn的原因。
现在是时候进行实际练习，以了解特征值和特征向量在 PCA 中的用法和重要性。 在下面的代码中，您将 PCA 与 NumPy 一起导入并应用，然后将您的结果与sklearn是用于验证的内置方法进行比较。 您将使用 sklearn 数据集库中的乳腺癌数据集。 让我们先导入所需的库和数据集，然后再对数据进行标准化。 标准化非常重要，有时甚至对于 scikit-learn 之类的某些 ML 库中的估计器来说也是必需的。 在我们的示例中，使用StandardScaler()方法来标准化特征； 主要目的是具有标准正态分布数据中的特征（均值= 0和单位方差）。fit_transform()方法用于转换原始数据，其形式为分布的平均值为 0，标准差为 1。它将计算所需的参数并应用转换。 当我们使用StandardScaler()时，此参数将为 
和 
。 请记住，标准化不会从我们的原始数据集中产生正态分布的数据。 只是重新缩放了数据，我们的平均值为零，标准差为 1：
In [15]: import numpy as np
from sklearn import decomposition, datasets
from sklearn.preprocessing import StandardScaler
data = datasets.load_breast_cancer()
cancer = data.data
cancer = StandardScaler().fit_transform(cancer)
cancer.shape
Out[15]: (569, 30)
在前面的代码中检查数据形状时，它显示数据由569行和30列组成。 您也可以在标准化开始和标准化之后比较数据，以便更清楚地了解原始数据是如何转换的。 以下代码块为您展示了癌症数据中列转换的示例：
In [16]: before_transformation = data.data
before_transformation[:10,:1]
Out[16]: array([[ 17.99],
[ 20.57],
[ 19.69],
[ 11.42],
[ 20.29],
[ 12.45],
[ 18.25],
[ 13.71],
[ 13\. ],
[ 12.46]])
In [17]: cancer[:10,:1]
Out[17]: array([[ 1.09706398],
[ 1.82982061],
[ 1.57988811],
[-0.76890929],
[ 1.75029663],
[-0.47637467],
[ 1.17090767],
[-0.11851678],
[-0.32016686],
[-0.47353452]])
从结果中可以看到，原始数据已转换为标准化版本的形式。 通过以下公式计算得出：
![]()
：正在标准化的值（原始数据中的值）![]()
：分布的平均值![]()
：分布的标准差![]()
：标准化值
：正在标准化的值（原始数据中的值）
：分布的平均值
：分布的标准差
：标准化值转换数据后，您将如下计算协方差矩阵，以便使用np.linalg.eig()方法计算特征值和特征向量，然后将其用于分解：
In [18]: covariance_matrix = np.cov(cancer,rowvar=False)
covariance_matrix.shape
Out[18]: (30, 30)
In [19]: eig_val_cov, eig_vec_cov = np.linalg.eig(covariance_matrix)
eig_pairs = [(np.abs(eig_val_cov[i]), eig_vec_cov[:,i]) for i in
range(len(eig_val_cov))]
正如您在前面的代码中已经注意到的那样，我们计算为所有要素构造的协方差矩阵。 由于我们在数据集中具有 30 个特征，因此协方差矩阵是形状为(30, 30)的 2D 数组。 以下代码块按降序对特征值进行排序：
In [20]: sorted_pairs = eig_pairs.sort(key=lambda x: x[0], reverse=True)
for i in eig_pairs:
print(i[0])
Out[20]:
13.3049907944
5.70137460373
2.82291015501
1.98412751773
1.65163324233
1.2094822398
0.676408881701
0.47745625469
0.417628782108
0.351310874882
0.294433153491
0.261621161366
0.241782421328
0.157286149218
0.0943006956011
0.0800034044774
0.0595036135304
0.0527114222101
0.049564700213
0.0312142605531
0.0300256630904
0.0274877113389
0.0243836913546
0.0180867939843
0.0155085271344
0.00819203711761
0.00691261257918
0.0015921360012
0.000750121412719
0.000133279056663
您需要以降序对特征值进行排序，以确定要降低维度工作空间要删除的特征向量。 正如您在前面的排序列表中所看到的那样，具有较高特征值的前两个特征向量具有有关数据分布的最多信息。 因此，其余的将针对低维子空间删除：
In [21]: matrix_w = np.hstack((eig_pairs[0][1].reshape(30,1), eig_pairs[1] [1].reshape(30,1)))
matrix_w.shape
transformed = matrix_w.T.dot(cancer.T)
transformed = transformed.T
transformed[0]
Out[21]: array([ 9.19283683, 1.94858307])
In [22]: transformed.shape
Out[22]: (569, 2)
在前面的代码块中，前两个特征向量将水平堆叠，因为它们将用于矩阵乘法，以将新维度的数据转换到新的子空间上。 最终数据从(569,30)转换为(569,2)矩阵，这意味着28特征在 PCA 过程中下降了：
In [23]: import numpy as np
from sklearn import decomposition
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
pca = decomposition.PCA(n_components=2)
x_std = StandardScaler().fit_transform(cancer)
pca.fit_transform(x_std)[0]
Out[23]: array([ 9.19283683, 1.94858307])
另一方面，其他库中有内置函数可以执行与您相同的操作。 在 scikit-learn 中，有许多内置方法可用于 ML 算法。 如您在前面的代码块中看到的，同一 PCA 是通过三行代码和两种方法执行的。 不过，本示例的目的是向您展示特征值和特征向量的重要性，因此，本书展示了使用 NumPy 进行 PCA 的普通方法。
计算范数和行列式
本小节将介绍线性代数中的两个重要值，即范数和行列式。 简而言之，范数给出向量的长度。 最常用的范数是 L2 范数，也称为欧几里得范数。 正式地，Lp 的范数计算如下：
L0 范数实际上是向量的基数。 您可以通过仅计算非零元素的总数来计算它。 例如，向量A = [2,5,9,0]包含三个非零元素，因此||A||[0] = 3。 以下代码块显示了与 numpy 相同的范数计算：
In [24]: import numpy as np
x = np.array([2,5,9,0])
np.linalg.norm(x,ord=0)
Out[24]: 3.0
在 NumPy 中，您可以使用linalg.norm()方法来计算向量的范数。 第一个参数是输入数组，而ord参数用于范数的顺序。 L1 范数也称为出租车范数或曼哈顿范数。 它通过计算直线距离来计算向量的长度，因此||A||[1] = (2 + 5 + 9)，如此||A||[1] = 16。以下代码块显示了相同的 numpy 范数计算：
In [25]: np.linalg.norm(x,ord=1)
Out[25]: 16.0
L1 范数的用途之一是用于计算平均绝对误差（MAE），如下式所示：
L2 范数是最受欢迎的范数。 它通过应用勾股定理来计算向量的长度，因此||A||[2] = 
，所以||A||[2] = 10.48：
In [26]: np.linalg.norm(x,ord=2)
Out[26]: 10.488088481701515
L2 范数的众所周知的应用之一是均方误差（MSE）计算。 类似地，由于矩阵由向量组成，因此范数给出了长度或大小，但解释和计算略有不同。 在前面的章节中，您学习了向量矩阵乘法。 在结果中将矩阵与向量相乘时，将拉伸向量。 矩阵的范数揭示了该矩阵可以拉伸向量的程度。 让我们看看如何为矩阵A计算 L1 和 L∞ 范数的模数：
假设您有一个m x n矩阵。||A||[1]的计算如下：
||A||[1] = 
因此结果将如下所示：
||A||[1] = max(3 + |-2| + 1; 7 + |-5| + 3; 6 + 4 + |-14|) = max(6, 15, 24) = 24
让我们使用linalg.norm()方法，用 NumPy 计算同一数组的范数：
In [27]: a = np.array([3,7,6,-2,-5,4,1,3,-14]).reshape(3,3)
In [28]: a
Out[28]: array([[ 3, 7, 6],
[ -2, -5, 4],
[ 1, 3, -14]])
In [29]: np.linalg.norm(a, ord=1)
Out[29]: 24.0
In [30]: np.linalg.norm(a, np.inf)
Out[30]: 18.0
在前面的计算中||A||[1]首先按列计算最大元素，并在所有列中给出最大结果，这将是一阶矩阵的范数。 对于 L∞ 元素计算按行执行，并在所有获得无穷范数范数矩阵的行中给出最大值。 计算结果如下：
||A||[∞] = max(3 + 7 + 6; |-2| + |-5| + 4; 1 + 3 + |-14|) = max(16, 11, 18) = 18
为了验证和使用 NumPy 函数，与我们在向量范数计算中所做的一样，使用linalg.norm()方法进行了相同的计算。 一阶和无限阶的计算比 Lp 范数（其中p > 2）的常规计算相对更直接，例如欧几里德/弗罗宾纽斯范数（其中p = 2）。 下面的公式显示了p范数的形式公式，只需将p替换为 2，就可以为给定数组制定欧几里得/弗罗比纽斯范数：
对于p = 2的特殊情况，它变为：
以下代码块显示了数组a的 L2 范数计算：
In [31]: np.linalg.norm(a, ord=2)
Out[31]: 15.832501006406099
在 NumPy 中，可以通过在linalg.norm()方法中将order参数设置为 2 来计算欧几里德/弗罗比纽斯范数，如您在前面的代码中所见。 在 ML 算法中，特征在特征空间的距离计算中大量使用。 例如，在模式识别中， 范数计算用于 K 最近邻算法（KNN），以便为连续变量或离散变量创建距离度量。 同样，范数对于 K 均值聚类中的距离度量也非常重要。 最受欢迎的顺序是曼哈顿范数和欧几里得范数。
线性代数中的另一个重要概念是计算矩阵的行列式。 根据定义，行列式是线性变换中给定矩阵的比例因子。 在上一节中，当您计算特征值和特征向量时，将给定矩阵乘以特征向量，并假定它们的行列式将为零。 换句话说，您假设将特征向量乘以矩阵后，会将矩阵展平为较低的维度。2 x 2和3 x 3矩阵的行列式如下：
在 NumPy 中，您可以使用linalg.det()方法来计算矩阵的行列式。 让我们用上式计算一个2 x 2和3 x 3矩阵行列式的示例，然后用 NumPy 交叉验证我们的结果：
因此计算如下：
对于我们有3 x 3矩阵的情况，它变为：
让我们使用linalg.det()方法，用 NumPy 计算相同数组的行列式：
In [32]: A= np.array([2,3,1,4]).reshape(2,2)
In [33]: A
Out[33]: array([[2, 3],
[1, 4]])
In [34]: np.linalg.det(A)
Out[34]: 5.0
In [35]: B= np.array([2,3,5,1,4,8,5,6,2]).reshape(3,3)
In [36]: B
Out[36]: array([[2, 3, 5],
[1, 4, 8],
[5, 6, 2]])
In [37]: np.linalg.det(B)
Out[37]: -36.0
转换的决定因素实际上显示了将扩展或压缩多少体积的因素。 如果矩阵的行列式等于2，则意味着此矩阵的变换将使体积增加2。 如果要进行链乘法，则还可以计算变换后音量的变化量。 假设您将两个矩阵相乘，这意味着将有两个转换。 如果det(A)= 2并且det(B)= 3，则总转换因子将乘以 6 作为det(AB) = det(A)det(B)。
最后，本章将为您的 ML 模型介绍一个非常有用的值，该值称为迹。 根据定义，迹线是矩阵对角线元素的总和。 在 ML 模型中，大多数情况下，您将使用几种回归模型来解释数据。 这些模型中的某些很有可能在某种程度上解释了您的数据质量，因此在这种情况下，您总是倾向于使用更简单的模型。 每当您必须进行权衡时，跟踪值对于量化复杂性就变得非常有价值。 以下代码块显示了使用 numpy 进行 2D 和 3D 矩阵的跟踪计算：
In [38]: a = np.arange(9).reshape(3,3)
In [39]: a
Out[39]: array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [40]: np.trace(a)
Out[40]: 12
In [41]: b = np.arange(27).reshape(3,3,3)
In [42]: b
Out[42]: array([[[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8]],
[[ 9, 10, 11],
[12, 13, 14],
[15, 16, 17]],
[[18, 19, 20],
[21, 22, 23],
[24, 25, 26]]])
In [43]: np.trace(b)
Out[43]: array([36, 39, 42])
在 NumPy 中，可以使用trace()方法来计算矩阵的跟踪值。 如果您的矩阵是二维的，那么轨迹就是对角线的总和。 如果矩阵的维数大于2，则迹线是沿对角线的和数组。 在前面的示例中，矩阵B具有三个维度，因此跟踪数组的构造如下：
在本小节中，您学习了如何在 ML 算法中计算范数，行列式，跟踪和用法。 主要目的是学习这些概念并熟悉 NumPy 线性代数库。
求解线性方程
在本节中，您将学习如何使用linalg.solve()方法求解线性方程。 当您有线性方程要求解时，例如形式 
，在简单情况下，您只需计算A的反函数即可。然后将其乘以B即可得到解，但是当A具有高维数时，这使得计算起来很难进行计算A的倒数。 让我们从具有三个未知数的三个线性方程式的示例开始，如下所示：
因此，这些方程式可以用矩阵形式化如下：
然后，我们的问题是解决 
。 我们可以不使用linalg.solve()，而使用普通 NumPy 计算解决方案。颠倒A矩阵后，您将与B相乘以获得x的结果。在下面的代码块中，我们计算A和B的逆矩阵的点积，以便计算 
：
In [44]: A = np.array([[2, 1, 2], [3, 2, 1], [0, 1, 1]])
A
Out[44]: array([[2, 1, 2],
[3, 2, 1],
[0, 1, 1]])
In [45]: B = np.array([8,3,4])
B
Out[45]: array([8, 3, 4])
In [46]: A_inv = np.linalg.inv(A)
A_inv
Out[46]: array([[ 0.2, 0.2, -0.6],
[-0.6, 0.4, 0.8],
[ 0.6, -0.4, 0.2]])
In [47]: np.dot(A_inv,B)
Out[47]: array([-0.2, -0.4, 4.4])
最后，得到a = -0.2，b = -0.4和c = 4.4的结果。 现在，让我们使用linalg.solve()执行相同的计算，如下所示：
In [48]: A = np.array([[2, 1, 2], [3, 2, 1], [0, 1, 1]])
B = np.array([8,3,4])
x = np.linalg.solve(A, B)
x
Out[48]: array([-0.2, -0.4, 4.4])
为了检查我们的结果，我们可以使用allclose()函数，该函数用于按元素比较两个数组：
In [49]: np.allclose(np.dot(A, x), B)
Out[49]: True
linalg.lstsq()方法是求解线性方程组的另一个重要函数，它返回最小二乘解。 此函数将返回回归线的参数。 回归线的主要思想是最小化每个数据点到回归线的距离平方和。 距离的平方和实际上量化了回归线的总误差。 距离越大，误差越大。 因此，我们正在寻找可以最大程度减少此误差的参数。 为了可视化我们的线性回归模型，我们将使用一个非常流行的 Python 二维绘图库matplotlib。 以下代码块在矩阵中运行最小二乘解，并返回权重和偏差：
In [50]: from numpy import arange,array,ones,linalg
from pylab import plot,show
x = np.arange(1,11)
A = np.vstack([x, np.ones(len(x))]).T
A
Out[50]: array([[ 1., 1.],
[ 2., 1.],
[ 3., 1.],
[ 4., 1.],
[ 5., 1.],
[ 6., 1.],
[ 7., 1.],
[ 8., 1.],
[ 9., 1.],
[ 10., 1.]])
In [51]: y = [5, 6, 6.5, 7, 8,9.5, 10, 10.4,13.1,15.5]
w = linalg.lstsq(A,y)[0]
w
Out[51]: array([ 1.05575758, 3.29333333])
In [52]: line = w[0]*x+w[1]
plot(x,line,'r-',x,y,'o')
show()
前面代码的输出如下：
上图给出了适合我们数据的回归线拟合结果。 该模型生成可用于预测或预测的线性线。 尽管线性回归有很多假设（恒定方差，误差的独立性，线性等），但它仍然是建模数据点之间线性关系的最常用方法。
计算梯度
当您有一条直线时，您可以使用导数，以便导数显示该线的斜率。 当函数中有多个变量时，梯度是导数的泛化，因此，梯度的结果实际上是向量函数，而不是导数中的标量值。 ML 的主要目标实际上是找到适合您数据的最佳模型。 您可以通过将损失函数或目标函数最小化来求值最佳均值。 梯度用于查找系数或函数，以最大程度地减少损失或成本函数。 查找最佳点的一种众所周知的方法是采用目标函数的导数，然后将其设置为零以找到模型系数。 如果系数不止一个，则它将成为梯度而不是导数，并且将成为向量方程式而不是标量值。 您可以在每个点上将梯度解释为向量，该梯度将指向函数的下一个局部最小值。 有一种非常流行的优化技术，可以通过计算每个点的梯度并沿该方向移动系数来寻找函数的最小值，以寻求最小值。 该方法称为梯度下降。 在 NumPy 中，gradient()函数用于计算数组的梯度：
In [53]: a = np.array([1,3, 6, 7, 11, 14])
gr = np.gradient(a)
gr
Out[53]: array([ 2\. , 2.5, 2\. , 2.5, 3.5, 3\. ])
让我们看看如何计算此梯度：
gr[0] = (a[1]-a[0])/1 = (3-1)/1 = 2
gr[1] = (a[2]-a[0])/2 = (6-1)/2 = 2.5
gr[2] = (a[3]-a[1])/2 = (7-3)/2 = 2
gr[3] = (a[4]-a[2])/2 = (11-6)/2 = 2.5
gr[4] = (a[5]-a[3])/2 = (14-7)/2 = 3.5
gr[5] = (a[5]-a[4])/1 = (14-11)/1 = 3
在前面的代码块中，我们计算了一维数组的梯度。 让我们添加另一个维度，并查看计算如何更改，如下所示：
[54]: a = np.array([1,3, 6, 7, 11, 14]). reshape(2,3)
gr = np.gradient(a)
gr
Out[54]: [array([[ 6., 8., 8.],
[ 6., 8., 8.]]), array([[ 2\. , 2.5, 3\. ],
[ 4\. , 3.5, 3\. ]])]
对于二维数组，与前面的代码一样，按列和按行计算梯度。 因此，结果将是两个数组的返回。 第一个数组代表行方向，第二个数组代表列方向。
总结
在本章中，我们介绍了线性代数的向量和矩阵运算。 我们研究了高级矩阵运算，尤其是特征点运算。 您还了解了的特征值和特征向量，然后在主成分分析（PCA）中实践了它们的用法。 此外，我们介绍了范数和行列式计算，并提到了它们在 ML 中的重要性和用法。 在最后两个小节中，您学习了如何将线性方程式转换为矩阵并求解它们，并研究了梯度的计算和重要性。
在下一章中，我们将使用 NumPy 统计数据进行解释性数据分析，以探索 2015 年美国住房数据。
三、使用 NumPy 统计函数对波士顿住房数据进行探索性数据分析
探索性数据分析（EDA）是数据科学项目的重要组成部分（如图数据科学过程所示）。 尽管在将任何统计模型或机器学习算法应用于数据之前这是非常重要的一步，但许多从业人员经常忽略或低估了它：
John Wilder Tukey 于 1977 年通过他的著作探索性数据分析促进了探索性数据分析。 在他的书中，他指导统计学家通过使用几种不同的视觉手段对统计数据集进行统计分析，这将帮助他们制定假设。 此外，在确定关键数据特征并了解有关数据的哪些问题后，EDA 还可以用于为高级建模准备分析。 在较高的层次上，EDA 是对您的数据的无假设探索，它利用了定量方法，使您可以可视化结果，从而可以识别模式，异常和数据特征。 在本章中，您将使用 NumPy 的内置统计方法执行探索性数据分析。
本章将涵盖以下主题：
- 加载和保存文件
- 探索数据集
- 查看基本统计量
- 计算平均值和方差
- 计算相关性
- 计算直方图
加载和保存文件
在本节中，您将学习如何加载/导入数据并保存。 加载数据的方式有很多，正确的方式取决于您的文件类型。 您可以加载/导入文本文件，SAS/Stata 文件，HDF5 文件以及许多其他文件。 HDF（分层数据格式）是一种流行的数据格式，用于存储和组织大量数据，并且在处理多维同构数组时非常有用。 例如，Pandas 库有一个非常方便的类，名为HDFStore，您可以在其中轻松使用 HDF5 文件。 在从事数据科学项目时，您很可能会看到许多这类文件，但是在本书中，我们将介绍最受欢迎的文件，例如 NumPy 二进制文件，文本文件（.txt）和逗号分隔值（.csv）文件。
如果内存和磁盘上有大量数据要管理，则可以使用bcolz库。 它提供了列式和压缩数据容器。 bcolz对象称为chunks，您在其中将整个数据压缩为位，然后在查询时部分解压缩。 压缩数据后，它会非常高效地使用存储。bcolz对象也提高了数据获取性能。 如果您对该库的性能感兴趣。 您可以在其官方 GitHub 存储库 上检查查询和速度比较。
使用数组时，通常在完成使用后将它们另存为 NumPy 二进制文件。 原因是您还需要存储数组形状和数据类型。 重新加载数组时，您希望 NumPy 记住它，并且您可以从上次中断的地方继续工作。 此外，即使您在具有不同架构的另一台计算机上打开文件， NumPy 二进制文件也可以存储有关数组的信息。 在 NumPy 中，load()，save()，savez()和savez_compressed()方法可帮助您加载和保存 NumPy 二进制文件，如下所示：
In [1]: import numpy as np
In [2]: example_array = np.arange(12).reshape(3,4)
In [3]: example_array
Out[3]: array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
In [4]: np.save('example.npy',example_array)
In [5]: d = np.load('example.npy')
In [6]: np.shape(d)
Out[6]: (3, 4)
In [7]: d
Out[7]: array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
在前面的代码中，我们执行以下步骤来练习将数组保存为二进制文件，以及如何在不影响其形状的情况下将其加载回：
- 创建形状为
(3, 4)的数组 - 将数组另存为二进制文件
- 加载回数组
- 检查形状是否仍然相同
同样，您可以使用savez()函数将多个数组保存到单个文件中。 如果要将文件另存为压缩的 NumPy 二进制文件，则可以按以下方式使用savez_compressed()：
In [8]: x = np.arange(10)
y = np.arange(12)
np.savez('second_example.npz',x, y)
npzfile = np.load('second_example.npz')
npzfile.files
Out[8]: ['arr_0', 'arr_1']
In [9]: npzfile['arr_0']
Out[9]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [10]: npzfile['arr_1']
Out[10]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
In [11]: np.savez_compressed('compressed_example.npz', first_array = x , second_array = y)
npzfile = np.load('compressed_example.npz')
npzfile.files
Out[11]: ['first_array', 'second_array']
In [12]: npzfile['first_array']
Out[12]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [13]: npzfile['second_array']
Out[13]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
当您将多个数组保存在单个文件中时，如果您提供关键字参数如first_array=x，则将使用该名称保存数组。 否则，默认情况下，您的第一个数组将被赋予一个变量名，例如arr_0。 NumPy 具有称为loadtxt()的内置函数，用于从文本文件加载数据。 让我们加载一个.txt文件，该文件由一些整数和一个标头组成，该标头是文件顶部的字符串：
在前面的代码中，由于无法将字符串转换为浮点数，因此会出现误差，这实际上意味着您正在读取非数字值。 这样做的原因是该文件在顶部包含作为标题的字符串以及数值。 如果知道标题有多少行，则可以使用skiprows参数跳过这些行，如下所示：
In [15]: a = np.loadtxt("My_file.txt", delimiter='\t', skiprows=4)
a
Out[15]: array([ 0., 1., 2., 3., 4., 5., 6., 7., 8.,
9., 10., 11., 12., 13., 14., 15., 16., 17.,
18., 19., 20., 21., 22., 23., 24., 25., 26.,
27., 28., 29., 30., 31., 32., 33., 34., 35.,
36., 37., 38., 39., 40., 41., 42., 43., 44.,
45., 46., 47., 48., 49., 50., 51., 52., 53.,
54., 55., 56., 57., 58., 59., 60., 61., 62.,
63., 64., 65., 66., 67., 68., 69., 70., 71.,
72., 73., 74., 75., 76., 77., 78., 79., 80.,
81., 82., 83., 84., 85., 86., 87., 88., 89.,
90., 91., 92., 93., 94., 95., 96., 97., 98.,
99., 100., 101., 102., 103., 104., 105., 106., 107.,
108., 109., 110., 111., 112., 113., 114., 115., 116.,
117., 118., 119., 120., 121., 122., 123., 124., 125.,
126., 127., 128., 129., 130., 131., 132., 133., 134.,
135., 136., 137., 138., 139., 140., 141., 142., 143.,
144., 145., 146., 147., 148., 149., 150., 151., 152.,
153., 154., 155., 156., 157., 158., 159., 160., 161.,
162., 163., 164., 165., 166., 167., 168., 169., 170.,
171., 172., 173., 174., 175., 176., 177., 178., 179.,
180., 181., 182., 183., 184., 185., 186., 187., 188.,
189., 190., 191., 192., 193., 194., 195., 196., 197.,
198., 199., 200., 201., 202., 203., 204., 205., 206.,
207., 208., 209., 210., 211., 212., 213., 214., 215.,
216., 217., 218., 219., 220., 221., 222., 223., 224.,
225., 226., 227., 228., 229., 230., 231., 232., 233.,
234., 235., 236., 237., 238., 239., 240., 241., 242.,
243., 244., 245., 246., 247., 248., 249.])
另外，您可以使用genfromtxt()并将标头的转换默认设置为nan。 然后，您可以检测到标题占用了多少行，并使用skip_header参数跳过它们，如下所示：
In [16]: b = np.genfromtxt("My_file.txt", delimiter='\t')
b
Out[16]: array([ nan, nan, nan, nan, 0., 1., 2., 3., 4.,
5., 6., 7., 8., 9., 10., 11., 12., 13.,
14., 15., 16., 17., 18., 19., 20., 21., 22.,
23., 24., 25., 26., 27., 28., 29., 30., 31.,
32., 33., 34., 35., 36., 37., 38., 39., 40.,
41., 42., 43., 44., 45., 46., 47., 48., 49.,
50., 51., 52., 53., 54., 55., 56., 57., 58.,
59., 60., 61., 62., 63., 64., 65., 66., 67.,
68., 69., 70., 71., 72., 73., 74., 75., 76.,
77., 78., 79., 80., 81., 82., 83., 84., 85.,
86., 87., 88., 89., 90., 91., 92., 93., 94.,
95., 96., 97., 98., 99., 100., 101., 102., 103.,
104., 105., 106., 107., 108., 109., 110., 111., 112.,
113., 114., 115., 116., 117., 118., 119., 120., 121.,
122., 123., 124., 125., 126., 127., 128., 129., 130.,
131., 132., 133., 134., 135., 136., 137., 138., 139.,
140., 141., 142., 143., 144., 145., 146., 147., 148.,
149., 150., 151., 152., 153., 154., 155., 156., 157.,
158., 159., 160., 161., 162., 163., 164., 165., 166.,
167., 168., 169., 170., 171., 172., 173., 174., 175.,
176., 177., 178., 179., 180., 181., 182., 183., 184.,
185., 186., 187., 188., 189., 190., 191., 192., 193.,
194., 195., 196., 197., 198., 199., 200., 201., 202.,
203., 204., 205., 206., 207., 208., 209., 210., 211.,
212., 213., 214., 215., 216., 217., 218., 219., 220.,
221., 222., 223., 224., 225., 226., 227., 228., 229.,
230., 231., 232., 233., 234., 235., 236., 237., 238.,
239., 240., 241., 242., 243., 244., 245., 246., 247.,
248., 249.])
同样，您可以使用loadtxt()，genfromtxt()和savetxt()函数加载和保存.csv文件。 您唯一需要记住的是使用逗号作为定界符，如下所示：
In [17]: data_csv = np.loadtxt("MyData.csv", delimiter=',')
In [18]: data_csv[1:3]
Out[18]: array([[ 0.21982, 0.31271, 0.66934, 0.06072, 0.77785, 0.59984,
0.82998, 0.77428, 0.73216, 0.29968],
[ 0.78866, 0.61444, 0.0107 , 0.37351, 0.77391, 0.76958,
0.46845, 0.76387, 0.70592, 0.0851 ]])
In [19]: np.shape(data_csv)
Out[19]: (15, 10)
In [20]: np.savetxt('MyData1.csv',data_csv, delimiter = ',')
In [21]: data_csv1 = np.genfromtxt("MyData1.csv", delimiter = ',')
In [22]: data_csv1[1:3]
Out[22]: array([[ 0.21982, 0.31271, 0.66934, 0.06072, 0.77785, 0.59984,
0.82998, 0.77428, 0.73216, 0.29968],
[ 0.78866, 0.61444, 0.0107 , 0.37351, 0.77391, 0.76958,
0.46845, 0.76387, 0.70592, 0.0851 ]])
In [23]: np.shape(data_csv1)
Out[23]: (15, 10)
加载.txt文件时，它们会默认返回带有值的numpy数组，如后面的代码中所示：
In [24]: print (type(a))
print (type(b))
<class 'numpy.ndarray'>
<class 'numpy.ndarray'>
您可以使用tolist()方法将数据从数组更改为列表，然后使用savetxt()使用不同的分隔符将其保存为新文件，如下所示：
In [25]: c = a.tolist()
c
Out[25]: [0.0,
1.0,
2.0,
3.0,
4.0,
5.0,
6.0,
7.0,
8.0,
9.0,
10.0,
11.0,
12.0,
13.0,
14.0,
15.0,
16.0,
17.0,
18.0,
...
In [26]: np.savetxt('My_List.txt',c, delimiter=';')
In [27]: myList = np.loadtxt("My_List.txt", delimiter=';')
type(myList)
Out[27]: numpy.ndarray
将列表保存到My_List.txt后，可以使用loadtxt()加载此文件，它将再次作为numpy数组返回。 如果要返回数组的字符串表示形式，可以使用array_str()，array_repr()或array2string()方法，如下所示：
In [28]: d = np.array_str(a,precision=1)
d
Out[28]: '[ 0\. 1\. 2\. 3\. 4\. 5\. 6\. 7\. 8\. 9\. 10\. 11.\n 12\. 13\. 14\. 15\. 16\. 17\. 18\. 19\. 20\. 21\. 22\. 23.
24\. 25\. 26\. ... ]'
即使array_str()和array_repr()看起来一样，array_str()返回数组内部数据的字符串表示形式，而array_repr()实际上返回数组的字符串表示形式。 因此，array_repr()返回有关数组类型及其数据类型的信息。 这两个函数在内部都使用array2string()，这实际上是字符串格式化函数最灵活的版本，因为它具有更多可选参数。 以下代码块使用load_boston()方法加载波士顿房屋数据集：
In [29]: from sklearn.datasets import load_boston
dataset = load_boston()
dataset
在本章中，您将对sklearn.datasets包中的示例数据集进行探索性数据分析。 该数据集是关于波士顿房价的。 在前面的代码中，load_boston()函数是从sklearn.datasets包中导入的，如您所见，它返回具有属性DESCR，data，feature_names和target的类似字典的对象 ]。 这些属性的详细信息如下：
| 属性 | 解释 |
|---|---|
DESCR |
数据集的完整描述 |
data |
特征列 |
feature_names |
特征名称 |
target |
标签数据 |
在本节中，您学习了如何使用numpy函数加载和保存文件。 在下一部分中，您将探索波士顿住房数据集。
探索我们的数据集
在本节中，您将探索数据集并对其进行质量检查。 您将检查数据形状是什么，以及它的数据类型，所有缺失值/ NaN，拥有多少个特征列以及每个列代表什么。 让我们开始加载数据并进行探索：
In [30]: from sklearn.datasets import load_boston
dataset = load_boston()
samples,label, feature_names = dataset.data , dataset.target , dataset.feature_names
In [31]: samples.shape
Out[31]: (506, 13)
In [32]: label.shape
Out[32]: (506,)
In [33]: feature_names
Out[33]: array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'],
dtype='<U7')
在前面的代码中，您加载了数据集并解析了数据集的属性。 这表明我们拥有具有13函数的506样本，并且具有506标签（回归目标）。 如果要阅读数据集的描述，可以使用print(dataset.DESCR)。 由于此代码的输出太长，无法放置在此处，因此您可以在以下屏幕截图中查看各种特征及其描述：
如第一章所示，您可以使用dtype检查数组的数据类型。 如下面的代码所示，示例的每一列中都有一个数字float64数据类型和一个标签。 检查数据类型是非常重要的一步，您可能会意识到类型和列描述之间存在一些不一致，或者，如果您认为使用不太精确的目标值仍然可以实现，则可以通过更改其数据类型来减少数组的内存大小：
In [35]: print(samples.dtype)
print(label.dtype)
float64
float64
缺失值处理在 Python 软件包中略有不同。numpy库没有缺失值。 如果您的数据集中缺少值，则导入后它们将转换为NaN。 NumPy 中非常常见的方法是使用掩码数组以忽略NaN值，这在第一章中已向您展示：
In [36]: np.isnan(samples)
Out[36]: array([[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]], dtype=bool)
In [37]: np.isnan(np.sum(samples))
Out[37]: False
In [38]: np.isnan(np.sum(label))
Out[38]: False
要针对数据中的NaN值逐元素进行测试，可以使用isnan()方法。 此方法将返回一个布尔数组。 对于大型数组而言，检测它是否返回true可能很麻烦。 在这种情况下，您可以将数组的np.sum()用作isnan()的参数，这样它将为结果返回单个布尔值。
在本节中，您将深入研究数据并进行常规质量检查。 在下一节中，我们将继续进行基本统计。
查看基本统计量
在本节中，您将从统计分析的第一步开始，计算数据集的基本统计信息。 即使 NumPy 的内置统计函数有限，我们也可以将其与 SciPy 结合使用。 在开始之前，让我们描述一下我们的分析将如何进行。 所有特征列和标签列都是数字，但您可能已经注意到 Charles River 虚拟变量（CHAS）列具有二进制值(0, 1)，这意味着它实际上是根据分类数据编码的。 分析数据集时，可以将列分为Categorical和Numerical。 为了一起分析它们，应将一种类型转换为另一种类型。 如果您具有分类值，并且想要将其转换为数字值，则可以通过将每个类别转换为数字值来实现。 该过程称为编码。 另一方面，可以通过将数值转换为类别对应物来执行装箱，类别对应物是通过将数据分成间隔来创建的。
我们将通过逐一探讨其特征来开始分析。 在统计中，此方法称为单变量分析。 单变量分析的目的主要围绕描述。 我们将计算最小值，最大值，范围，百分位数，均值和方差，然后绘制一些直方图并分析每个特征的分布。 我们将介绍偏度和峰度的概念，然后看一下修剪的重要性。 完成单变量分析后，我们将继续进行双变量分析，这意味着同时分析两个特征。 为此，我们将探索两组特征之间的关系：
In [39]: np.set_printoptions(suppress=True, linewidth=125)
minimums = np.round(np.amin(samples, axis=0), decimals=1)
maximums = np.round(np.amax(samples, axis=0), decimals=1)
range_column = np.round(np.ptp(samples, axis=0), decimals=1)
mean = np.round(np.mean(samples, axis=0), decimals=1)
median = np.round(np.median(samples, axis=0), decimals=1)
variance = np.round(np.var(samples, axis=0), decimals=1)
tenth_percentile = np.round(np.percentile(samples, 10, axis=0), decimals = 1)
ninety_percentile = np.round(np.percentile(samples, 90 ,axis=0), decimals = 1)
In [40]: range_column
Out[40]: array([ 89\. , 100\. , 27.3, 1\. , 0.5, 5.2, 97.1, 11\. , 23\. , 524\. , 9.4, 396.6, 36.2])
在前面的代码中，我们首先使用set_printoptions()方法设置打印选项，以便查看四舍五入的小数位数，并且其线宽足以容纳所有列。 要计算基本统计信息，我们使用numpy函数，例如amin()，amax()，mean()，median()，var()，percentile()和ptp()。 所有计算都是按列进行的，因为每一列都代表一个要素。 最终的数组看起来有点草率且无用，因为您看不到哪行显示了哪些统计信息：
为了打印对齐的 numpy 数组，可以使用zip()函数添加行标签并在数组之前打印列标签。 在 SciPy 中，您可以使用更多统计函数来计算基本统计信息。 SciPy 提供describe()函数，该函数返回给定数组的一些描述性统计信息。 您可以使用一个函数来计算点，最小值，最大值，均值，方差，偏度和峰度，然后分别解析它们，如下所示：
以下代码块分别计算基本统计信息并将其堆叠到最终数组中：
In [45]: minimum = arr.minmax[0]
maximum = arr.minmax[1]
mean = arr.mean
median = np.round(np.median(samples,axis = 0), decimals = 1)
variance = arr.variance
tenth_percentile = stats.scoreatpercentile(samples, per = 10, axis = 0)
ninety_percentile = stats.scoreatpercentile(samples, per =90, axis = 0)
rng = stats.iqr(samples, rng = (20,80), axis = 0)
np.set_printoptions(suppress = True, linewidth = 125)
Basic_Statistics1 = np.round(np.vstack((minimum,maximum,rng, mean,median,variance,tenth_percentile,ninety_percentile)), Basic_Statistics1.shape
Out[45]: (8, 13)
In [46]: stat_labels1 = ['minm', 'maxm', 'rang', 'mean', 'medi', 'vari', '50%t', '90%t']
与 NumPy 相比，我们可以使用iqr()函数来计算范围。 此函数计算沿指定轴和范围（rng参数）的数据的四分位数范围。 默认情况下，rng = (25, 75)，这意味着该函数将计算数据的第 75 和第 25 个百分位数之间的差。 为了返回与numpy.ptp()函数相同的输出，可以使用rng = (0, 100)，因为这将计算所有给定数据的范围。 我们使用stat.scoreatpercentile()作为与numpy.percentile()方法的等效物，来计算特征的第 10 个和第 90 个百分位值。 如果您看一眼结果，您会发现几乎每个特征的差异都很大。 您可以看到，由于我们通过将参数传递为(20,80)来限制范围计算，因此范围值显着减小。 它也向您展示了我们在分布式要素的两侧都有许多极值。 从我们的结果可以得出结论，对于大多数特征，平均值均高于中位数，这表明我们这些特征的分布偏向右侧。 在下一节中，当我们绘制直方图然后深入分析这些特征的偏度和峰度时，您将清楚地看到这一点。
计算直方图
直方图是数字数据分布的直观表示。 一百多年前，卡尔·皮尔森（Karl Pearson）首次提出了这一概念。 直方图是一种用于连续数据的条形图，而条形图直观地表示类别变量。 第一步，您需要将整个值范围划分为一系列间隔（仓位）。 每个桶必须相邻，并且它们不能重叠。 通常，桶大小相等，要包含的桶数量的经验法则是 5–20。 这意味着，如果您有 20 个以上的容器，则图形将很难阅读。 相反，如果您的箱数少于 5 个，则您的图形将很少了解数据的分布：
In [48]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
NOX = samples[:,5:6]
plt.hist(NOX,bins ='auto')
plt.title("Distribution nitric oxides concentration (parts per 10 million)")
plt.show()
在前面的代码中，我们绘制了特征NOX的直方图。 箱柜计算自动完成，如下所示：
您可以通过将切片数组作为第一个参数，使用pyplot.hist()绘制直方图。 y 轴显示数据集中每个间隔（桶）中有多少值的计数， x 轴表示其值。 通过设置normed=True，您可以看到箱柜的百分比，如下所示：
In [49]: plt.hist(NOX,bins ='auto', normed = True)
plt.title("Distribution nitric oxides concentration (parts per 10 million)")
plt.show()
当您查看直方图时，可能很难计算每个桶及其边界的大小。 当pyplot.hist()返回包含此信息的元组时，您可以轻松地获得这些值，如下所示：
In [50]: import matplotlib.pyplot as plt
NOX = samples[:,5:6]
n, bins, patches = plt.hist(NOX, bins='auto')
print('Bin Sizes')
print(n)
print('Bin Edges')
print(bins)
在前面的代码中，我们打印出每个桶中包含多少个数字，以及桶的边界是什么，如下所示：
让我们解释一下前面的输出。 第一个桶的大小为1，这意味着此桶中仅存在一个值下降。 第一仓的间隔在3.561和3.74096552之间。 为了使它们更整洁，可以使用以下代码块，该代码块有意义地堆叠了这两个数组（bin_interval，bin_size）：
In [51]: bins_string = bins.astype(np.str)
n_string = n.astype(np.str)
lists = []
for i in range(0, len(bins_string)-1):
c = bins_string[i]+ "-" + bins_string[i+1]
lists.append(c)
new_bins = np.asarray(lists)
Stacked_Bins = np.vstack((new_bins, n_string)).T
Stacked_Bins
确定桶的数量及其宽度非常重要。 一些理论家对此进行试验以确定最佳拟合。 下表显示了最受欢迎的估计器。 您可以将估计器设置为numpy.histogram()中的bins参数，以相应地更改桶计算。 这些方法由pyplot.hist()函数隐式支持，因为其参数传递给numpy.histogram()：
| 估计器 | 公式 |
|---|---|
| Freedman-Diaconis 估计 |  |
| Doane 公式 |  |
| Rice 法则 |  |
| Scott 的常规参考法则 |  |
| Sturges 公式 |  |
| 平方根选项 |  |
IQR = 四分位间距
g[1] = 分布的估计三阶偏度
所有这些方法都有不同的优势。 例如，Strurges 公式对于高斯数据是最佳的。 Rice 的规则是 Sturges 公式的简化版本，并假定近似正态分布，因此，如果数据不是正态分布，则可能效果不佳。 Doane 的公式是 Sturges 公式的改进版本，尤其是在具有非正态分布的情况下。 Freedman–Diaconis 估计器是 Scott 规则的修改版本，他用IQR=2代替了 3.5 个标准差，这使公式对异常值的敏感性降低。 平方根选择是一种非常常见的方法，许多工具都以其快速简便的方式使用它。 在numpy.histogram()中，还有一个名为'auto'的选项，它为我们提供了最大的 Sturges 和 Freedman–Diaconis 估计器方法。 让我们看看这些方法将如何影响直方图：
In [52]: fig, ((ax1, ax2, ax3),(ax4,ax5,ax6)) = plt.subplots(2,3,sharex=True)
axs = [ax1,ax2,ax3,ax4,ax5,ax6]
list_methods = ['fd','doane', 'scott', 'rice', 'sturges','sqrt']
plt.tight_layout(pad=1.1, w_pad=0.8, h_pad=1.0)
for n in range(0, len(axs)):
axs[n].hist(NOX,bins = list_methods[n])
axs[n].set_title('{}'.format(list_methods[n]))
在前面的代码中，我们编译了六个直方图，它们全部共享相同的 x 轴，如下所示：
所有直方图都表示相同的特征，但具有不同的合并方法。 例如，fd方法显示数据看起来接近其正态分布，而另一方面，doane方法看起来更像学生的 T 分布。 另外，sturges方法会创建两个小容器，因此很难分析数据。 当我们查看基本统计数据时，我们指出大多数特征的均值均高于其中位数，因此他们将数据向右倾斜。 但是，如果您查看sturges，rice和sqrt方法，则很难这么说。 但是，当我们使用自动箱绘制直方图时，这是显而易见的：
In [53]: import numpy as np
samples_new = np.delete(samples, 3 , axis=1)
samples_new.shape
Out[53]: (506, 12)
In [54]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
fig,((ax1, ax2 , ax3),(ax4, ax5, ax6), (ax7, ax8, ax9), (ax10, ax11, ax12)) = plt.subplots(4,3, figsize = (10,15))
axs =[ax1, ax2 , ax3, ax4, ax5, ax6, ax7, ax8, ax9, ax10, ax11, ax12]
feature_names_new = np.delete(feature_names,3)
for n in range(0, len(axs)):
axs[n].hist(samples_new[:,n:n+1], bins ='auto', normed = True)
axs[n].set_title('{}'.format(feature_names[n]))
在前面的代码中，我们在单个布局中编译所有特征直方图。 这将帮助我们轻松比较它们。 前面代码的输出如下：
在前面的代码中，我们删除了CHAS，因为它包含二进制值，这意味着直方图将无助于我们进一步了解此特征。 我们还从特征列表中取出了特征名称，以便正确绘制其余特征。
从这些图表我们可以得出结论，在大多数城镇中，人均犯罪率非常低，但是在某些城镇中，这一比例非常高。 通常，住宅用地面积低于 25,000 英尺。在许多情况下，非零售业务用地所占土地不到整个城镇的 10% 。 另一方面，我们可以看到许多城镇的非零售业务用地约为 20% 。 一氧化氮的浓度非常偏斜，尽管此处存在一些离平均值很远的异常值。 每个住宅的平均房间数在 5 至 7 间之间。 在 1940 年之前建造的架构物中，超过 50% 的架构物被其所有者占用。 他们中的大多数人离波士顿的就业中心都不远。 超过 10% 的人无法轻松到达径向公路。 对于相当多的人来说，财产税很高。 通常，教室规模为 15 至 20 个孩子。 大多数城镇中居住的黑人比例非常相似。 大多数城镇居民的经济地位较低。 通过查看这些直方图，您可以得出许多其他结论。 如您所见，直方图在描述数据的分布方式时非常有用，并且可以看到平均值，方差和离群值。 在下一节中，我们将重点介绍偏度和峰度。
解释偏度和峰度
在统计分析中，力矩是一种定量度量，用于描述距参考点的预期距离。 如果期望参考点，则将其称为中心矩。 在统计中，中心矩是与均值相关的矩。 第一个和第二个矩分别是平均值和方差。 平均值是您的数据点的平均值。 方差是每个数据点与平均值的总偏差。 换句话说，方差表明您的数据如何与均值分散。 第三个中心矩是偏度，它测量均值分布的不对称性。 在标准正态分布中，偏度是对称的，因此等于零。 另一方面，如果均值
In [55]: %matplotlib notebook
%matplotlib notebook
from scipy.stats import skewnorm
fig, (ax1, ax2, ax3) = plt.subplots(1,3 ,figsize=(10,2.5))
x1 = np.linspace(skewnorm.ppf(0.01,-3), skewnorm.ppf(0.99,-3),100)
x2 = np.linspace(skewnorm.ppf(0.01,0), skewnorm.ppf(0.99,0),100)
x3 = np.linspace(skewnorm.ppf(0.01,3), skewnorm.ppf(0.99,3),100)
ax1.plot(skewnorm(-3).pdf(x1),'k-', lw=4)
ax2.plot(skewnorm(0).pdf(x2),'k-', lw=4)
ax3.plot(skewnorm(3).pdf(x3),'k-', lw=4)
ax1.set_title('Left Skew')
ax2.set_title('Normal Dist')
ax3.set_title('Right Skew')
您可以使用skewnorm()函数轻松创建具有偏斜度的正态分布。 如您在前面的代码中看到的，我们生成一个具有 100 个值的百分点函数（与累积分布函数百分位数相反），然后创建具有不同偏斜度的线。 您不能直接得出结论，左偏斜比右偏斜好，反之亦然。 当分析数据的偏斜度时，您需要考虑这些尾巴会导致什么。 例如，如果您正在分析股票收益的时间序列，并将其绘制后看到一个正确的偏斜，则不必担心获得比预期更高的收益，因为这些离群值不会给您的股票带来任何风险。 交易策略。 当您分析服务器的响应时间时，可能会出现类似的示例。 当您绘制响应时间的概率密度函数时，您对左尾巴并不真正感兴趣，因为它们代表了快速的响应时间。
关于均值的第四个中心矩是峰度。 峰度描述了与分布曲线的平坦程度或变薄程度有关的尾部（离群值）。 在正态分布中，峰度等于 3。峰度主要有三种类型：Leptokurtic（薄），Mesokurtic和Platykurtic（平坦）：
In [56]: %matplotlib notebook
%matplotlib notebook
import scipy
from scipy import stats
import matplotlib.pyplot as plt
fig, (ax1, ax2, ax3) = plt.subplots(1, 3 , figsize=(10,2))
axs= [ax1, ax2, ax3]
Titles = ['Mesokurtic', 'Lebtokurtic', 'Platykurtic']
#Mesokurtic Distribution - Normal Distribution
dist = scipy.stats.norm(loc=100, scale=5)
sample_norm = dist.rvs(size = 10000)
#leptokurtic Distribution
dist2 = scipy.stats.laplace(loc= 100, scale= 5)
sample_laplace = dist2.rvs(size= 10000)
#platykurtic Distribution
dist3 = scipy.stats.cosine(loc= 100, scale= 5)
sample_cosine = dist3.rvs(size= 10000)
samples = [sample_norm, sample_laplace, sample_cosine]
for n in range(0, len(axs)):
axs[n].hist(samples[n],bins= 'auto', normed= True)
axs[n].set_title('{}'.format(Titles[n]))
print ("kurtosis of" + Titles[n])
print(scipy .stats.describe(samples[n])[5])
在前面的代码中，您可以更清楚地看到分布形状的差异。 所有峰度值均使用stats.describe()方法进行归一化，因此正态分布的实际峰度值约为 3，而归一化版本为 0.02。 让我们快速检查一下我们特征的偏度和峰度值是什么：
In [57]: samples,label, feature_names = dataset.data , dataset.target , dataset.feature_names
for n in range(0, len(feature_names_new)):
kurt = scipy.stats.describe(samples[n])[5]
skew = scipy.stats.describe(samples[n])[4]
print (feature_names_new[n] + "-Kurtosis: {} Skewness: {}" .format(kurt, skew))
CRIM-Kurtosis: 2.102090573040533 Skewness: 1.9534138515494224
ZN-Kurtosis: 2.8706349006925134 Skewness: 2.0753333576721893
INDUS-Kurtosis: 2.9386308786131767 Skewness: 2.1061627843164086
NOX-Kurtosis: 3.47131446484547 Skewness: 2.2172838215060517
RM-Kurtosis: 3.461596258869246 Skewness: 2.2086627738768234
AGE-Kurtosis: 3.395079726813977 Skewness: 2.1917520072643533
DIS-Kurtosis: 1.9313625761956317 Skewness: 1.924572804475305
RAD-Kurtosis: 1.7633603556547106 Skewness: 1.8601991629604233
TAX-Kurtosis: 1.637076772210217 Skewness: 1.8266096199819994
PTRATIO-Kurtosis: 1.7459544645159752 Skewness: 1.8679592455694167
B-Kurtosis: 1.7375702020429316 Skewness: 1.8566444885400044
LSTAT-Kurtosis: 1.8522036606250456 Skewness: 1.892802610207445
从结果中可以看出，所有特征都具有正偏度，这表明它们向右偏斜。 就峰度而言，它们都具有正值，尤其是 NOX 和 RM，它们具有非常高的色氨酸。 我们可以得出结论，它们全都是瘦腿型的，这表明数据更加集中在均值上。 在下一节中，我们将着重于数据修整并使用修整后的数据计算统计数据。
整理统计量
正如您在上一节中已经注意到的那样，我们特征的分布非常分散。 处理模型中的异常值是分析中非常重要的部分。 当您查看描述性统计数据时，它也非常重要。 由于这些极端值，您很容易混淆并曲解分布。 SciPy 具有非常广泛的统计函数，可以计算有关修剪数据的描述性统计数据。 使用调整后的统计数据的主要思想是删除异常值（尾部），以减少其对统计计算的影响。 让我们看看如何使用这些函数，以及它们如何影响我们的特征分布：
In [58]: np.set_printoptions(suppress= True, linewidth= 125)
samples = dataset.data
CRIM = samples[:,0:1]
minimum = np.round(np.amin(CRIM), decimals=1)
maximum = np.round(np.amax(CRIM), decimals=1)
variance = np.round(np.var(CRIM), decimals=1)
mean = np.round(np.mean(CRIM), decimals=1)
Before_Trim = np.vstack((minimum, maximum, variance, mean))
minimum_trim = stats.tmin(CRIM, 1)
maximum_trim = stats.tmax(CRIM, 40)
variance_trim = stats.tvar(CRIM, (1,40))
mean_trim = stats.tmean(CRIM, (1,40))
After_Trim = np.round(np.vstack((minimum_trim,maximum_trim,variance_trim,mean_trim)), decimals=1)
stat_labels1 = ['minm', 'maxm', 'vari', 'mean']
Basic_Statistics1 = np.hstack((Before_Trim, After_Trim))
print (" Before After")
for stat_labels1, row1 in zip(stat_labels1, Basic_Statistics1):
print ('%s [%s]' % (stat_labels1, ''.join('%07s' % a for a in row1)))
Before After
minm [ 0.0 1.0]
maxm [ 89.0 38.4]
vari [ 73.8 48.1]
mean [ 3.6 8.3]
要计算调整后的统计信息，我们使用了tmin()，tmax()，tvar()和tmean()，如前面的代码所示。 所有这些函数都将极限值作为第二个参数。 在CRIM函数中，我们可以在右侧看到很多尾巴，因此我们将数据限制为(1, 40)，然后计算统计信息。 在修整值之前和之后，您可以通过比较计算的统计信息来观察差异。 作为tmean()的替代方法，如果要从两侧按比例分割数据，可以使用trim_mean()函数。 您可以看到我们的均值和方差在修剪数据后如何显着变化。 当我们删除了 40 至 89 之间的许多极端离群值时，方差显着减小。相同的修整对均值有不同的影响，此后均值增加了一倍以上。 原因是先前的分布中有很多零，并且通过将计算限制在 1 和 40 的值之间，我们删除了所有这些零，这导致了更高的平均值。 请注意，上述所有函数只是在计算这些值时即时修剪数据，因此CRIM数组不会被修剪。 如果要从两侧修剪数据，则可以在一侧使用trimboth()和trim1()。 在这两种方法中，应该使用比例作为参数，而不是使用极限值。 例如，如果您传递proportiontocut =0.2，它将把您最左边和最右边的值切成 20% ：
In [59]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
CRIM_TRIMMED = stats.trimboth(CRIM, 0.2)
fig, (ax1, ax2) = plt.subplots(1,2 , figsize =(10,2))
axs = [ax1, ax2]
df = [CRIM, CRIM_TRIMMED]
list_methods = ['Before Trim', 'After Trim']
for n in range(0, len(axs)):
axs[n].hist(df[n], bins = 'auto')
axs[n].set_title('{}'.format(list_methods[n]))
从两侧修剪掉 20% 的值后，我们可以更好地解释分布，实际上可以看到大多数值在 0 到 1.5 之间。 仅查看左方的直方图就很难获得这种见解，因为极值主导了直方图，并且我们只能看到 0 周围的一条线。因此，trimmed函数在探索性数据分析中非常有用。 在下一部分中，我们将重点讨论箱形图，它们对于数据的描述性分析和异常值检测非常有用且非常流行。
箱形图
探索性数据分析中的另一个重要视觉效果是箱形图，也称为箱须图。 它基于五位数摘要建立，该摘要是最小值，第一四分位数，中位数，第三四分位数和最大值。 在标准箱形图中，这些值表示如下：
这是比较几种分布的一种非常方便的方法。 通常，情节的胡须通常延伸到极限点。 或者，您可以将其切入 1.5 个四分位范围。 让我们检查一下CRIM和RM函数：
In [60]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
from scipy import stats
samples = dataset.data
fig, (ax1,ax2) = plt.subplots(1,2, figsize =(8,3))
axs = [ax1, ax2]
list_features = ['CRIM', 'RM']
ax1.boxplot(stats.trimboth(samples[:,0:1],0.2))
ax1.set_title('{}'.format(list_features[0]))
ax2.boxplot(stats.trimboth(samples[:,5:6],0.2))
ax2.set_title('{}'.format(list_features[1]))
正如您在前面的代码中看到的那样，RM 值一直沿最小值和最大值连续分散，以便胡须看起来像一条线。 您可以轻松地检测到中位数几乎位于第一和第三四分位数之间的中间位置。 在CRIM中，极值不会沿着最小值继续； 它们更像是单个外围数据点，因此表示就像一个圆圈。 您可以看到这些离群值大多位于第三个四分位数之后，并且中位数非常接近第一个四分位数。 由此，您还可以得出结论，分布是右偏的。 因此，箱形图是直方图的非常有用的替代方法，因为您可以轻松分析分布并一次比较多个分布。 在下一节中，我们将继续进行双变量分析，并研究标签数据与要素的相关性。
计算相关性
本节专门用于双变量分析，您可以在其中分析两列。 在这种情况下，我们通常研究这两个变量之间的关联，这称为相关性。 相关性显示两个变量之间的关系，并回答一些问题，例如，如果变量 B 增加 10% ，变量 A 将会发生什么？ 在本节中，我们将说明如何计算数据的相关性并以二维散点图表示它。
通常，相关是指任何统计依赖性。 相关系数是计算相关度量的定量值。 您可以将相关性和相关系数之间的关系视为湿度计和湿度之间的相似关系。 相关系数最流行的类型之一是皮尔逊积矩相关系数。 相关系数的值在 +1 和 -1 之间，其中 -1 表示两个变量之间的强负线性关系，而 +1 表示两个变量之间的强正线性关系。 在 NumPy 中，我们可以使用corrcoef()方法来计算系数相关矩阵，如下所示：
In [61]: np.set_printoptions(suppress= True, linewidth = 125)
CorrelationCoef_Matrix = np.round(np.corrcoef(samples, rowvar= False), decimals= 1)
CorrelationCoef_Matrix
Seaborn 是一个基于 matplotlib 的统计数据可视化库，您可以使用它创建非常吸引人的美观统计图形。 它是一个非常受欢迎的库，具有精美的可视化效果，并且与流行的软件包（尤其是 Pandas）具有完美的兼容性。 您可以使用seaborn程序包中的热图来形象化相关系数矩阵。 当您具有数百种特征时，它对于检测高相关系数非常有用：
In [62]: CorrelationCoef_Matrix1 = np.round(np.corrcoef(samples, rowvar= False), decimals= 1)
CorrelationCoef_Matrix1
import seaborn as sns; sns.set()
ax = sns.heatmap(CorrelationCoef_Matrix1, cmap= "YlGnBu")
在这里，我们具有标签列特征的相关系数。 您可以在corrcoef()函数中添加其他变量集作为第二个参数，就像我们对label列所做的那样，如下面的屏幕快照所示。 只要形状相同，该函数将在最后堆叠此列并计算相关系数矩阵：
如您所见，除了F13之外，大多数特征都具有弱或中等的负线性关系。 另一方面，F6具有很强的正线性关系。 让我们绘制此特征并用散点图查看这种关系。 以下代码块借助matplotlib显示了三个不同的特征散点图（RM，DIS和LSTAT）和一个标签列：
In [64]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
from scipy import stats
fig, (ax1, ax2, ax3) = plt.subplots(1,3 ,figsize= (10,4))
axs =[ax1,ax2,ax3]
feature_list = [samples[:,5:6], samples[:,7:8], samples[:,12:13]]
feature_names = ["RM", "DIS", "LSTAT"]
for n in range(0, len(feature_list)):
axs[n].scatter(label, feature_list[n], edgecolors=(0, 0, 0))
axs[n].set_ylabel(feature_names[n])
axs[n].set_xlabel('label')
在前面的代码中，RM 和标签值分散为一条正线性线，这与相关系数值为 0.7 一致，如前面的屏幕快照所示。 散点图表明样品的 RM 值越高，标签的值越高。 在中间的散点图中，您可以看到系数相关性等于 0.25 的位置。 这显示出非常弱的正线性相关性，这在散点图中也可以看到。 我们可以得出结论，由于值分散在各处，因此没有明确的关系。 第三个散点图显示了强大的线性关系，相关系数为 -0.7。 随着LSTAT下降，标签值开始增加。 所有这些相关矩阵和散点图均使用未修剪的数据对此进行了计算。 让我们看看如何针对每个要素和标签从两侧修剪数据 10% 来改变数据集的线性关系结果：
In [65]: %matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
from scipy import stats
fig, (ax1, ax2, ax3) = plt.subplots(1,3 ,figsize= (9,4))
axs = [ax1, ax2, ax3]
RM_tr = stats.trimboth(samples[:,5:6],0.1)
label_tr = stats.trimboth(label, 0.1)
LSTAT_tr = stats.trimboth(samples[:,12:13],0.1)
DIS_tr = stats.trimboth(samples[:,7:8],0.1)
feature_names = ["RM", "DIS", "LSTAT"]
feature_list = [RM_tr, DIS_tr, LSTAT_tr]
for n in range(0, len(feature_list)):
axs[n].scatter(label_tr,feature_list[n], edgecolors=(0, 0, 0))
axs[n].set_ylabel(feature_names[n])
axs[n].set_xlabel('label')
修剪数据后，您可以看到所有特征都与标签具有很强的正线性相关性，尤其是DIS和LSTAT的相关性，其中标签发生了巨大变化。 这显示了修剪的力量。 如果您不知道如何处理异常值，则很容易误解您的数据。 离群值可以更改分布的形状以及其他变量之间的相关性，最终它们可以影响模型的性能。
总结
在本章中，我们介绍了使用 NumPy，SciPy，matplotlib 和 Seaborn 软件包进行探索性数据分析的过程。 首先，我们学习了如何加载和保存文件以及探索数据集。 然后，我们解释并计算了重要的统计中心矩，例如均值，方差，偏度和峰度。 四个重要的可视化分别用于单变量和变量分析的图形表示; 这些是直方图，箱形图，散点图和热图。 还通过示例强调了数据修整的重要性。
在下一章中，我们将更进一步，并开始使用线性回归预测房价。
四、使用线性回归预测房价
在本章中，我们将通过实现线性回归来介绍监督学习和预测建模。 在上一章中，您了解了探索性分析的知识，但尚未研究建模。 在本章中，我们将创建一个线性回归模型来预测房地产市场价格。 广义上讲，我们将借助目标变量与其他变量的关系来预测目标变量。 线性回归被广泛使用，并且是监督机器学习算法的简单模型。 本质上是关于为观察到的数据拟合一条线。 我们将从解释监督学习和线性回归开始我们的旅程。 然后，我们将分析线性回归的关键概念，例如自变量和因变量，超参数，损失和误差函数以及随机梯度下降。 对于建模，我们将使用与上一章相同的数据集。
本章将介绍以下主题：
- 监督学习和线性回归
- 自变量和因变量
- 超参数
- 损失和误差函数
- 为单个变量实现我们的算法
- 计算随机梯度下降
- 使用线性回归建模房价
监督学习和线性回归
机器学习使计算机系统无需显式编程即可学习。 监督学习是最常见的机器学习类型之一。 监督学习由一组不同的算法组成，这些算法提出学习问题并通过使用历史数据映射输入和输出来解决它们。 该算法分析输入和相应的输出，然后将它们链接在一起以找到关系（学习）。 最后，对于新的给定数据集，它将使用此学习来预测输出。
为了区分监督学习和非监督学习，我们可以考虑基于输入/输出的建模。 在监督学习中，将对计算机系统的每组输入数据使用标签进行监督。 在无监督学习中，计算机系统将仅使用没有任何标签的输入数据。
例如，假设我们有 100 万张猫和狗的照片。 在监督学习中，我们标记输入数据并声明给定的照片是猫还是狗。 假设每张照片（输入数据）有 20 个特征。 标有照片的计算机系统将知道照片是猫还是狗（输出数据）。 当我们向计算机系统显示一张新照片时，它将通过分析新照片的 20 个特征来确定它是猫还是狗，并根据其先前的学习进行预测。 在无监督学习中，我们将只拥有 100 万张猫和狗的照片，而没有任何标签说明是猫还是狗的照片，因此该算法将在没有我们监督的情况下通过分析其特征来对数据进行聚类。 聚类完成后，会将新照片输入到无监督学习算法中，系统会告诉我们该照片属于哪个聚类。
在这两种情况下，系统都将具有简单或复杂的决策算法。 唯一的区别是是否有任何初始监督。 监督学习方法的概述方案如下：
监督学习可以分为分类和回归两种类型，如上图所示。 分类模型预测标签。 例如，可以将前面的示例视为监督分类问题。 为了执行分类，我们需要训练分类算法，例如支持向量分类器（SVC），随机森林， K 近邻（KNN），依此类推。 基本上，分类器是指用于对数据进行分类（分类）的算法。
当我们的目标变量是分类变量时使用分类方法，而当我们的目标变量是连续变量时，则应用回归模型，因为目标是预测数值而不是类别。
考虑一下我们在上一章中使用的数据集：波士顿房屋价格数据集。 在该数据集中，我们的目的是对特征值进行统计分析，因为我们需要了解特征值的分布方式，其基本统计量以及之间的相互关系。 最后，我们想知道每个特征如何导致房价上涨。 它对正面，负面还是根本没有影响？ 如果特征 x 与房屋价格 A 之间存在潜在的影响（关系），则该关系的强弱是多少？
我们试图通过建立模型来预测具有给定特征的房价来回答这些问题。 结果，在为模型提供新特征后，我们期望模型将输出变量生成为连续值（150k，120，154 美元等）。 让我们看一个非常基本的工作流程：
如上图所示，分析从为模型准备数据开始。 在此阶段，我们应该清理数据，处理缺失的值，并提取将要使用的特征。 数据干净后，我们应该将其分为训练数据和测试数据两部分，以测试模型的性能。
模型验证中的重要部分是过拟合的概念。 用外行的术语来说，过度拟合意味着从训练数据集中学到太多，因此我们的模型过度拟合并为训练数据集产生近乎完美的结果。 但是，对于从未见过的数据，它的灵活性不足以产生良好的结果，因此无法很好地泛化。
可以将数据集分为训练，验证（强烈推荐）和测试数据集来解决此问题。 训练数据是算法最初学习算法参数（权重）并在其中将误差最小化的地方建立模型的地方。 当您有几种算法并且需要调整超参数时，或者当算法中有许多参数并且需要调整参数时，平移数据集非常有用。 测试数据集用于性能评估。
简而言之，您可以使用训练数据来训练算法，然后在验证数据集中微调算法的参数或权重，最后一步，在测试数据集中测试调整后的算法的性能。
与过度拟合相反的情况是欠拟合，这意味着该算法从数据中学习得更少，并且我们的算法与我们的观察结果不太吻合。 让我们以图形方式查看过度拟合，欠拟合和最佳拟合的样子：
您可能已经注意到，在上图中，即使过度拟合看起来非常合适，它也会生成一个回归线，该回归线对于该数据集非常独特，并且无法正确捕获特征。 第二个图，即欠拟合图，实际上无法捕获数据的形状。 当我们的数据是非线性的时，它没有从中学习并产生了线性回归线。 第三个图是我们的最佳拟合图，拟合得很好，掌握了分布特征并产生了一条曲线。 我们可以预期，对于此数据集的延续，它不会有令人失望的性能指标。
在本章中，我们将使用线性回归作为有监督的学习方法。 我们将首先解释非常重要的概念，例如自变量和因变量，超参数以及损失和误差函数。 我们还将介绍线性回归的实际示例。 在下一章中，我们将介绍线性回归模型的最重要组成部分：独立变量和因变量。
自变量和因变量
正如我们在前面的小节中提到的，线性回归用于基于其他变量来预测变量的值。 我们正在研究输入变量X与输出变量Y之间的关系。
在线性回归中，因变量是我们要预测的变量。 之所以称其为因变量，是因为线性回归背后的假设。 该模型假设这些变量取决于方程另一侧的变量，这些变量称为独立变量。
在简单回归模型中，模型将解释因变量如何基于自变量而变化。
例如，假设我们要分析基于给定产品的价格变化如何影响销售值。 如果您仔细阅读此句子，则可以轻松检测出我们的因变量和自变量。 在我们的示例中，我们假设销售值受价格变化的影响，换句话说，销售值取决于产品的价格。 结果，销售值是从属值，价格是独立的值。 这不一定意味着给定产品的价格不依赖于任何东西。 当然，它取决于许多因素（变量），但是在我们的模型中，我们假设价格是给定的，并且给定的价格会改变销售值。 线性回归线的公式如下：
哪里：
Yi = 估计值或因变量
B0 = 截距
B1 = 斜率
Xi = 自变量或探索变量：
线性回归线的斜率（B1）实际上显示了这两个变量之间的关系。 假设我们将斜率计算为 0.8。 这意味着自变量增加 1 个单位可能会使估计值增加 0.8 个单位。 前面的线性回归线仅生成估计，这意味着它们仅是给定X的Y的预测。 如下图所示，每个观测值和线性线之间有一段距离。 该距离称为误差，这是预期的，在回归线拟合和模型评估中也非常重要：
拟合线性回归线的最常见方法是使用最小二乘方法。 该方法通过最小化这些误差的平方和来拟合回归线。 计算公式如下：
这些误差均方根的原因是我们不希望负误差和正误差彼此抵消。 在模型评估中，使用 R 平方，F 检验和均方根误差（RMSE）。 他们都使用平方和（SST）和平方和误差（SSE）作为基本度量。 如您所见，在 SSE 中，我们再次计算预测值和实际值之间的差，取平方并求和，以评估回归线对数据的拟合程度。
如前所述，最小二乘方法旨在最小化平方误差（残差），并找到最适合数据点的斜率和截距值。 由于这是一种封闭形式的解决方案，因此您可以轻松地手动计算它，以便在后台查看此方法的作用。 让我们使用一个包含少量数据集的示例：
| 牛奶消耗量（每周升） | 身高 |
|---|---|
| 14 | 175 |
| 20 | 182 |
| 10 | 170 |
| 15 | 185 |
| 12 | 164 |
| 15 | 173 |
| 22 | 181 |
| 25 | 193 |
| 12 | 160 |
| 13 | 165 |
假设我们有 10 个观察值，如上表所示，用于每周食用牛奶和食用牛奶的人的身高值。 如果我们绘制这些数据，我们可以看到每日牛奶消耗量与身高之间存在正相关：
现在，我们要使用最小二乘法拟合线性回归线，该方法通过对斜率和截距使用以下公式来完成：
然后，我们需要创建一个表格来帮助我们进行计算，例如x，y，xy，x^2和y^2的总和：
x（牛奶消耗量） |
y（高度） |
xy |
x^2 |
y^2 |
|
|---|---|---|---|---|---|
| 1 | 14 | 175 | 2,450 | 196 | 30,625 |
| 2 | 20 | 182 | 3,640 | 400 | 33,124 |
| 3 | 10 | 170 | 1,700 | 100 | 28,900 |
| 4 | 15 | 185 | 2,775 | 225 | 34,225 |
| 5 | 12 | 164 | 1,968 | 144 | 26,896 |
| 6 | 15 | 173 | 2,595 | 225 | 29,929 |
| 7 | 22 | 181 | 3,982 | 484 | 32,761 |
| 8 | 25 | 193 | 4,800 | 625 | 37,249 |
| 9 | 12 | 160 | 1,920 | 144 | 25,600 |
| 10 | 13 | 165 | 2,145 | 169 | 27,225 |
| ∑ | 158 | 1,748 | 27,975 | 2,712 | 306,534 |
B0 = (1748*2712) - (158*27975) / (10*2712) - (158)^2 = 320526/2156 = 148.66
B1 = (10*27975) - (158*1748) / (10*2712) - (158)^2 = 1.65
然后，我们可以将回归线公式如下：
在本小节中，我们提到了自变量和因变量，并介绍了线性回归线和拟合方法。 在下一节中，我们将介绍超参数，这些参数在回归模型调整中非常有用。
超参数
在开始之前，也许最好解释一下为什么我们称它们为超参数而不是参数。 在机器学习中，可以从数据中学习模型参数，这意味着在训练模型时，您可以拟合模型的参数。 另一方面，我们通常在开始训练模型之前先设置超参数。 为了举例说明，您可以将回归模型中的系数视为模型参数。 以超参数为例，我们可以说许多不同模型中的学习率或 K 均值聚类中的聚类数（k）。
另一个重要的事情是模型参数和超参数之间的关系，以及它们如何塑造我们的机器学习模型，换句话说，就是我们模型的假设。 在机器学习中，参数用于配置模型，此配置将为我们的数据集专门定制算法。 我们需要处理的是如何优化超参数。 另外，如上所述，可以在验证期间执行该优化。 在许多情况下，优化超参数将带来性能提升。
您还可以将超参数视为模型参数之上的高级参数。 想象一下您使用无监督学习的 K 均值聚类的情况。 如果您使用误差的群集号（K）作为超参数，则可以确保您的数据不适合。
现在，您应该问的是，如果在训练模型之前手动设置超参数，我们该如何调整它们。 有几种方法可以调整超参数。 此优化的底线是使用一组不同的超参数测试算法，并在选择性能更好的超参数集之前针对每种情况计算误差函数或损失函数。
在本节中，我们简要介绍了参数，超参数及其差异。 在下一节中，我们将介绍损失和误差函数，这对于超参数优化非常重要。
损失和误差函数
在前面的小节中，我们解释了有监督和无监督的学习。 无论使用哪种机器学习算法，我们的主要挑战都是关于优化的问题。 在优化函数中，我们实际上是在尝试使损失函数最小化。 设想一下您试图优化每月储蓄的情况。 在关闭状态下，您要做的就是最小化支出，换句话说，将损失函数最小化。
建立损失函数的一种非常常见的方法是从预测值与实际值之差开始。 通常，我们尝试估计模型的参数，然后进行预测。 我们可以用来评估我们的预测水平的主要度量包括计算实际值之间的差：
在不同的模型中，使用了不同的损失函数。 例如，您可以在回归模型中使用均方误差，但是将其用作分类模型的损失函数可能不是一个好主意。 例如，您可以按以下方式计算均方误差：
其中回归模型如下：
有许多不同的损失函数可用于不同的机器学习模型。 以下是一些重要的说明，并简要说明了它们的用法：
| 损失函数 | 解释 |
|---|---|
| 交叉熵 | 这用于分类模型，其中模型的输出为 0-1 之间的概率。 这是对数损失函数，也称为对数损失。 当理想的模型的概率接近 1.0 时， 交叉熵损失降低。 |
| MAE（L1） | 计算误差绝对值的均值。 由于它仅使用绝对值，因此不会将较大误差的权重放大。 当较大误差与较小误差相比可以容忍时，这很有用。 |
| MSE（L2） | 计算误差的平方根的均值。 放大较大误差的权重。 当不希望出现较大误差时，这很明智。 |
| Hinge | 这是用于线性分类器模型（例如支持向量机）的损失函数。 |
| Huber | 这是回归模型的损失函数。 它与 MSE 非常相似，但对异常值不敏感。 |
| Kullback-Leibler | Kullback-Leibler 散度衡量两个概率分布之间的差异 。 KL 损失函数在 T 分布随机邻居嵌入算法中大量使用。 |
在机器学习算法中，损失函数在更新变量权重时至关重要。 假设您使用反向传播训练神经网络。 在每次迭代中，都会计算出总误差。 然后，权重被更新以便最小化总误差。 因此，使用正确的损失函数会直接影响您的机器学习算法的性能，因为它直接影响模型参数。 在下一章中，我们将从住房数据中的单个变量开始简单的线性回归。
使用梯度下降的单变量线性回归
在本小节中，我们将为波士顿住房数据集实现单变量线性回归，该数据已在上一章中用于探索性数据分析。 在拟合回归线之前，让我们导入必要的库并按以下方式加载数据集：
In [1]: import numpy as np
import pandas as pd
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
%matplotlib inline
In [2]: from sklearn.datasets import load_boston
dataset = load_boston()
samples , label, feature_names = dataset.data, dataset.target, dataset.feature_names
In [3]: bostondf = pd.DataFrame(dataset.data)
bostondf.columns = dataset.feature_names
bostondf['Target Price'] = dataset.target
bostondf.head()
Out[3]: CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Target Price
0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0
1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6
2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7
3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4
4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2
与上一章相比，我们使用 Pandas 数据帧代替了 numpy 数组，以便向您展示数据帧的用法，因为它也是一个非常方便的数据结构。 从技术上讲，在大多数情况下，如果仅将数值存储在 numpy 数组或 Pandas 数据帧中，则没有任何区别。 让我们将目标值添加到数据框中，并使用散点图查看RM函数与目标值之间的关系：
In [4]: import matplotlib.pyplot as plt
bostondf.plot(x='RM', y='Target Price', style= 'o')
plt.title('RM vs Target Price')
plt.ylabel('Target Price')
plt.show()
从该图可以看出，每个住宅的平均房间数（RM）与房价之间存在正相关，正如预期的那样。 现在，我们将看到这种关系的强度，并尝试使用这种关系来预测房价。
想象一下，您对线性回归的了解非常有限。 假设您只是熟悉方程式，但不知道误差函数是什么，为什么我们需要迭代，什么是梯度下降，以及为什么要在线性回归模型中使用它。 在这种情况下，您要做的就是简单地开始传递系数的一些初始值，并在计算预测值之前截取方程。
在计算了几个预测值之后，您可以将它们与实际值进行比较，看看您与实际情况有多远。 下一步将是更改系数或截距，或者同时进行这两项操作，以查看是否可以使结果更接近实际值。 如果您对此过程感到满意，这就是我们的算法将以更智能的方式进行的操作。
为了更好地理解线性回归模型步骤，我们将代码分成几个块。 首先，让我们创建一个函数，该函数作为回归线的结果返回预测值：
In [5]: def prediction(X, coefficient, intercept):
return X*coefficient + intercept
前面的函数计算线性回归模型如下：
然后，我们需要一个成本函数，该函数将在每次迭代中进行计算。 作为cost_function，我们将使用均方误差，它是预测值与实际值之间总平方差的平均值：
In [6]: def cost_function(X, Y, coefficient, intercept):
MSE = 0.0
for i in range(len(X)):
MSE += (Y[i] - (coefficient*X[i] + intercept))**2
return MSE / len(X)
我们的最后一个代码块将用于更新权重。 当我们讨论权重时，不仅涉及自变量的系数，还涉及截距。 拦截也称为偏差。 为了逻辑上更新权重，我们需要一种迭代优化算法，该算法可以找到给定函数的最小值。 在此示例中，我们将使用梯度下降方法来最小化每次迭代中的损失函数。 让我们一步一步地发现梯度下降的作用。
首先，我们应该初始化权重（截距和系数）并计算均方误差。 然后，我们需要计算梯度，这意味着查看当改变权重时均方误差如何变化。
为了以更智能的方式更改权重，我们需要了解我们必须更改系数和截距的方向。 这意味着我们在更改权重时应计算误差函数的梯度。 我们可以通过对系数和截距采用损失函数的偏导数来计算梯度。
在单变量线性回归中，我们只有一个系数。 在计算偏导数之后，该算法将调整权重并重新计算均方误差。 此过程将重复进行，直到更新的权重不再减小均方误差为止：
In [7]: def update_weights(X, Y, coefficient, intercept, learning_rate):
coefficient_derivative = 0
intercept_derivative = 0
for i in range(len(X)):
coefficient_derivative += -2*X[i] * (Y[i] - (coefficient*X[i] + intercept))
intercept_derivative += -2*(Y[i] - (coefficient*X[i] + intercept))
coefficient -= (coefficient_derivative / len(X)) * learning_rate
intercept -= (intercept_derivative / len(X)) * learning_rate
return coefficient, intercept
前面的代码块定义了更新权重的函数，然后返回更新的系数并进行拦截。 此函数中的另一个重要参数是learning_rate。 学习速度将决定变化的幅度：
在上图中，您可以看到我们的损失函数为x = y^2的形状，因为它是平方误差的总和。 通过梯度下降，我们试图找到最小损耗，如图所示，其中偏导数非常接近零。 如前所述，我们通过初始化权重来开始算法，这意味着我们从远离最小值的点开始。 在每次迭代中，我们将更新权重，从而减少损失。 这意味着给定足够的迭代次数，我们将收敛到全局最小值。 学习速度将决定这种融合发生的速度。
换句话说，当我们更新权重时，高学习率看起来就像是从一个点到另一个点的巨大跳跃。 低学习率将慢慢接近全局最小值（所需的最小损失点）。 由于学习率是一个超参数，因此我们需要在运行它之前进行设置：
那么，我们需要设定学习率的大还是小？ 答案是我们应该找到最佳速率。 如果我们设定较高的学习率，我们的算法将超出最小值。 我们很容易错过最小值，因为这些跳跃永远不会让我们的算法收敛到全局最小值。 另一方面，如果我们将学习率设置得太小，则可能需要进行大量迭代才能收敛。
具有先前的代码块，是时候编写主函数了。 主函数应遵循我们前面讨论的流程：
然后，主函数的代码块应如下所示：
In [8]: def train(X, Y, coefficient, intercept, LearningRate, iteration):
cost_hist = []
for i in range(iteration):
coefficient, intercept = update_weights(X, Y, coefficient, intercept, learning_cost = cost_function(X, Y, coefficient, intercept)
cost_hist.append(cost)
return coefficient, intercept, cost_hist
现在，我们已经定义了所有代码块，并且可以运行单变量模型了。 在运行train()函数之前，我们需要设置系数和截距的超参数和初始值。 您可以创建变量，设置值，然后将这些变量作为函数的参数，如下所示：
In [9]: learning_rate = 0.01
iteration = 10001
coefficient = 0.3
intercept = 2
X = bostondf.iloc[:, 5:6].values
Y = bostondf.iloc[:, 13:14].values
coefficient, intercept, cost_history = train(X, Y, coefficient, intercept, learning_rate, iteration)
或者，可以在调用函数时将这些值作为关键字参数传递：
coefficient, intercept, cost_history = train(X, Y, coefficient, intercept = 2, learning_rate = 0.01, iteration = 10001)
我们的主函数将返回两个数组，这将为我们提供拦截和系数的最终值。 此外，主函数将返回损失值列表，这是每次迭代的结果。 这些是每次迭代中均方误差的结果。 拥有此列表以跟踪每次迭代中损耗的变化非常有用：
In [10]: coefficient
Out[10]: array([8.57526661])
In [11]: intercept
Out[11]: array([-31.31931428])
In [12]: cost_history
array([54.18545801]),
array([54.18036786]),
array([54.17528017]),
array([54.17019493]),
array([54.16511212]),
array([54.16003177]),
array([54.15495385]),
array([54.14987838]),
array([54.14480535]),
array([54.13973476]),
array([54.13466661]),
array([54.12960089]),
array([54.12453761]),
array([54.11947677]),
array([54.11441836]),
array([54.10936238]),
array([54.10430883]),
...]
从系数值中可以看出，RM的每增加 1 单位，住房价格就会升至 8,575 美元左右。 现在，让我们通过将计算出的截距和系数注入回归公式中来计算预测值。 然后，我们可以绘制线性回归线并查看其如何拟合我们的数据：
In [13]: y_hat = X*coefficient + intercept
plt.plot(X, Y, 'bo')
plt.plot(X, y_hat)
plt.show()
在本节中，我们通过选择一个变量来应用单变量模型。 在以下小节中，我们将通过在模型中添加更多自变量来执行多元线性回归模型，这意味着我们将有多个系数可进行优化以实现最佳拟合。
使用线性回归建模房价
在本节中，我们将对同一数据集执行多元线性回归。 与上一节相反，我们将使用 sklearn 库向您展示执行线性回归模型的几种方法。 在开始线性回归模型之前，我们将使用trimboth()方法从两侧按比例修剪数据集。 通过这样做，我们将消除异常值：
In [14]: import numpy as np
import pandas as pd
from scipy import stats
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import LinearRegression
In [15]: from sklearn.datasets import load_boston
dataset = load_boston()
In [16]: samples , label, feature_names = dataset.data, dataset.target, dataset.feature_names
In [17]: samples_trim = stats.trimboth(samples, 0.1)
label_trim = stats.trimboth(label, 0.1)
In [18]: print(samples.shape)
print(label.shape)
(506, 13)
(506,)
In [19]: print(samples_trim.shape)
print(label_trim.shape)
(406, 13)
(406,)
正如您在前面的代码块中看到的那样，由于我们将左右数据修剪了 10% ，每个属性和标签列的样本量从 506 减少到 406：
In [20]: from sklearn.model_selection import train_test_split
samples_train, samples_test, label_train, label_test = train_test_split(samples_trim, In [21]: In [21]: print(samples_train.shape)
print(samples_test.shape)
print(label_train.shape)
print(label_test.shape)
(324, 13)
(82, 13)
(324,)
(82,)
In [22]: regressor = LinearRegression()
regressor.fit(samples_train, label_train)
Out[22]: LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
In [23]: regressor.coef_
Out[23]: array([ 2.12924665e-01, 9.16706914e-02, 1.04316071e-01, -3.18634008e-14,
5.34177385e+00, -7.81823481e-02, 1.91366342e-02, 2.81852916e-01,
3.19533878e-04, -4.24007416e-03, 1.94206366e-01, 3.96802252e-02,
3.81858253e-01])
In [24]: regressor.intercept_
Out[24]: -6.899291747292615
然后，我们使用train_test_split()方法将数据集拆分为训练和测试。 这种方法在机器学习算法中非常普遍。 您将数据分为两组，然后让模型训练（学习），然后再使用数据的另一部分测试模型。 使用这种方法的原因是为了减少验证模型时的偏差。 每个系数代表当我们将这些样本增加一个单位时，预期目标值将改变多少：
In [25]: label_pred = regressor.predict(samples_test)
In [26]: plt.scatter(label_test, label_pred)
plt.xlabel("Prices")
plt.ylabel("Predicted Prices")
plt.title("Prices vs Predicted Prices")
plt.axis("equal")
Out[26]: (11.770143369175626, 34.22985663082437, 10.865962968036989, 34.20549738482051)
让我们在前面的代码块中测试我们的模型。 为了在数据集上运行我们的模型，我们可以使用predict()方法。 此方法将在模型上运行给定的数据集并返回结果。 理想情况下，我们期望该点以x = y线形分布，换言之，例如具有 45 度角的直线。 这是一个完美的一对一匹配预测。 我们不能说我们的预测是完美的，但是仅通过查看散点图，我们可以得出结论，该预测不是一个坏的预测。 我们可以查看有关模型性能的两个重要指标，如下所示：
In [27]: from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
mse = mean_squared_error(label_test, label_pred)
r2 = r2_score(label_test, label_pred)
print(mse)
print(r2)
2.032691267250478
0.9154474686142619
第一个是均方误差。 您可能已经注意到，当我们在模型中添加更多自变量时，它会大大降低。 第二个是 R 方，它是确定系数，或者换句话说，是回归得分。 确定系数可以解释您的自变量解释了因变量有多少方差。 在我们的模型中，结果为 0.91，这说明我们的 13 个特征可以解释 91% 的房价差异。 如您所见，sklearn 具有许多有用的线性回归内置函数，这些函数可以加快模型构建过程。 另一方面，您也可以仅使用 numpy 轻松构建线性回归模型，这可以为您提供更大的灵活性，因为您可以控制算法的每个部分。
总结
线性回归是对连续变量之间的关系进行建模的最常用技术之一。 这种方法的应用在工业上非常广泛。 我们开始对本书的一部分进行线性回归建模，这不仅是因为它非常流行，还因为它是一种相对简单的技术，并且包含了几乎每种机器学习算法都具有的大多数要素。
在本章中，我们了解了监督学习和无监督学习，并使用波士顿房屋数据集建立了线性回归模型。 我们谈到了不同的重要概念，例如超参数，损失函数和梯度下降。 本章的主要目的是为您提供足够的知识，以便您可以建立和调整线性回归模型并逐步了解它的作用。 我们研究了两个使用单变量和多元线性回归的实际案例。 您还体验了 numpy 和 sklearn 在线性回归中的用法。 我们强烈建议您对不同的数据集进行进一步练习，并检查更改超参数时结果如何变化。
在下一章中，我们将学习聚类方法，并通过批发分销商数据集上的示例进行实践。
五、使用 NumPy 对批发分销商的客户进行聚类
通过查看 NumPy 在各种使用案例中的使用，您肯定可以提高自己的技能。 本章介绍的分析类型与迄今为止所见不同。 聚类是一种无监督的学习技术，用于理解和捕获数据集中的各种形式。 由于您没有标签来监督您的学习算法，因此在很多情况下，可视化是关键，这就是为什么您也会看到各种可视化技术的原因。
在本章中，我们将介绍以下主题：
- 无监督学习和聚类
- 超参数
- 扩展简单算法来聚类批发分销商的客户
无监督学习和聚类
让我们用一个例子快速回顾一下监督学习。 在训练机器学习算法时，您可以通过提供标签来观察和指导学习。 考虑下面的数据集，其中每一行代表一个客户，每一列代表一个不同的特征，例如年龄，性别，收入，职业 ，任期和城市。 看看这个表：
您可能需要执行不同类型的分析。 其中一项可能是预测哪些客户可能会离开，即客户流失分析。 为此，您需要根据每个客户的历史记录为其标记标签，以指示哪些客户已离开或留下，如下表所示：
您的算法将根据客户的标签了解客户的特征。 算法将学习离开或留下的客户的特征，并且，当您希望根据这些特征为新客户评分时，算法将根据该学习进行预测。 这称为监督学习。
关于此数据集可能要问的另一个问题可能是：此数据集中存在多少个客户组，以便其组中的每个客户都与其他客户相似，而与属于其他组的客户不同。 流行的聚类算法（例如 K 均值）可以帮助您解决这一问题。 例如，一旦 K 均值将客户分配到不同的集群，则一个集群可以大体上包括 30 岁以下且以 IT 作为职业的客户，而另一个集群可以大体上包括 60 岁以上的客户并且职业为老师。 您无需标记数据集即可执行此分析，因为算法足以查看记录并确定它们之间的相似性。 由于没有监督，因此这种学习称为无监督学习。
进行此类分析时，首先可视化数据集会很有帮助。 您可以从可用的数据集开始构建您的处理和建模工作流程。 以下代码段显示了如何使用plotly可视化三维数据集。plotly是一个库，可让您绘制许多不同的交互式图表以进行探索性分析，并使数据探索更加容易。
首先，您需要使用以下代码段安装必要的库：
## Installing necessary libraries with pip
!pip install plotly --user
!pip install cufflinks -user
然后，使用以下代码导入必要的库：
## Necessary imports
import os
import sys
import numpy as np
import pandas
import matplotlib.pyplot as plt
%matplotlib inline
import plotly.plotly as py
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import cufflinks as cf
import plotly.graph_objs as go
init_notebook_mode(connected=True)
sys.path.append("".join([os.environ["HOME"]]))
您将使用sklearn.datasets模块中可用的iris数据集，如下所示：
from sklearn.datasets import load_iris
iris_data = load_iris()
iris数据具有四个特征； 它们如下：
iris_data.feature_names
['sepal length (cm)',
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)']
首先，让我们检查一下前两个特征：
x = [v[0] for v in iris_data.data]
y = [v[1] for v in iris_data.data]
创建一个trace，然后创建数据和图形，如下所示：
trace = go.Scatter(
x = x,
y = y,
mode = 'markers'
)
layout= go.Layout(
title= 'Iris Dataset',
hovermode= 'closest',
xaxis= dict(
title= 'sepal length (cm)',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'sepal width (cm)',
ticklen= 5,
gridwidth= 2,
),
showlegend= False
)
data = [trace]
fig= go.Figure(data=data, layout=layout)
plot(fig)
如下图所示，这将为您提供以下输出：
您可以继续查看其他变量，但是要更好地理解一张图表中要素之间的关系，可以使用scatterplot矩阵。 创建pandas.DataFrame与plotly结合使用会更加方便：
import pandas as pd
df = pd.DataFrame(iris_data.data,
columns=['sepal length (cm)',
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)'])
df['class'] = [iris_data.target_names[i] for i in iris_data.target]
df.head()
使用plotly图形工厂，可以绘制scatterplot矩阵，如下所示：
import plotly.figure_factory as ff
fig = ff.create_scatterplotmatrix(df, index='class', diag='histogram', size=10, height=800, width=800)
plot(fig)
这将为您提供以下图表，如下图所示：
乍一看，花瓣长度，花瓣宽度和萼片长度似乎是建模的不错选择。 您可以通过以下代码使用 3D 图表进一步检查该数据集：
## Creating data for the plotly
trace1 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
trace2 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
mode='markers',
marker=dict(
color='rgb(#3742fa)',
size=12,
symbol='circle',
line=dict(
color='rgb(204, 204, 204)',
width=1
),
opacity=0.9
)
)
trace3 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
mode='markers',
marker=dict(
color='rgb(#ff4757)',
size=12,
symbol='circle',
line=dict(
color='rgb(104, 74, 114)',
width=1
),
opacity=0.9
)
)
data = [trace1, trace2, trace3]
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'sepal length (cm)'),
yaxis = dict(
title= 'petal length (cm)'),
zaxis = dict(
title= 'petal width (cm)'),),
)
fig = go.Figure(data=data, layout=layout)
plot(fig)
这将为您提供以下交互式图表，如下图所示：
Interactive plotting of petal length, petal width, and sepal length
使用这些图表，您可以更好地理解数据并为建模做准备。
超参数
超参数可以被视为确定模型各种属性之一的高级参数，例如复杂性，训练行为和学习率。 这些参数自然与模型参数有所不同，因为它们需要在训练开始之前设置。
例如，K 均值或 K 最近邻中的 k 是这些算法的超参数。 K 均值中的 K 表示要找到的聚类数，K 最近邻中的 k 表示用于进行预测的最近记录数。
调整超参数是任何机器学习项目中提高预测性能的关键步骤。 有多种调优技术，例如网格搜索，随机搜索和贝叶斯优化，但这些技术不在本章范围之内。
让我们通过以下屏幕快照快速浏览 scikit-learn 库中的 K 均值算法参数：
如您所见，有许多参数可以使用，并且至少应查看算法的函数签名，以查看运行算法之前的选项。
让我们一起玩吧。 基线模型将与样本数据一起使用，几乎具有默认设置，如下所示：
from sklearn.datasets.samples_generator import make_blobs
X, y = make_blobs(n_samples=20, centers=3, n_features=3, random_state=42)
k_means = KMeans(n_clusters=3)
y_hat = k_means.fit_predict(X)
y_hat保留集群的成员身份信息，这与原始标签相同，如您在此处看到的：
y_hat
## array([0, 2, 1, 1, 1, 0, 2, 0, 0, 0, 2, 0, 1, 2, 1, 2, 2, 1, 0, 1],
dtype=int32)
y
## array([0, 2, 1, 1, 1, 0, 2, 0, 0, 0, 2, 0, 1, 2, 1, 2, 2, 1, 0, 1])
您可以使用不同的选项来查看它如何影响训练和预测。
损失函数
损失函数通过测量误差来帮助算法在训练过程中更新模型参数，这是预测性能的指标。 损失函数通常表示为：
其中L测量预测值与实际值之间的差。 在训练过程中，此误差被最小化。 不同的算法具有不同的损失函数，迭代次数将取决于收敛条件。
例如，K 均值的损失函数使点与最接近的簇均值之间的平方距离最小化，如下所示：
您将在以下部分中看到详细的实现。
为单个变量实现我们的算法
让我们为单个变量实现 K 均值算法。 您将从一维向量开始，该向量具有 20 条记录，如下所示：
data = [1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27]
trace1 = go.Scatter(
x=data,
y=[0 for x in data],
mode='markers',
name='Data',
marker=dict(
size=12
)
)
layout = go.Layout(
title='1D vector',
)
traces = [trace1]
fig = go.Figure(data=traces, layout=layout)
plot(fig)
如下图所示，将输出以下图表：
我们的目标是找到在数据中可见的3簇。 为了开始实施 K 均值算法，您需要通过选择随机索引来初始化集群中心，如下所示：
n_clusters = 3
c_centers = np.random.choice(X, n_clusters)
print(c_centers)
## [ 1 22 26]
接下来，您需要计算每个点与聚类中心之间的距离，因此请使用以下代码：
deltas = np.array([np.abs(point - c_centers) for point in X])
deltas
array([[ 7, 26, 10],
[ 6, 25, 9],
[ 5, 24, 8],
[ 6, 25, 9],
[ 7, 26, 10],
[ 5, 24, 8],
[ 1, 18, 2],
[ 0, 19, 3],
[ 3, 16, 0],
[ 4, 15, 1],
[ 2, 17, 1],
[ 3, 16, 0],
[ 6, 13, 3],
[17, 2, 14],
[18, 1, 15],
[16, 3, 13],
[22, 3, 19],
[14, 5, 11],
[16, 3, 13],
[19, 0, 16]])
现在，可以使用以下代码来计算集群成员资格：
deltas.argmin(1)
## array([0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1])
现在，您需要使用以下代码来计算记录和群集中心之间的距离：
c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(3)])
print(c_centers)
## [ 3.625 25.42857143 11.6 ]
这是一次迭代； 您可以继续计算新的群集中心，直到没有任何改善为止。
您可以编写一个函数来包装所有这些功能，如下所示：
def Kmeans_1D(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([np.abs(point - c_centers) for point in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([np.abs(point - new_c_centers) for point in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.all(c_centers == new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
c_centers, labels = Kmeans_1D(X, 3)
print(c_centers, labels)
## [11.16666667 25.42857143 2.85714286] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
让我们使用以下代码绘制集群中心的图表：
trace1 = go.Scatter(
x=X,
y=[0 for num in X],
mode='markers',
name='Data',
marker=dict(
size=12
)
)
trace2 = go.Scatter(
x = c_centers,
y = [0 for num in X],
mode='markers',
name = 'Cluster centers',
marker = dict(
size=12,
color = ('rgb(122, 296, 167)')
)
)
layout = go.Layout(
title='1D vector',
)
traces = [trace1, trace2]
fig = go.Figure(data=traces, layout=layout)
plot(fig)
看下图。 给定的代码将输出以下内容：
您可以清楚地看到，可以为每个元素分配一个聚类中心。
修改我们的算法
现在，您已经了解了单个变量的 K 均值内部，可以将该实现扩展到多个变量，并将其应用于更实际的数据集。
本节中使用的数据集来自 UCI 机器学习存储库，它包括批发分销商的客户信息。 有 440 个具有八项特征的客户。 在下面的列表中，前六个特征与相应产品的年度支出相关，第七个特征显示了购买该产品的渠道，第八个特征显示了该地区：
FRESHMILKGROCERYFROZENDETERGENTS_PAPERDELICATESSENCHANNELREGION
首先下载数据集并将其读取为numpy数组：
from numpy import genfromtxt
wholesales_data = genfromtxt('Wholesale customers data.csv', delimiter=',', skip_header=1)
您可以快速查看数据。 这里是：
print(wholesales_data[:5])
[[2.0000e+00 3.0000e+00 1.2669e+04 9.6560e+03 7.5610e+03 2.1400e+02
2.6740e+03 1.3380e+03]
[2.0000e+00 3.0000e+00 7.0570e+03 9.8100e+03 9.5680e+03 1.7620e+03
3.2930e+03 1.7760e+03]
[2.0000e+00 3.0000e+00 6.3530e+03 8.8080e+03 7.6840e+03 2.4050e+03
3.5160e+03 7.8440e+03]
[1.0000e+00 3.0000e+00 1.3265e+04 1.1960e+03 4.2210e+03 6.4040e+03
5.0700e+02 1.7880e+03]
[2.0000e+00 3.0000e+00 2.2615e+04 5.4100e+03 7.1980e+03 3.9150e+03
1.7770e+03 5.1850e+03]]
选中shape将显示行数和变量数，如下所示：
wholesales_data.shape
## (440, 8)
该数据集具有440个记录，每个记录都有8个特征。
通过使用以下代码来规范化数据集是一个好主意：
wholesales_data_norm = wholesales_data / np.linalg.norm(wholesales_data)
print(wholesales_data_norm[:5])
[[ 1\. 0\. 0.30168043 1.06571214 0.32995207 -0.46657183
0.50678671 0.2638102 ]
[ 1\. 0\. -0.1048095 1.09293385 0.56599336 0.08392603
0.67567015 0.5740085 ]
[ 1\. 0\. -0.15580183 0.91581599 0.34441798 0.3125889
0.73651183 4.87145892]
[ 0\. 0\. 0.34485007 -0.42971408 -0.06286202 1.73470839
-0.08444172 0.58250708]
[ 1\. 0\. 1.02209184 0.3151708 0.28726 0.84957326
0.26205579 2.98831445]]
您可以使用以下代码将数据集读取到pandas.DataFrame中：
import pandas as pd
df = pd.DataFrame(wholesales_data_norm,
columns=['Channel',
'Region',
'Fresh',
'Milk',
'Grocery',
'Frozen',
'Detergents_Paper',
'Delicatessen'])
df.head(10)
让我们创建一个scatterplot矩阵以更仔细地查看数据集。 看一下代码：
fig = ff.create_scatterplotmatrix(df, diag='histogram', size=7, height=1200, width=1200)
plot(fig)
您还可以通过运行以下命令来检查要素之间的相关性：
df.corr()
这将为您提供一个关联表，如下所示：
您还可以使用seaborn热图，如下所示：
import seaborn as sns; sns.set()
ax = sns.heatmap(df.corr(), annot=True)
Correlations between the features
您可以看到特征之间存在一些很强的相关性，例如Grocery和Detergents_Paper之间的相关性。
让我们使用以下代码绘制Grocery，Detergents_Paper和Milk三个特征：
## Creating data for the plotly
trace1 = go.Scatter3d(
# Extracting data based on label
x=df['Grocery'],
y=df['Detergents_Paper'],
z=df['Milk'],
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'Grocery'),
yaxis = dict(
title= 'Detergents_Paper'),
zaxis = dict(
title= 'Milk'),),
)
data = [trace1]
fig = dict(data=data, layout=layout)
plot(fig)
现在，您将继续扩展已为更高维度实现的 K 均值算法。 首先，您可以从数据集中删除Channel和Region，如下所示：
df = df[[col for col in df.columns if col not in ['Channel', 'Region']]]
df.head(10)
在实现方面，您还可以使用np.linalg.norm来计算距离，这实际上取决于您使用哪种距离函数。 另一个替代方法是scipy.spatial中的distance.euclidean，如下所示：
def Kmeans_nD(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([[np.linalg.norm(i - c) for c in c_centers] for i in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean(axis=0) for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([[np.linalg.norm(i - c) for c in new_c_centers] for i in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.array_equal(c_centers, new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
Grocery和Detergents_Paper将用于聚类，并且k将设置为3。 通常，应使用外观检查或弯头方法确定k，如下所示：
centers, labels = Kmeans_nD(df[['Grocery', 'Detergents_Paper']].values, 3)
现在，您可以使用以下命令在数据集中再添加一列：
df['labels'] = labels
您可以使用以下代码来首先可视化结果，以查看结果是否有意义：
## Creating data for the plotly
trace1 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 0]['Grocery'],
y=df[df['labels'] == 0]['Detergents_Paper'],
mode='markers',
name='clust_1',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
trace2 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 1]['Grocery'],
y=df[df['labels'] == 1]['Detergents_Paper'],
mode='markers',
name='clust_2',
marker=dict(
color='rgb(#3742fa)',
size=12,
symbol='circle',
line=dict(
color='rgb(204, 204, 204)',
width=1
),
opacity=0.9
)
)
trace3 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 2]['Grocery'],
y=df[df['labels'] == 2]['Detergents_Paper'],
mode='markers',
name='clust_3',
marker=dict(
color='rgb(#ff4757)',
size=12,
symbol='circle',
line=dict(
color='rgb(104, 74, 114)',
width=1
),
opacity=0.9
)
)
data = [trace1, trace2, trace3]
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'Grocery'),
yaxis = dict(title= 'Detergents_Paper'),
)
)
fig = go.Figure(data=data, layout=layout)
plot(fig)
Plotting the clusters
乍一看，集群看起来很合理，并将最终取决于领域知识支持的解释。
您可以使用以下代码轻松查看每个集群每个特征的平均支出：
df.groupby('labels').mean()
从这个简单的集群中可以看出，第三集群在Milk，Grocery和Detergents_Paper上的支出最高。 第二个集群的支出较低，第一个集群倾向于Milk，Grocery和Detergents_Paper，因此k=2也可以使用。
总结
在本章中，您学习了无监督学习的基础知识，并使用 K 均值算法进行聚类。
有许多聚类算法显示不同的行为。 当涉及到无监督学习算法时，可视化是关键，并且您已经看到了几种不同的方式来可视化和检查数据集。
在下一章中，您将学习 NumPy 常用的其他库，例如 SciPy，Pandas 和 scikit-learn。 这些都是从业者工具包中的重要库，它们相互补充。 您会发现自己将这些库与 NumPy 一起使用，因为每个库都会使某些任务变得更容易。 因此，重要的是要了解有关 Python 数据科学堆栈的更多信息。
六、NumPy，SciPy，Pandas 和 Scikit-Learn
到目前为止，您应该能够使用 NumPy 编写小型实现。 在整个章节中，我们旨在提供使用其他库的示例，在本章中，我们应退后一步，看看可以与 NumPy 一起用于项目的周围库。
本章将介绍其他 Python 库如何对 NumPy 进行补充。 我们将研究以下主题：
- NumPy 和 SciPy
- NumPy 和 Pandas
- SciPy 和 Scikit-learn
NumPy 和 SciPy
到目前为止，您已经看到了许多有关 NumPy 用法的示例，而只有少数 SciPy。 NumPy 具有数组数据类型，它允许您执行各种数组操作，例如排序和整形。
NumPy 具有一些数值算法，可用于执行诸如计算范数，特征值和特征向量之类的任务。 但是，如果数值算法是您的重点，则理想情况下应使用 SciPy，因为它包含更全面的算法集以及最新版本的算法。 SciPy 有许多有用的子程序包可用于某些类型的分析。
以下列表将使您对子软件包有一个整体的了解：
Cluster：此子程序包包含聚类算法。 它具有两个子模块vq和hierarchy。vq模块提供用于 K 均值聚类的函数。 层次结构模块包括用于层次结构聚类的函数。Fftpack：此子程序包包含用于快速傅立叶变换的函数和算法，以及差分和伪差分算符。Interpolate：此子程序包提供用于单变量和多变量插值的函数：1D 和 2D 样条曲线。Linalg：此子程序包提供用于线性代数的函数和算法，例如matrix运算和函数，特征值和-向量计算，矩阵分解，矩阵方程求解器和特殊矩阵。Ndimage：此子程序包提供用于多维图像处理的函数和算法，例如滤镜，插值，测量和形态。Optimize：此子程序包提供函数和算法，用于函数局部和全局优化，函数拟合，求根和线性编程。Signal：此子程序包提供信号处理的函数和算法，例如卷积，B 样条，滤波，连续和离散时间线性系统，波形，小波和频谱分析。Stats：此子程序包提供概率分布，例如连续分布，多元分布和离散分布，以及可以找到均值，众数，方差，偏度，峰度和相关系数的统计函数。
让我们来看一下其中的一个子包。 以下代码显示了用于群集分析的cluster程序包：
Scipy.cluster
%matplotlib inline
import matplotlib.pyplot as plt
## Import ndimage to read the image
from scipy import ndimage
## Import cluster for clustering algorithms
from scipy import cluster
## Read the image
image = ndimage.imread("cluster_test_image.jpg")
## Image is 1000x1000 pixels and it has 3 channels.
image.shape
(1000, 1000, 3)
这将为您提供以下输出：
array([[[30, 30, 30],
[16, 16, 16],
[14, 14, 14],
...,
[14, 14, 14],
[16, 16, 16],
[29, 29, 29]],
[[13, 13, 13],
[ 0, 0, 0],
[ 0, 0, 0],
...,
[ 0, 0, 0],
[ 0, 0, 0],
[12, 12, 12]],
[[16, 16, 16],
[ 3, 3, 3],
[ 1, 1, 1],
...,
[ 0, 0, 0],
[ 2, 2, 2],
[16, 16, 16]],
...,
[[17, 17, 17],
[ 3, 3, 3],
[ 1, 1, 1],
...,
[34, 26, 39],
[27, 21, 33],
[59, 55, 69]],
[[15, 15, 15],
[ 2, 2, 2],
[ 0, 0, 0],
...,
[37, 31, 43],
[34, 28, 42],
[60, 56, 71]],
[[33, 33, 33],
[20, 20, 20],
[17, 17, 17],
...,
[55, 49, 63],
[47, 43, 57],
[65, 61, 76]]], dtype=uint8)
在这里，您可以看到该图：
plt.figure(figsize = (15,8))
plt.imshow(image)
您可以从前面的代码块中获得以下图表：
使用以下代码将图像数组转换为二维数据集：
x, y, z = image.shape
image_2d = image.reshape(x*y, z).astype(float)
image_2d.shape
(1000000, 3)
image_2d
array([[30., 30., 30.],
[16., 16., 16.],
[14., 14., 14.],
...,
[55., 49., 63.],
[47., 43., 57.],
[65., 61., 76.]])
## kmeans will return cluster centers and the distortion
cluster_centers, distortion = cluster.vq.kmeans(image_2d, k_or_guess=2)
print(cluster_centers, distortion)
[[179.28653454 179.30176248 179.44142117]
[ 3.75308484 3.83491111 4.49236356]] 26.87835069294931
image_2d_labeled = image_2d.copy()
labels = []
from scipy.spatial.distance import euclidean
import numpy as np
for i in range(image_2d.shape[0]):
distances = [euclidean(image_2d[i], center) for center in cluster_centers]
labels.append(np.argmin(distances))
plt.figure(figsize = (15,8))
plt.imshow(cluster_centers[labels].reshape(x, y, z))
您从前面的代码中获得以下输出：
SciPy 和 NumPy 和线性回归
您已经了解了如何使用 NumPy 从头开始编写线性回归算法。Scipy.stats模块具有linregress函数，用于计算斜率，截距，相关系数（r 值），两侧 p 值以及标准差估计，如下所示：
from sklearn import datasets
%matplotlib inline
import matplotlib.pyplot as plt
## Boston House Prices dataset
boston = datasets.load_boston()
x = boston.data
y = boston.target
boston.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')
x.shape
(506, 13)
y.shape
(506,)
## We will consider "lower status of population" as independent variable for its importance
lstat = x[0:,-1]
lstat.shape
(506,)
from scipy import stats
slope, intercept, r_value, p_value, std_err = stats.linregress(lstat, y)
print(slope, intercept, r_value, p_value, std_err)
-0.9500493537579909 34.55384087938311 -0.737662726174015 5.081103394387796e-88 0.03873341621263942
print("r-squared:", r_value**2)
r-squared: 0.5441462975864798
plt.plot(lstat, y, 'o', label='original data')
plt.plot(lstat, intercept + slope*lstat, 'r', label='fitted line')
plt.legend()
plt.show()
我们从前面代码的输出中获得以下图表，如下图所示：
您还可以查看平均房间数与房价之间的关系。 以下代码块打印出性能指标：
rm = x[0:,5]
slope, intercept, r_value, p_value, std_err = stats.linregress(rm, y)
print(slope, intercept, r_value, p_value, std_err)
print("r-squared:", r_value**2)
## 9.102108981180308 -34.670620776438554 0.6953599470715394 2.48722887100781e-74 0.4190265601213402
## r-squared: 0.483525455991334
以下代码块绘制了拟合线：
plt.plot(rm, y, 'o', label='original data')
plt.plot(rm, intercept + slope*rm, 'r', label='fitted line')
plt.legend()
plt.show()
我们从前面的代码中获得以下输出，如下图所示：
NumPy 和 Pandas
考虑一下时，NumPy 是一个相当低级的数组操作库，大多数其他 Python 库都写在它的顶部。
这些库之一是pandas，它是一个高级数据处理库。 浏览数据集时，通常会执行诸如计算描述性统计数据，按特定特征分组以及合并之类的操作。pandas库具有许多友好的函数来执行这些各种有用的操作。
在此示例中，我们使用糖尿病数据集。sklearn.datasets中的糖尿病数据集使用零均值和 L2 范数标准化。
该数据集包含 442 条记录，这些记录具有 10 个特征：年龄，性别，体重指数，平均血压和 6 个血清测量值。
目标代表采取这些基准措施后的疾病进展。 您可以在 web 和相关论文 中查看数据描述。
我们从操作开始，如下所示：
import pandas as pd
from sklearn import datasets
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
diabetes = datasets.load_diabetes()
df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)
diabetes.feature_names
['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']
df.head(10)
我们从前面的代码中获得以下输出，如下表所示：
这段代码向您展示了如何在数据帧中创建目标列：
df['Target'] = diabetes.target
df.head(10)
Pandas 帮助我们轻松地处理表格数据，并通过各种辅助方法和可视化支持我们的分析。 看一下代码：
## Descriptive statistics
df.describe()
我们从前面的代码中获得以下输出，如下表所示：
通过使用以下代码行，看看目标的分布方式：
plt.hist(df['Target'])
下图显示了上一行的输出：
您可以看到目标变量向右倾斜。 看一下这段代码：
## Since 'sex' is categorical, excluding it from numerical columns
numeric_cols = [col for col in df.columns if col != 'sex']
numeric_cols
## ['age', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6', 'Target']
## You can have a look at variable distributions individually, but there's a better way
df[numeric_cols].hist(figsize=(20, 20), bins=30, xlabelsize=12, ylabelsize=12)
## You can also choose create dataframes for numerical and categorical variables
前一个代码块的输出：
Feature distributions
您可以检查其中一些特征的分布，并确定其中哪些看起来相似。 对于此示例，特征s1，s2和s6似乎具有相似的分布，如从此代码中可以看到的：
## corr method will give you the correlation between features
df[numeric_cols].corr()
下图显示了上一行的输出：
您可以使用heatmap更好地表示这种关系，如下所示：
plt.figure(figsize=(15, 15))
sns.heatmap(df[numeric_cols].corr(), annot=True)
下图是由前面的代码块生成的热图：
Correlations heatmap
您还可以通过以下方式过滤相关性：
plt.figure(figsize=(18, 15))
sns.heatmap(df[numeric_cols].corr()
[(df[numeric_cols].corr() >= 0.3) & (df[numeric_cols].corr() <= 0.5)],
annot=True)
此图显示了过滤后的相关性：
Filtered correlations heatmap
您还可以使用其他有用的可视化来检查统计关系，如下所示：
fig, ax = plt.subplots(3, 3, figsize = (18, 12))
for i, ax in enumerate(fig.axes):
if i < 9:
sns.regplot(x=df[numeric_cols[i]],y='Target', data=df, ax=ax)
该图显示了前面代码的以下输出：
Regression Plots
您可以看到，使用pandas使探索性数据分析相对简单。 使用pandas，您可以检查特征及其关系。
Pandas 和股票价格的定量建模
pandas最初是为在金融数据集中使用而编写的，它包含许多用于处理时间序列数据的便捷函数。 在本节中，您将看到如何使用pandas库处理股票价格序列。
您将使用quandl Python 库获取公司的财务数据。 看一下这段代码：
import quandl msft = quandl.get('WIKI/MSFT')
msft.columns
## Index(['Open', 'High', 'Low', 'Close', 'Volume', 'Ex-Dividend', 'Split Ratio', 'Adj. Open', 'Adj. High', 'Adj. Low', 'Adj. Close', 'Adj. Volume'], dtype='object')
msft.tail()
下表显示了msft.tail()的输出：
让我们通过导入以下设置来自定义绘图：
matplotlib.font_manager as font_manager font_path = '/Library/Fonts/Cochin.ttc' font_prop = font_manager.FontProperties(fname=font_path, size=24) axis_font = {'fontname':'Arial', 'size':'18'} title_font = {'fontname':'Arial', 'size':'22', 'color':'black', 'weight':'normal', 'verticalalignment':'bottom'} plt.figure(figsize=(10, 8)) plt.plot(msft['Adj. Close'], label='Adj. Close') plt.xticks(fontsize=22) plt.yticks(fontsize=22) plt.xlabel("Date", **axis_font) plt.ylabel("Adj. Close", **axis_font) plt.title("MSFT", **title_font) plt.legend(loc='upper left', prop=font_prop, numpoints=1) plt.show()
该图显示了来自先前设置的图：
您可以使用以下代码来计算每日变化：
msft['Daily Pct. Change'] = (msft['Adj. Close'] - msft['Adj. Open']) / msft['Adj. Open']
msft.tail(10)
下图显示了msft.tail(10)的输出：
尝试使用此每日收益的直方图：
plt.figure(figsize=(22, 8))
plt.hist(msft['Daily Pct. Change'], bins=100)
如下图所示，它将为您提供以下图表：
Distribution of daily returns
返回的分布具有较长的尾巴，尤其是在负侧，这是财务分析中的已知现象。 它产生的风险称为尾部风险，它与市场收益服从正态分布的假设相矛盾。 这基本上告诉您，极端事件发生的可能性比更正态分布的可能性更大。
在可视化方面，使它们具有交互性很有帮助。 为此，plotly提供了一个很好的替代当前绘图库的方法，如下所示：
import plotly.plotly as py
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot init_notebook_mode(connected=True)
from datetime import datetime
import pandas_datareader.data as web
import quandl
msft = quandl.get('WIKI/MSFT')
msft['Daily Pct. Change'] = (msft['Adj. Close'] - msft['Adj. Open']) / msft['Adj. Open']
data = [go.Scatter(x=msft.index, y=msft['Adj. Close'])]
plot(data)
我们从前面的代码中获得以下图表，如下图所示：
您可以创建开盘-最高-最低-收盘价（OHLC）的图表，其中每个日期都有 4 个不同的价格值，它们是开，高，低和关。 看一下这段代码：
charts trace = go.Ohlc(x=msft.index, open=msft['Adj. Open'], high=msft['Adj. High'], low=msft['Adj. Low'], close=msft['Adj. Close'])
data = [trace]
plot(data)
该图显示了先前代码的图：
您可以通过在图表上选择自定义范围来检查特定区域。 看一下这个图：
同样，您可以使用以下代码创建Candlestick图表：
trace = go.Candlestick(x=msft.index, open=msft['Adj. Open'], high=msft['Adj. High'], low=msft['Adj. Low'], close=msft['Adj. Close'])
data = [trace]
plot(data)
下图显示了此代码的输出：
您也可以在Candlestick图表中选择特定范围。 看一下这个图：
分布图如下：
import plotly.figure_factory as ff
fig = ff.create_distplot([msft['Daily Pct. Change'].values], ['MSFT Daily Returns'], show_hist=False)
plot(fig)
下图显示了前面代码的输出：
我们可以创建三个移动平均线，如下所示：
msft['200MA'] = msft['Adj. Close'].rolling(window=200).mean()
msft['100MA'] = msft['Adj. Close'].rolling(window=100).mean()
msft['50MA'] = msft['Adj. Close'].rolling(window=50).mean()
msft.tail(10)
下表显示了msft.tail(10)的输出：
根据切片的数据，将包括最近 2,000 天。 看一下这段代码：
trace_adjclose = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['Adj. Close'], name = "Adj. Close", line = dict(color = '#000000'), opacity = 0.8)
trace_200 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['200MA'], name = "200MA", line = dict(color = '#FF0000'), opacity = 0.8)
trace_100 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['100MA'], name = "100MA", line = dict(color = '#0000FF'), opacity = 0.8)
trace_50 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['50MA'], name = "50MA", line = dict(color = '#FF00FF'), opacity = 0.8)
data = [trace_adjclose, trace_200, trace_100, trace_50]
layout = dict( title = "MSFT Moving Averages: 200, 100, 50 days", )
fig = dict(data=data, layout=layout)
plot(fig)
下图显示了前面代码块中的图：
移动平均线用于监视金融市场的趋势。 在此示例中，有三个移动平均线，每个移动线均具有不同的周期。 您可以设置分析天数，以进行短期，中期和长期趋势监视。
当您开始使用财务时间序列时，您会很快意识到您需要基于不同期间的汇总，并且在pandas中创建这些汇总非常容易。 以下代码段将通过计算平均值每月汇总记录：
msft_monthly = msft.resample('M').mean()
msft_monthly.tail(10)
下图显示了msft_monthly.tail(10)的输出：
这是一个简单的时间序列图：
data = [go.Scatter(x=msft_monthly[-24:].index, y = msft_monthly[-24:]['Adj. Close'])]
plot(data)
如下图所示，这将为您提供以下图表：
在检查要素之间的关系时，可以使用我们在前面的示例中已经看到的相关矩阵。 在时间序列中，从业者对自相关感兴趣，自相关显示了时间序列与其自身的相关性。 例如，理想情况下，您希望时间序列中出现周期性的峰值，以显示您的季节性。 通过使用以下代码，让我们看看每日百分比变化是否有任何明显的峰值：
plt.figure(figsize=(22, 14)) pd.plotting.autocorrelation_plot(msft_monthly['Daily Pct. Change'])
我们从前面的代码中获得以下图表，如下图所示：
Monthly autocorrelation plot
在这个系列中没有有意义的显着滞后，但是如果您使用宏观经济变量（例如 GDP，通货膨胀率和失业水平）进行尝试，则可能会看到显着的季度或年度峰值。
SciPy 和 Scikit-learn
Scikit-learn 是用于机器学习的 SciKit 库之一，它建立在 SciPy 之上。 您可以使用它执行回归分析，就像在前几章中使用 scikit-learn 库所做的那样。 看一下这段代码：
from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score
diabetes = datasets.load_diabetes()
linreg = linear_model.LinearRegression()
linreg.fit(diabetes.data, diabetes.target)
## You can inspect the results by looking at evaluation metrics
print('Coeff.: n', linreg.coef_)
print("MSE: {}".format(mean_squared_error(diabetes.target, linreg.predict(diabetes.data)))) print('Variance Score: {}'.format(r2_score(diabetes.target, linreg.predict(diabetes.data))))
scikit-learn 和住房数据的 K 均值聚类
在本节中，我们将使用 scikit-learn 的 K 均值算法对房屋数据进行聚类，如下所示：
from sklearn.cluster import KMeans from sklearn.datasets import load_boston boston = load_boston()
## As previously, you have implemented the KMeans from scracth and in this example, you use sklearns API k_means = KMeans(n_clusters=3) # Training k_means.fit(boston.data)
KMeans(algorithm='auto', copy_x=True, init='K 均值++', max_iter=300, n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto', random_state=None, tol=0.0001, verbose=0)
print(k_means.labels_)
前面代码的输出如下：
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 0 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 0 0 0 0 0 0 0 0 0 0 0 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2
2 0 2 2 2 2 0 2 2 2 0 0 0 0 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1]
您可以使用以下代码行检查群集中心：
print(k_means.cluster_centers_)
这是控制台的输出：
[[ 1.49558803e+01 -5.32907052e-15 1.79268421e+01 2.63157895e-02
6.73710526e-01 6.06550000e+00 8.99052632e+01 1.99442895e+00
2.25000000e+01 6.44736842e+02 1.99289474e+01 5.77863158e+01
2.04486842e+01]
[ 3.74992678e-01 1.57103825e+01 8.35953552e+00 7.10382514e-02
5.09862568e-01 6.39165301e+00 6.04133880e+01 4.46074481e+00
4.45081967e+00 3.11232240e+02 1.78177596e+01 3.83489809e+02
1.03886612e+01]
[ 1.09105113e+01 5.32907052e-15 1.85725490e+01 7.84313725e-02
6.71225490e-01 5.98226471e+00 8.99137255e+01 2.07716373e+00
2.30196078e+01 6.68205882e+02 2.01950980e+01 3.71803039e+02
1.78740196e+01]]
关于聚类算法的评估，通常将使用诸如轮廓分析或弯头方法之类的技术来评估聚类的质量并确定正确的超参数（例如 K 均值）。 使用 scikit-learn 提供的简单 API，您还将发现这种分析易于执行。 强烈建议您通过实践从这些示例中学到的知识为基础，这将提高您的知识和技能。
总结
在本章中，您使用各种示例（主要用于机器学习任务）练习了 NumPy，SciPy，Pandas 和 scikit-learn。 使用 Python 数据科学库时，通常有不止一种执行给定任务的方法，而且通常有助于了解不止一种方法。
您可以使用替代方法以获得更好的实现，也可以出于比较的目的。 在为给定任务尝试不同的方法时，您可能会找到不同的选项，这些选项将允许您进一步自定义实现，或者只是观察到一些性能改进。
本章的目的是向您展示这些不同的选项，以及 Python 语言由于其丰富的分析库生态系统而具有的灵活性。 在下一章中，您将了解有关 NumPy 内部的更多信息，例如 numpy 如何管理数据结构和内存，代码概要分析以及有效编程的技巧。
七、高级 NumPy
许多库都具有易于使用的 API。 您需要做的就是调用提供的 API 函数，库将为您处理其余的函数。 幕后发生的事情与您无关，您只对输出感兴趣。 在大多数情况下，这很好，但是，至少要了解所使用库的基本内部结构很重要。 了解内部结构将帮助您掌握代码的最新动态以及开发应用时应避免的危险信号。
在本章中，将回顾 NumPy 的内部结构，例如 NumPy 的类型层次结构和内存使用情况。 在本章的最后，您将了解用于逐行检查程序的代码配置文件。
NumPy 内部
如您在前几章中所见，NumPy 数组使数值计算变得高效，其 API 直观且易于使用。 NumPy 数组也是其他科学图书馆的核心，因为其中许多库都是基于 NumPy 数组构建的。
为了编写更好，更高效的代码，您需要了解数据处理的内部。 NumPy 数组及其元数据位于数据缓冲区中，该缓冲区是带有某些数据项的专用内存块。
NumPy 如何管理内存？
初始化 NumPy 数组后，其元数据和数据将存储在随机存取存储器（RAM）中分配的存储位置上。
import numpy as np
array_x = np.array([100.12, 120.23, 130.91])
首先，Python 是一种动态类型化的语言； 不需要显式声明变量类型，例如int或double。 可以推断变量类型，在这种情况下，您希望array_x的数据类型为np.float64：
print(array_x.dtype)
float64
使用numpy库而不是 Python 的优势在于numpy支持许多不同的数值数据类型，例如bool_，int_，intc，intp，int8，int16，int32，int64，uint8，uint16，uint32，uint64，float_，float16，float32，float64，complex_，complex64和complex128。
您可以通过检查sctypes查看这些类型：
np.sctypes
{'complex': [numpy.complex64, numpy.complex128, numpy.complex256],
'float': [numpy.float16, numpy.float32, numpy.float64, numpy.float128],
'int': [numpy.int8, numpy.int16, numpy.int32, numpy.int64],
'others': [bool, object, bytes, str, numpy.void],
'uint': [numpy.uint8, numpy.uint16, numpy.uint32, numpy.uint64]}
下图显示了数据类型树：
您可以通过调用mro方法来检查诸如np.float64之类的数据类型的父类：
np.float64.mro()
[numpy.float64,
numpy.floating,
numpy.inexact,
numpy.number,
numpy.generic,
float,
object]
和np.int64的父类：
np.int64.mro()
[numpy.int64,
numpy.signedinteger,
numpy.integer,
numpy.number,
numpy.generic,
object]
mro方法代表“方法解析顺序”。 为了更好地理解继承，应该首先了解继承的概念。 在可以使用面向对象范例的编程语言中，可以将对象的属性和方法基于之前创建的另一个对象，这就是继承。 在前面的示例中，np.int64保留了np.signedinteger及其之后的属性和行为。
让我们看一个简单的例子：
class First:
def firstmethod(self):
print("Call from First Class, first method.")
class Second:
def secondmethod(self):
print("Call from Second Class, second method.")
class Third(First, Second):
def thirdmethod(self):
print("Call from Third Class, third method.")
这里，有 3 个类别，而First和Second类别是独立的，Third类别是从First和Second继承的。 您可以创建Third类的实例，并使用dir方法检查其内容：
myclass = Third()
dir(myclass)
[...
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'firstmethod',
'secondmethod',
'thirdmethod']
dir表示myclass的方法中有firstmethod，secondmethod和thirdmethod。
您可以调用这些方法，如下所示：
myclass.firstmethod()
myclass.secondmethod()
myclass.thirdmethod()
## Call from First Class, first method.
## Call from Second Class, second method.
## Call from Third Class, third method.
现在，让我们向Second类添加firstmethod，看看会发生什么：
class First:
def firstmethod(self):
print("Call from First Class, first method.")
class Second:
def firstmethod(self):
print("Call from Second Class, first method.")
def secondmethod(self):
print("Call from Second Class, second method.")
class Third(First, Second):
def thirdmethod(self):
print("Call from Third Class, third method.")
像以前一样检查方法输出：
myclass = Third()
myclass.firstmethod()
myclass.secondmethod()
myclass.thirdmethod()
## Call from First Class, first method.
## Call from Second Class, second method.
## Call from Third Class, third method.
如您所见，已添加到Second类的方法无效，因为Third类的实例从First类继承了该实例。
您可以按以下方式检查类的mro：
Third.__mro__
这将为您提供以下输出：
(__main__.Third, __main__.First, __main__.Second, object)
这是使用继承机制时解析属性和方法的方式，并且现在您应该或多或少地了解mro的工作方式。 现在，您可以再次查看我们之前拥有的 numpy 数据类型的mro示例。
您可以使用nbytes来查看存储数据类型所需的内存。
首先，让我们看看单个float64的大小：
np.float64(100.12).nbytes
8
np.str_('n').nbytes
4
np.str_('numpy').nbytes
20
array_x具有 3 个float64，其大小将是元素数乘以商品大小，即24，如以下代码段所示：
np.float64(array_x).nbytes
24
例如，如果您需要较低的计算精度，则可以使用np.float32，它将占用float64占用的一半内存：
array_x2 = array_x.astype(np.float32)
array_x2
array([100.12, 120.23, 130.91], dtype=float32)
np.float32(array_x2).nbytes
12
简单来说，8 个字节的内存将保存 1 float64或 2 float32。
Python 的动态性质引入了一种处理数据类型的新方法，因为 Python 应该包含有关其存储的数据的更多信息。 虽然典型的 C 变量将具有有关内存位置的信息，但 Python 变量应具有存储为 C 结构的信息，该结构包含引用计数，对象的类型，对象的大小以及变量本身。
这是提供灵活的环境来处理不同数据类型所必需的。 如果诸如列表之类的数据结构可以容纳不同类型的对象，这是由于该信息对于列表中的每个元素的存储。
但是，由于 NumPy 数组中的数据类型是固定的，由于使用了连续的内存块，因此内存使用效率可能更高。
您可以通过检查 NumPy 数组的__array_interface__属性来查看地址和其他相关信息。
编写此接口是为了允许开发人员共享数组内存和信息：
array_x.__array_interface__
{'data': (140378873611440, False),
'descr': [('', '<f8')],
'shape': (3,),
'strides': None,
'typestr': '<f8',
'version': 3}
__array_interface__是具有 6 个键的 python 字典：
shape的工作方式类似于 NumPy 数组或pandas数据帧的常规shape属性。 它显示每个大小的大小。 由于array_x具有1大小和3元素，因此它是具有3大小的元组。typestr具有3值，第一个显示字节顺序，第二个显示字符代码，其余字符显示字节数。 在此示例中，其值为'<f8'，这表示字节顺序为低位字节序，字符代码为浮点，并且使用的字节数为 8。descr可能会提供有关内存布局的更多详细信息。 默认值为[('', typestr)]。data显示数据的存储位置。 这是一个元组，其中第一个元素显示 NumPy 数组的存储块地址，第二个元素是指示其是否为只读的标志。 在此示例中，内存块地址为140378873611440，它不是只读的。strides指示给定的数组是否为 C 样式的连续内存缓冲区。 在此示例中，None 表示这是 C 样式的连续数组。 否则，它将包含跨步元组，以了解跳转到给定维度中的下一个数组元素所要跳转的位置。 步幅是重要的属性，因为当您使用不同的切片（例如X[::4]）时，步幅将引导数组视图。version表示在此示例中版本号为 3。
以下片段显示了一个简单的示例：
import numpy as np
X = np.array([1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27])
X[::4]
## array([ 1, 1, 11, 14, 30])
这一点很重要，因为当您使用基于现有ndarrays的切片创建新的ndarrays时，可能会降低性能。 让我们看一个简单的例子； 以下代码段创建了 3D ndarray：
nd_1 = np.random.randn(4, 6, 8)
nd_1
## array([[[ 0.64900179, -0.00494884, -0.97565618, -0.78851039],
[ 0.05165607, 0.068948 , 1.54542042, 1.68553396],
[-0.80311258, 0.95298682, -0.85879725, 0.67537715]],
[[ 0.24014811, -0.41894241, -0.00855571, 0.43483418],
[ 0.43001636, -0.75432657, 1.16955535, -0.42471807],
[ 0.6781286 , -1.87876591, 1.02969921, 0.43215107]]])
您可以对其进行切片并创建另一个数组：
nd_2 = nd_1[::, ::2, ::2]
将会选择：
- 首先，第一维的所有项目
- 然后，第二维的每两个项目
- 然后，第三维的每两个项目
它将具有以下数组：
print(nd_2)
[[[ 0.64900179 -0.97565618]
[-0.80311258 -0.85879725]]
[[ 0.24014811 -0.00855571]
[ 0.6781286 1.02969921]]]
您可以看到nd_1和nd_2的内存地址相同：
nd_1.__array_interface__
{'data': (140547049888960, False),
'descr': [('', '<f8')],
'shape': (2, 3, 4),
'strides': None,
'typestr': '<f8',
'version': 3}
nd_2.__array_interface__
{'data': (140547049888960, False),
'descr': [('', '<f8')],
'shape': (2, 2, 2),
'strides': (96, 64, 16),
'typestr': '<f8',
'version': 3}
nd_2大步前进，了解如何沿nd_1数组的不同维度移动。
为了强调这些跨步在数值计算中的作用，下面的示例将为数组维和切片使用更大的大小：
nd_1 = np.random.randn(400, 600)
nd_2 = np.random.randn(400, 600*20)[::, ::20]
nd_1和nd_2具有相同的大小：
print(nd_1.shape, nd_2.shape)
(400, 600) (400, 600)
您可以测量用于计算nd_1和nd_2的数组元素的累积乘积的时间：
%%timeit
np.cumprod(nd_1)
## 802 µs ± 20.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
np.cumprod(nd_2)
## 12 ms ± 71.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
两次操作之间存在明显的时间差； 这是为什么？ 如您所料，nd_2中的步幅过大会导致此问题：
nd_1.__array_interface__
{'data': (4569473024, False),
'descr': [('', '<f8')],
'shape': (400, 600),
'strides': None,
'typestr': '<f8',
'version': 3}
nd_2.__array_interface__
{'data': (4603252736, False),
'descr': [('', '<f8')],
'shape': (400, 600),
'strides': (96000, 160),
'typestr': '<f8',
'version': 3}
从存储器向 CPU 读取数据时，nd_2中存在跨步会导致跳转到不同的存储器位置。 如果将数组元素顺序存储为连续的内存块，那么从时间测量来看，此操作会更快。 步伐越小越好，可以更好地利用 CPU 缓存来提高性能。
有一些变通办法可以缓解与 CPU 缓存相关的问题。 您可以使用的一种库是numexpr库，它是 NumPy 的快速数值表达式求值器。 库使内存使用效率更高，并且还可以使多线程编程受益，以充分利用可用的内核。 您也可以将其与 Intel 的 VML 结合使用以进行进一步的改进。
在下面的示例中，让我们看看它是否对nd_2有帮助：
import numexpr as ne
%%timeit
2 * nd_2 + 48
## 4 ms ± 10.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
ne.evaluate("2 * nd_2 + 48")
## 843 µs ± 8.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
您应该使用其他示例进行尝试，以查看性能提升。
如果从头开始索引数组到某个元素，您将看到它具有相同的内存地址：
array_x[:2].__array_interface__['data'][0]
## 140378873611440
array_x[:2].__array_interface__['data'][0] == array_x.__array_interface__['data'][0]
## True
但是，如果您在0以外的位置开始索引，则将为您提供不同的内存地址：
array_x[1:].__array_interface__['data'][0]
## 140378873611448
array_x[1:].__array_interface__['data'][0] == array_x.__array_interface__['data'][0]
## False
ndarray的另一个属性称为标志，它提供有关给定 NumPy 数组的内存布局的信息：
array_f = np.array([[100.12, 120.23, 130.91], [90.45, 110.32, 120.32]])
print(array_f)
## [[100.12 120.23 130.91]
## [ 90.45 110.32 120.32]]
array_f.flags
## C_CONTIGUOUS : True
## F_CONTIGUOUS : False
## OWNDATA : True
## WRITEABLE : True
## ALIGNED : True
## WRITEBACKIFCOPY : False
## UPDATEIFCOPY : False
您可以使用类似于字典的符号或小写的属性名称来获得单个标志：
array_f.flags['C_CONTIGUOUS']
## True
array_f.flags.c_contiguous
## True
让我们看一下每个属性：
C_CONTIGUOUS：C 样式的连续内存的单个块F_CONTIGUOUS：连续内存的单个块，Fortran 风格
您的数据可以使用不同的布局存储在内存中。 这里有 2 种不同的内存布局要考虑：对应于C_CONTIGUOUS的行主要顺序和对应于F_CONTIGUOUS的列主要顺序。
在该示例中，array_f是二维的，array_f的行项目存储在相邻的存储位置中。 类似地，在F_CONTIGUOUS情况下，每列的值存储在相邻的存储位置中。
某些numpy函数将使用参数order将此顺序指示为'C'或'F'。 以下示例显示了具有不同顺序的reshape函数：
np.reshape(array_f, (3, 2), order='C')
## array([[100.12, 120.23],
## [130.91, 90.45],
## [110.32, 120.32]])
np.reshape(array_f, (3, 2), order='F')
## array([[100.12, 110.32],
## [ 90.45, 130.91],
## [120.23, 120.32]])
其余的：
OWNDATA：数组是否与另一个对象共享其内存块或拥有所有权WRITEABLE：False表示它是只读的； 否则可以将该区域写入。ALIGNED：数据是否针对硬件对齐WRITEBACKIFCOPY：该数组是否是另一个数组的副本UPDATEIFCOPY：（不建议使用WRITEBACKIFCOPY）
了解内存管理很重要，因为它会影响性能。 根据您执行计算的方式，计算速度会有所不同。 您可能没有意识到某些计算涉及现有数组的隐式副本，这会减慢计算速度。
以下代码块显示了两个示例，其中第一个不需要复制，而第二个具有隐式复制操作：
shape = (400,400,400)
array_x = np.random.random_sample(shape)
import cProfile
import re
cProfile.run('array_x *= 2')
## 3 function calls in 0.065 seconds
## Ordered by: standard name
## ncalls tottime percall cumtime percall filename:lineno(function)
## 1 0.065 0.065 0.065 0.065 <string>:1(<module>)
## 1 0.000 0.000 0.065 0.065 {built-in method builtins.exec}
## 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
import cProfile
import re
cProfile.run('array_y = array_x * 2')
## 3 function calls in 0.318 seconds
## Ordered by: standard name
## ncalls tottime percall cumtime percall filename:lineno(function)
## 1 0.318 0.318 0.318 0.318 <string>:1(<module>)
## 1 0.000 0.000 0.318 0.318 {built-in method builtins.exec}
## 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
首次运行比第二次慢 5 倍。 您需要了解隐式复制操作，并熟悉它在哪种情况下发生。 重塑数组不需要隐式复制，除非它是矩阵转置。
许多数组操作会返回一个新数组以获取结果。 此行为是预期的，但会破坏迭代任务的性能，在迭代任务中，您可能会有数百万或数十亿次迭代。 某些numpy函数具有out参数，该参数创建输出数组，并使用其写入迭代结果。 通过这种方式，您的程序可以更好地管理内存，并且需要更少的时间：
shape_x = (8000,3000)
array_x = np.random.random_sample(shape_x)
%%timeit
np.cumprod(array_x)
## 176 ms ± 2.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
output_array的类型和大小应与操作的预期输出相同：
output_array = np.zeros(array_x.shape[0] * array_x.shape[1])
%%timeit
np.cumprod(array_x, out=output_array)
## 86.4 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
对 NumPy 代码进行性能分析以了解性能
有几个有用的库可以监视给定 python 脚本的性能指标。 您已经看到cProfile库的用法。 本节将介绍vprof，它是可视分析器库。
它将为您提供给定 python 程序的运行时统计信息和内存利用率。
第 5 章“使用 NumPy 聚类批发分销商的客户”的一维特征将在此处使用，以下代码段应保存到名为to_be_profiled.py的文件中：
import numpy as np
X = np.array([1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27])
n_clusters = 3
def Kmeans_1D(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([np.abs(point - c_centers) for point in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([np.abs(point - new_c_centers) for point in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.all(c_centers == new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
c_centers, labels = Kmeans_1D(X, 3)
print(c_centers, labels)
保存此文件后，您可以使用命令行开始对其进行性能分析。
可以通过 4 种不同的方式配置vprof来获取：
- CPU 火焰图（
vprof -c c to_be_profiled.py） - 内置的探查器统计信息（
vprof -c p to_be_profiled.py） - 运行程序中的行后，Cpython 垃圾收集器跟踪的对象和进程内存的内存图（
vprof -c m to_be_profiled.py） - 用
runtime编码heatmap并为每条执行的行计数（vprof -c h to_be_profiled.py）
这 4 种配置可用于单个源文件或包。 让我们看一下p，m和h配置的输出：
探查器的配置：
$ vprof -c p to_be_profiled.py
Running Profiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
一个新的选项卡将在您的浏览器中打开，并显示以下输出：
Time spent for each call
您可以看到文件名，函数名，行号和每次调用所花费的时间。
内存使用情况统计信息的配置：
$ vprof -c m to_be_profiled.py
Running MemoryProfiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
新标签页将在您的浏览器中打开，并显示以下输出：
内存使用情况
在左侧，您可以看到内存中的对象，并且图表显示了随着执行的行数增加，内存使用量（以兆字节为单位）。 如果将鼠标悬停在图表上，则每个点将具有以下信息：
- 执行的行
- 行号
- 函数名称
- 文件名称
- 内存使用情况
例如，to_be_profiled.py文件中的第 27 行是计算deltas的下一行：
deltas = np.array([np.abs(point - new_c_centers) for point in X])
它执行了很多次，因为如果您检查图表，这是列表理解。
代码heatmap的配置：
$ vprof -c h to_be_profiled.py
Running CodeHeatmapProfiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
“新”标签将在您的浏览器中打开，并显示以下输出：
Heatmap for the lines executed
在左侧，有一个已检查模块的列表，在这种情况下，只有一个文件要检查。
右边是程序中每行的热图。 如果将鼠标悬停在行上，它将为您提供以下信息：
- 花费的时间
- 总运行时间
- 百分比
- 运行计数
如果将鼠标悬停在27行上，它将为您提供以下统计信息：
总结
在进行科学操作时，了解 NumPy 的内部结构至关重要。 效率是关键，因为许多科学计算都是计算和内存密集型的。 因此，如果您的代码编写效率不高，则计算所需的时间将远远超过所需的时间，这将损害您的研发时间。
在本章中，您已经了解了 NumPy 库的一些内部和性能方面，还了解了vprof库，该库可帮助您检查 python 程序的性能。
代码概要分析将极大地帮助您逐行检查程序，并且如您先前所见，查看相同数据的方式也有所不同。 确定了程序中最苛刻的部分后，便可以开始搜索更有效的方法或实现以提高性能并节省更多时间。
在下一章中，我们将概述高性能，低级的数值计算库。 NumPy 可以使用这些实现来获得可观的性能提升。
八、高性能数值计算库概述
在科学计算应用中可以执行许多数值运算，并且未经优化的代码或库实现会导致严重的性能瓶颈。
NumPy 库通过更有效地使用其内存布局来帮助提高 Python 程序的性能。
在实际应用中，最常用的数学分支之一是线性代数。 线性代数用于计算机图形学，密码学，计量经济学，机器学习，深度学习和自然语言处理，仅举几个例子。 具有高效的矩阵和向量运算至关重要。
高性能，低级框架（例如 BLAS，LAPACK 和 ATLAS）是 Netlib 库的一部分，用于密集的线性代数运算；其他框架（例如 Intel MKL）也可以在其中使用您的程序。 这些库在计算中具有很高的性能和准确性。 您可以通过其他高级编程语言（例如 Python 或 C++）调用它们来使用这些库。
当 NumPy 与不同的 BLAS 库链接时，您可以观察到性能差异而无需更改代码，因此了解哪种链接可以更好地提高性能非常重要。
让我们看一下这些低级库。
BLAS 和 LAPACK
BLAS 代表基本线性代数子程序，并且是处理线性代数运算的低级例程的标准。 低级例程包括向量和矩阵加/乘，线性组合等操作。 BLAS 为线性代数运算提供了三种不同的级别：
- BLAS1：标量向量和向量向量运算
- BLAS2：矩阵向量运算
- BLAS3：矩阵矩阵运算
LAPACK 代表线性代数软件包，并包含更高级的操作。 LAPACK 提供了用于矩阵分解（例如 LU，Cholesky 和 QR）以及解决特征值问题的例程。 LAPACK 主要取决于 BLAS 例程。
ATLAS
有许多优化的 BLAS 实现。 ATLAS 代表自动调谐线性代数软件，并且是与平台无关的项目，可以生成优化的 BLAS 实现。
英特尔 MKL
英特尔 MKL 为英特尔处理器优化了 BLAS。 改进了例程和函数，例如 1 级，2 级和 3 级 BLAS，LAPACK 例程，求解器，FFT 函数，其他数学和统计函数。 这些改进的例程和函数得益于共享内存多处理等改进，它们可用于在发行版（如 Anaconda 发行版）中加速科学 python 库（例如 NumPy 和 SciPy）。 如果您查看其发行说明， 可以看到每个发行版都进行了一些重要的改进，例如 LAPACK 函数的性能得到了提高。
OpenBLAS
OpenBLAS 是另一个优化的 BLAS 库，它为不同的配置提供了 BLAS3 级的优化。 作者报告说，与 BLAS 相比，性能增强和改进可与英特尔 MKL 的性能相媲美。
使用 AWS EC2 和底层库配置 NumPy
- 登录到 AWS。 如果您没有帐户，请创建一个：
-
选择 EC2 。
-
单击启动实例：
- 选择
Ubuntu Server 16.04 LTS (HVM), SSD Volume Type - ami-db710fa3：
- 选择
t2.micro实例类型：
- 点击启动：
- 单击启动。
- 选择创建一个新的密钥对：
- 给它命名，然后单击启动实例。 它需要一些时间才能运行：
- 一旦其状态为
running，点击实例 ID ，在这种情况下为i-00ccaeca61a24e042。 然后选择实例并单击Connect：
- 然后它将向您显示以下窗口，其中包含一些有用的信息：
- 打开终端，然后导航到保存所生成密钥的文件夹。 在此示例中，键名称为
aws_oregon。 运行以下命令：
$ chmod 400 aws_oregon.pem
- 然后，在上一个窗口的示例部分中复制该行并运行它：
$ ssh -i "aws_oregon.pem" ubuntu@ec2-34-219-121-1.us-west-2.compute.amazonaws.com
- 在第一个问题的答案中输入
yes将其添加到已知主机中，它将连接到您的实例：
您需要做的第一件事是通过运行以下命令来更新和升级预安装的软件包：
sudo apt-get update
sudo apt-get upgrade
sudo短语为您提供了更新和升级的必要权利，因为软件包的更改可能会对系统产生负面影响，并非所有人都可以对其进行授权。 您可以将apt-get视为 Ubuntu 的软件包管理器。
您可以创建许多虚拟环境，并链接到不同的低级库，但是，每次使用新的低级库配置 NumPy 时，您都将从一个新的预配实例开始。 这将为您提供有关配置过程的想法，以后使您可以轻松地设置虚拟环境。
安装 BLAS 和 LAPACK
为了设置您的开发环境，您需要在运行以下命令后安装所需的软件包，例如编译器，库和其他必要的部分，
$ sudo apt-get update $ sudo apt-get upgrade
对于此配置，很幸运，因为您可以运行以下命令来安装 Python 的 SciPy 软件包，它将安装所有必需的软件包，包括 NumPy，基本线性代数子程序（libblas3）和线性代数软件包（liblapack3）：
$ sudo apt-get install python3-scipy
控制台输出：
- 输入
Y并按Enter继续。 安装完成后，运行以下命令以打开python3解释器：
$ python3
启动 Python 控制台：
- 导入
numpy并使用show_config方法查看 NumPy 的配置：
控制台输出：
- 由于在安装 NumPy 时可以使用 BLAS 和 LAPACK 库，因此它使用它们来构建库。 您可以在
lapack_info和blas_info中看到它们； 否则，您将不会在输出中看到它们，如以下屏幕截图所示：
- 如果您使用的是 macOS，则可以使用 Accelerate/vecLib 框架。 以下命令将输出 macOS 系统的加速器选项：
安装 OpenBLAS
OpenBLAS 的步骤略有不同，如下所示：
- 在先前的配置中运行以下命令：
$ sudo apt-get update $ sudo apt-get upgrade
- 您需要通过运行以下命令来安装
build-essential，其中包括make命令和其他必要的库：
$ sudo apt-get install build-essential libc6 gcc gfortran
- 创建一个名为
openblas_setup.sh的文件，然后粘贴以下内容。 如果您搜索 GitHub，则可以找到不同的设置脚本，并且可以尝试一种满足您需要的脚本：
##!/bin/bash
set -e
pushd /root
git clone https://github.com/xianyi/OpenBLAS.git
pushd /root/OpenBLAS
make clean
make -j4
rm -rf /root/openblas-install
make install PREFIX=/root/openblas-install
popd
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/libblas.so
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/libblas.so.3
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/liblapack.so.3
- 保存此文件并运行以下命令：
$ chmod +777 openblas_setup.sh
$ sudo ./openblas_setup.sh
- 安装完成后，您可以按以下方式安装 numpy 和 scipy：
$ sudo apt-get install python3-pip
$ pip3 install numpy
$ pip3 install scipy
- 现在，您可以像以前一样检查 NumPy 配置：
安装英特尔 MKL
为了针对英特尔 MKL 构建 NumPy 和 SciPy，请按照以下说明进行操作
- 运行以下命令：
$ sudo apt-get update
$ sudo apt-get upgrade
- 您需要安装 Anaconda 发行版，因为 Anaconda 的安装随 Intel MKL 一起提供。 首先使用以下命令下载 Anaconda：
$ wget https://repo.continuum.io/archive/Anaconda3-5.2.0-Linux-x86_64.sh
- 安装完成后，将
cd插入anaconda3/bin并运行python：
$ cd anaconda3/bin
$ ./python
- 您可以像以前一样检查
numpy配置：
安装 ATLAS
为了针对 ATLAS 构建 NumPy，请按照以下说明进行操作
- 运行以下命令：
$ sudo apt-get update
$ sudo apt-get upgrade
- 您需要通过运行以下命令来安装
build-essential，其中包括make命令和其他必要的库：
$ sudo apt-get install build-essential libc6 gcc gfortran
- 然后，您需要安装
atlas：
$ sudo apt-get install libatlas-base-dev
- 您现在可以按以下方式安装
pip和numpy：
$ sudo apt-get install python3-pip
$ pip3 install --no-cache-dir Cython
$ git clone https://github.com/numpy/numpy.git
$ cd numpy
$ cp site.cfg.example site.cfg
$ vi site.cfg
在site.cfg内，您应注释掉地图集行，并将其设置为地图集安装，如下所示：
[atlas]
library_dirs = /usr/local/atlas/lib
include_dirs = /usr/local/atlas/include
然后运行：
$ sudo python3 setup.py install
- 安装完成后，安装
scipy：
$ pip3 install scipy
然后返回到您的主目录，启动python解释器并检查numpy配置，这将为您提供以下输出：
>>> import numpy as np
>>> np.show_config()
atlas_blas_info:
include_dirs = ['/usr/include/atlas']
language = c
library_dirs = ['/usr/lib/atlas-base']
define_macros = [('HAVE_CBLAS', None), ('ATLAS_INFO', '"\\"3.10.2\\""')]
libraries = ['f77blas', 'cblas', 'atlas', 'f77blas', 'cblas']
...
您已经介绍了上述所有低级库的配置。 是时候了解基准测试的计算密集型任务了。
用于基准测试的计算密集型任务
现在，您将能够使用不同的配置（例如是否使用 BLAS/LAPACK，OpenBLAS，ATLAS 和 Intel MKL）对 NumPy 性能进行基准测试。 让我们回顾一下要为基准计算的内容。
矩阵分解
矩阵分解或分解方法涉及计算矩阵的组成部分，以便可以使用它们简化要求更高的矩阵操作。 在实践中，这意味着将您拥有的矩阵分解为多个矩阵，这样，当您计算这些较小矩阵的乘积时，您将获得原始矩阵。 矩阵分解方法的一些示例是奇异值分解（SVD），特征值分解，Cholesky 分解，下上（LU）和 QR 分解。
奇异值分解
SVD 是线性代数中最有用的工具之一。 Beltrami 和 Jordan 发表了有关其使用的几篇论文。 SVD 用于各种应用，例如计算机视觉和信号处理。
如果您具有正方形或矩形矩阵（M），则可以将其分解为矩阵（U），矩阵（V）（计算中使用矩阵转置）和奇异值（d）。
您的最终公式将如下所示：
以下是奇异值分解的说明：
一种简单的数据精简方法是排除该公式中d小到可以忽略不计的部分。
让我们看看如何使用numpy来实现：
import numpy as np
M = np.random.randint(low=0, high=20, size=20).reshape(4,5)
print(M)
## Output
[[18 15 11 13 19]
[ 1 6 8 13 18]
[ 9 7 15 13 10]
[17 15 12 14 12]]
U, d, VT = np.linalg.svd(M)
print("U:n {}n".format(U))
print("d:n {}n".format(d))
print("VT:n {}".format(VT))
U:
[[-0.60773852 -0.22318957 0.5276743 -0.54990921]
[-0.38123886 0.86738201 0.19333528 0.25480749]
[-0.42657252 0.10181457 -0.82343563 -0.36003255]
[-0.55076919 -0.43297652 -0.07832665 0.70925987]]
d:
[56.31276456 13.15721839 8.08763849 2.51997135]
VT:
[[-0.43547429 -0.40223663 -0.40386674 -0.46371223 -0.52002929] [-0.72920427 -0.29835313 0.06197899 0.27638212 0.54682545]
[ 0.11733943 0.26412864 -0.73449806 -0.30022507 0.53557916]
[-0.32795351 0.55511623 -0.3571117 0.56067806 -0.3773643 ]
[-0.39661218 0.60932187 0.40747282 -0.55144258 0.03609177]]
## Setting full_matrices to false gives you reduced form where small values close to zero are excluded
U, d, VT = np.linalg.svd(M, full_matrices=False)
print("U:n {}n".format(U))
print("d:n {}n".format(d))
print("VT:n {}".format(VT))
## Output
U:
[[-0.60773852 -0.22318957 0.5276743 -0.54990921]
[-0.38123886 0.86738201 0.19333528 0.25480749]
[-0.42657252 0.10181457 -0.82343563 -0.36003255]
[-0.55076919 -0.43297652 -0.07832665 0.70925987]]
d:
[56.31276456 13.15721839 8.08763849 2.51997135]
VT:
[[-0.43547429 -0.40223663 -0.40386674 -0.46371223 -0.52002929]
[-0.72920427 -0.29835313 0.06197899 0.27638212 0.54682545]
[ 0.11733943 0.26412864 -0.73449806 -0.30022507 0.53557916]
[-0.32795351 0.55511623 -0.3571117 0.56067806 -0.3773643 ]]
Cholesky 分解
如果您有一个正方形矩阵，也可以应用 Cholesky 分解，将一个矩阵（M）分解为两个三角形矩阵（U和U^T）。 Cholesky 分解可帮助您简化计算复杂性。 可以将其总结为以下公式：
M = U^T U
以下是 Cholesky 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from scipy.linalg import cholesky
M = np.array([[1, 3, 4],
[2, 13, 15],
[5, 31, 33]])
print(M)
## Output
[[ 1 3 4]
[ 2 13 15]
[ 5 31 33]]
L = cholesky(M)
print(L)
## Output
[[1\. 3\. 4\. ]
[0\. 2\. 1.5 ]
[0\. 0\. 3.84057287]]
L.T.dot(L)
## Output
array([[ 1., 3., 4.],
[ 3., 13., 15.],
[ 4., 15., 33.]])
LU 分解
与 Cholesky 分解类似，LU 分解将矩阵（M）分解为下（L）和上（U）三角矩阵。 这也有助于我们简化计算密集型代数。 可以将其总结为以下公式：
M = LU
下面是 LU 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from scipy.linalg import lu
M = np.random.randint(low=0, high=20, size=25).reshape(5,5)
print(M)
## Output
[[18 12 14 15 2]
[ 4 2 12 18 3]
[ 9 19 5 16 8]
[15 19 6 16 11]
[ 1 19 2 18 17]]
P, L, U = lu(M)
print("P:n {}n".format(P))
print("L:n {}n".format(L))
print("U:n {}".format(U))
## Output
P:
[[1\. 0\. 0\. 0\. 0.]
[0\. 0\. 1\. 0\. 0.]
[0\. 0\. 0\. 0\. 1.]
[0\. 0\. 0\. 1\. 0.]
[0\. 1\. 0\. 0\. 0.]]
L:
[[ 1\. 0\. 0\. 0\. 0\. ]
[ 0.05555556 1\. 0\. 0\. 0\. ]
[ 0.22222222 -0.03636364 1\. 0\. 0\. ]
[ 0.83333333 0.49090909 -0.70149254 1\. 0\. ]
[ 0.5 0.70909091 -0.32089552 0.21279832 1\. ]]
U:
[[18\. 12\. 14\. 15\. 2\. ]
[ 0\. 18.33333333 1.22222222 17.16666667 16.88888889]
[ 0\. 0\. 8.93333333 15.29090909 3.16969697]
[ 0\. 0\. 0\. 5.79918589 3.26594301]
[ 0\. 0\. 0\. 0\. -4.65360318]]
P.dot(L).dot(U)
## Output
array([[18., 12., 14., 15., 2.],
[ 4., 2., 12., 18., 3.],
[ 9., 19., 5., 16., 8.],
[15., 19., 6., 16., 11.],
[ 1., 19., 2., 18., 17.]])
特征值分解
特征值分解也是一种适用于平方矩阵的分解技术。 使用特征值分解分解方阵（M）时，将得到三个矩阵。 这些矩阵之一（Q）在列中包含特征向量，另一个矩阵（L）在对角线中包含特征值，最后一个矩阵是特征向量矩阵（Q^(-1)）。
可以将其总结为以下公式：
M = QVQ^(-1)
特征值分解将为您提供矩阵的特征值和特征向量。
下面是特征值分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from numpy.linalg import eig
M = np.random.randint(low=0, high=20, size=25).reshape(5,5)
print(M)
## Output
[[13 9 5 0 12]
[13 6 11 8 15]
[16 17 15 12 1]
[17 8 5 7 5]
[10 6 18 5 19]]
V, Q = eig(M)
print("Eigenvalues:n {}n".format(V))
print("Eigenvectors:n {}".format(Q))
## Output
Eigenvalues:
[50.79415691 +0.j 5.76076687+11.52079216j
5.76076687-11.52079216j -1.15784533 +3.28961651j
-1.15784533 -3.28961651j]
Eigenvectors:
[[ 0.34875973+0.j -0.36831427+0.21725348j -0.36831427-0.21725348j
-0.40737336-0.19752276j -0.40737336+0.19752276j]
[ 0.46629571+0.j -0.08027011-0.03330739j -0.08027011+0.03330739j
0.58904402+0.j 0.58904402-0.j ]
[ 0.50628483+0.j 0.62334823+0.j 0.62334823-0.j
-0.27738359-0.22063552j -0.27738359+0.22063552j]
[ 0.33975886+0.j 0.14035596+0.39427693j 0.14035596-0.39427693j
0.125282 +0.46663129j 0.125282 -0.46663129j]
[ 0.53774952+0.j -0.18591079-0.45968785j -0.18591079+0.45968785j
0.20856874+0.21329768j 0.20856874-0.21329768j]]
from numpy import diag
from numpy import dot
from numpy.linalg import inv
Q.dot(diag(V)).dot(inv(Q))
## Output
array([[1.30000000e+01-2.88657986e-15j, 9.00000000e+00-2.33146835e-15j,
5.00000000e+00+2.38697950e-15j, 1.17683641e-14+1.77635684e-15j,
1.20000000e+01-4.99600361e-16j],
[1.30000000e+01-4.32986980e-15j, 6.00000000e+00-3.99680289e-15j,
1.10000000e+01+3.38618023e-15j, 8.00000000e+00+1.72084569e-15j,
1.50000000e+01-2.77555756e-16j],
[1.60000000e+01-7.21644966e-15j, 1.70000000e+01-6.66133815e-15j,
1.50000000e+01+5.71764858e-15j, 1.20000000e+01+2.99760217e-15j,
1.00000000e+00-6.66133815e-16j],
[1.70000000e+01-5.27355937e-15j, 8.00000000e+00-3.10862447e-15j,
5.00000000e+00+4.27435864e-15j, 7.00000000e+00+2.22044605e-15j,
5.00000000e+00-1.22124533e-15j],
[1.00000000e+01-3.60822483e-15j, 6.00000000e+00-4.21884749e-15j,
1.80000000e+01+2.27595720e-15j, 5.00000000e+00+1.55431223e-15j,
1.90000000e+01+3.88578059e-16j]])
QR 分解
您可以通过应用 QR 分解将正方形或矩形矩阵（M）分解为正交矩阵（Q）和上三角矩阵（R）。 可以用以下公式表示：
M = QR
以下是 QR 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from numpy.linalg import qr
M = np.random.randint(low=0, high=20, size=20).reshape(4,5)
print(M)
## Output
[[14 6 0 19 3]
[ 9 6 17 8 8]
[ 4 13 17 4 4]
[ 0 0 2 7 11]]
Q, R = qr(M, 'complete')
print("Q:n {}n".format(Q))
print("R:n {}".format(R))
## Output
Q:
[[-0.81788873 0.28364908 -0.49345895 0.08425845]
[-0.52578561 -0.01509441 0.83834961 -0.14314877]
[-0.2336825 -0.95880935 -0.15918031 0.02718015]
[-0\. -0\. 0.16831464 0.98573332]]
R:
[[-17.11724277 -11.09991852 -12.91095786 -20.68090082 -7.59468109]
[ 0\. -10.85319349 -16.5563638 1.43333978 -3.10504542]
[ 0\. 0\. 11.88250752 -2.12744187 6.4411599 ]
[ 0\. 0\. 0\. 7.4645743 10.05937231]]
array([[1.40000000e+01, 6.00000000e+00, 1.77635684e-15, 1.90000000e+01,
3.00000000e+00],
[9.00000000e+00, 6.00000000e+00, 1.70000000e+01, 8.00000000e+00,
8.00000000e+00],
[4.00000000e+00, 1.30000000e+01, 1.70000000e+01, 4.00000000e+00,
4.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 2.00000000e+00, 7.00000000e+00,
1.10000000e+01]])
处理稀疏线性系统
您将不会总是使用密集矩阵，并且当您需要使用稀疏矩阵时，有些库将帮助您优化稀疏矩阵运算。 即使这些可能没有 Python API，您仍可能需要通过使用其他编程语言来使用它们，例如 C 和 C++：
- Hypre：包含预处理器和求解器，以利用并行实现来处理稀疏线性方程组。
- SuperLU：处理大型，稀疏，不对称的线性方程组。
- UMFPACK：解决稀疏线性方程组。
- CUSP：带有并行实现的稀疏线性代数和图形计算的开源库。 通过使用 CUSP，您可以访问 NVIDIA GPU 提供的计算资源。
- cuSPARSE：包含用于处理稀疏矩阵的线性代数子例程。 与 CUSP 一样，您可以访问 Nvidia GPU 提供的计算资源。
总结
在本章中，您探索了可以与 NumPy 配对的各种低级库及其配置。 我们特意运行了 EC2 条款，以便您熟悉基本的 Linux 命令行操作。 您还研究了各种计算密集型，数值，线性代数运算，这些运算将在下一章中用作基准测试不同的配置。
在下一章中，我们将创建一个基准 python 脚本，以在每种配置上运行。 您将能够查看不同线性代数运算和不同矩阵大小的性能指标
九、性能基准
在本章中，您将研究上一章介绍的不同配置的性能统计信息。 当然，当前设置无法为您提供最准确的环境，因为您无法控制 EC2 实例，但是它将使您了解自己环境中所需的设置。
我们将涵盖以下主题：
- 基准的重要性
- BLAS，LAPACK，OpenBLAS，ATLAS 和 Intel MKL 的性能
- 最终结果
为什么我们需要基准？
随着编程技巧的提高，您将开始实施更高效的程序。 您将搜索数十个代码存储库，以了解其他人如何解决类似的问题，并且您会发现那些令您赞叹不已的稀有宝石。
在编写更好的软件和实施系统的整个过程中，您将需要测量和跟踪改进速度的方法。 通常，您将以起点为基准，并查看所做的改进将如何累加性能指标。
设置基准后，您将对几种不同的实现进行基准测试，并将有机会根据您选择的性能指标进行比较。 您可以选择各种指标，并且需要事先决定。
这些基准的性能指标将保持相当简单，并且仅使用所花费的时间指标。 您将使用不同的配置多次执行相同的操作，并首先计算平均花费的时间。 计算平均值的公式为：
这是一个计算均值的好公式； 在我们的示例中，公式解释如下：
基线将基于此公式创建。 第一组计算如下：
的加法和乘法：
- 向量向量
- 向量矩阵
- 矩阵矩阵
通常，您将运行这些计算给定次数并计算平均值。
以下代码段向您展示了一个自定义函数，而不是 Python 中可用的通用计时器。 使用自定义函数的原因是，您以后可以将其与其他统计函数一起扩展，并通过适当的日志记录更好地查看详细信息。 函数将在计算开始之前输出有用的信息，并在迭代完成之后输出结果。
import inspect
import time
from datetime import datetime
def timer(*args, operation, n):
"""
Returns average time spent
for given operation and arguments.
Parameters
----------
*args: list (of numpy.ndarray, numpy.matrixlib.defmatrix.matrix or both)
one or more numpy vectors or matrices
operation: function
numpy or scipy operation to be applied to given arguments
n: int
number of iterations to apply given operation
Returns
-------
avg_time_spent: double
Average time spent to apply given operation
std_time_spent: double
Standard deviation of time spent to apply given operation
Examples
--------
>>> import numpy as np
>>> vec1 = np.array(np.random.rand(1000))
>>> vec2 = np.array(np.random.rand(1000))
>>> args = (vec1, vec2)
>>> timer(*args, operation=np.dot, n=1000000)
8.942582607269287e-07
"""
# Following list will hold the
# time spent value for each iteration
time_spent = []
# Configuration info
print("""
-------------------------------------------
### {} Operation ###
Arguments Info
--------------
args[0] Dimension: {},
args[0] Shape: {},
args[0] Length: {}
""".format(operation.__name__,
args[0].ndim,
args[0].shape,
len(args[0])))
# If *args length is greater than 1,
# print out the info for second argument
args_len = 0
for i, arg in enumerate(args):
args_len += 1
if args_len > 1:
print("""
args[1] Dimension: {},
args[1] Shape: {},
args[1] Length: {}
""".format(args[1].ndim,
args[1].shape,
len(args[1])))
print("""
Operation Info
--------------
Name: {},
Docstring: {}
Iterations Info
---------------
# of iterations: {}""".format(
operation.__name__,
operation.__doc__[:100] +
"... For more info type 'operation?'",
n))
print("""
-> Starting {} of iterations at: {}""".format(n, datetime.now()))
if args_len > 1:
for i in range(n):
start = time.time()
operation(args[0], args[1])
time_spent.append(time.time()-start)
else:
for i in range(n):
start = time.time()
operation(args[0])
time_spent.append(time.time()-start)
avg_time_spent = np.sum(time_spent) / n
print("""
-> Average time spent: {} seconds,
-------------------------------------------
""".format(avg_time_spent))
return avg_time_spent
当此函数中包含Docstring时，可以显示它以查看函数参数，返回的内容以及用法示例：
print(timer.__doc__)
这将生成以下输出：
Returns average time spent
for given operation and arguments.
Parameters
----------
*args: list (of numpy.ndarray, numpy.matrixlib.defmatrix.matrix or both)
one or more numpy vectors or matrices
operation: function
numpy or scipy operation to be applied to given arguments
n: int
number of iterations to apply given operation
Returns
-------
avg_time_spent: double
Average time spent to apply given operation
Examples
--------
>>> import numpy as np
>>> vec1 = np.array(np.random.rand(1000))
>>> vec2 = np.array(np.random.rand(1000))
>>> args = [vec1, vec2]
>>> timer(*args, operation=np.dot, n=1000000)
8.942582607269287e-07
让我们开始测量两个向量的点积所花费的平均时间。 以下代码块定义向量，并创建要输入到计时器函数中的参数：
import numpy as np
vec1 = np.array(np.random.rand(1000))
vec2 = np.array(np.random.rand(1000))
args = [vec1, vec2]
您现在可以按以下方式调用计时器函数：
timer(*args, operation=np.dot, n=1000000)
-------------------------------------------
### dot Operation ###
Arguments Info
--------------
args[0] Dimension: 1,
args[0] Shape: (1000,),
args[0] Length: 1000
args[1] Dimension: 1,
args[1] Shape: (1000,),
args[1] Length: 1000
Operation Info
--------------
Name: dot,
Docstring: dot(a, b, out=None)
Dot product of two arrays. Specifically,
- If both `a` and `b` are 1-D... For more info type 'operation?'
Iterations Info
---------------
# of iterations: 1000000
-> Starting 1000000 of iterations at: 2018-06-09 21:02:51.711211
-> Average time spent: 1.0054986476898194e-06 seconds,
-------------------------------------------
1.0054986476898194e-06
我们的向量乘积平均需要 1 微秒。 让我们看看如何通过添加其他指标来改进此计算。 您可以轻松添加的另一个指标是标准差，如下公式所示：
您熟悉上图中的公式术语。 标准差只是告诉您所报告指标的可变性，这是在我们的示例中花费的平均时间。
通过计算std_time_spent，打印其值并返回以下内容来扩展计时器函数：
avg_time_spent = np.sum(time_spent) / n
std_time_spent = np.std(time_spent)
print("""
-> Average time spent: {} seconds,
-> Std. deviation time spent: {} seconds
""".format(avg_time_spent, std_time_spent))
return avg_time_spent, std_time_spent
您还可以如下更新Docstring：
Returns
-------
avg_time_spent: double
Average time spent to apply given operation
std_time_spent: double
Standard deviation of time spent to apply given operation.
您可以重新定义时间函数，然后再次运行先前的计算，如下所示：
timer(*args, operation=np.dot, n=1000000)
您将获得以下输出（为简洁起见，仅显示最后一部分）以及其他信息：
-> Starting {} of iterations at: {}".format(n, datetime.now())
-> Average time spent: 1.0006928443908692e-06 seconds,
-> Std. deviation time spent: 1.2182541822530471e-06 seconds
(1.0006928443908692e-06, 1.2182541822530471e-06)
大！ 您还添加了哪些其他指标？ 如何添加置信区间？ 该部分将留给您锻炼，但是对您来说应该很容易！
让我们继续向量矩阵乘积：
mat1 = np.random.rand(1000,1000)
args = [vec1, mat1]
timer(*args, operation=np.dot, n=1000000)
这将为您提供以下输出：
Arguments Info
--------------
args[0] Dimension: 1,
args[0] Shape: (1000,),
args[0] Length: 1000
args[1] Dimension: 2,
args[1] Shape: (1000, 1000),
args[1] Length: 1000
Operation Info
--------------
Name: dot,
Docstring: dot(a, b, out=None)
Dot product of two arrays. Specifically,
- If both `a` and `b` are 1-D... For more info type 'operation?'
Iterations Info
---------------
# of iterations: 1000000
-> Starting 1000000 of iterations at: 2018-06-09 19:13:07.013949
-> Average time spent: 0.00020063393139839174 seconds,
-> Std. deviation time spent: 9.579314466482879e-05 seconds
(0.00020063393139839174, 9.579314466482879e-05)
最后，矩阵矩阵乘法如下：
mat1 = np.random.rand(100,100)
mat2 = np.random.rand(100,100)
args = [mat1, mat2]
timer(*args, operation=np.dot, n=1000000)
这将为您提供类似于先前输出的输出。
现在，我们或多或少有了一个想法，即如何挑战如何在计算机上执行这些任务。 基准函数列表已完成，在上一章中您看到了将点积添加到矩阵分解中的信息。
您将要做的是创建一个包含这些计算和统计信息的 Python 脚本文件。 然后，您将使用在 AWS 上设置的不同配置运行此文件。
让我们看一下linalg_benchmark.py，您可以在这个页面 中找到它。
以下代码块向您展示了linalg_benchmark.py脚本的重要部分，该脚本将用于测试您先前在 AWS 上设置的不同配置：
## Seed for reproducibility
np.random.seed(8053)
dim = 100
n = 10000
v1, v2 = np.array(rand(dim)), np.array(rand(dim))
m1, m2 = rand(dim, dim), rand(dim, dim)
## Vector - Vector Product
args = [v1, v2]
timer(*args, operation=np.dot, n=n)
## Vector - Matrix Product
args = [v1, m1]
timer(*args, operation=np.dot, n=n)
## Matrix - Matrix Product
args = [m1, m2]
timer(*args, operation=np.dot, n=n)
## Singular-value Decomposition
args = [m1]
timer(*args, operation=np.linalg.svd, n=n)
## LU Decomposition
args = [m1]
timer(*args, operation=lu, n=n)
## QR Decomposition
args = [m1]
timer(*args, operation=qr, n=n)
## Cholesky Decomposition
M = np.array([[1, 3, 4],
[2, 13, 15],
[5, 31, 33]])
args = [M]
timer(*args, operation=cholesky, n=n)
## Eigenvalue Decomposition
args = [m1]
timer(*args, operation=eig, n=n)
print("""
NumPy Configuration:
--------------------
""")
np.__config__.show()
将有两个单独的运行：
- 1 st 与
dim = 100一起运行 - 2 和与
dim = 500一起运行
让我们看一下结果。
准备性能基准
对于每个实例和配置，导航到您的Home目录并创建一个名为py_scripts的文件夹：
使用以下命令创建名为linalg_benchmark.py的文件并粘贴内容：
粘贴内容后，键入:，然后键入wq!和Enter保存并退出：
现在，您可以使用以下命令运行该文件：
对于 Anaconda 分发，您将使用以下命令运行脚本：
BLAS 和 LAPACK 的性能
在这里，您将使用 BLAS 和 LAPACK 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后如上一节中所示运行脚本。
以下是dim = 100的运行结果：
以下是dim = 500的运行结果：
OpenBLAS 的性能
在这里，您将使用 OpenBLAS 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后运行上一节中显示的脚本。
以下是dim = 100的运行结果：
以下是dim = 500的运行结果：
ATLAS 的性能
在这里，您将使用 ATLAS 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，如上一节中所示运行脚本。
以下是dim = 100的运行结果：
以下是dim = 500的运行结果：
英特尔 MKL 的性能
在这里，您将使用英特尔 MKL 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后运行上一节中显示的脚本。
以下是dim = 100的运行结果：
以下是dim = 500的运行结果：
结果
当然，t2.micro实例相当薄弱，您应该更多地了解 Amazon 如何为 EC2 实例提供这种计算能力。 您可以在这个页面 上阅读有关它们的更多信息。
如果您使用功能更强大的计算机并具有更多的内核，则不同配置之间的性能差异将更加明显。
说到结果，毫不奇怪，默认安装的 BLAS 和 LAPACK 为我们提供了基准性能，而经过优化的版本（如 OpenBLAS，ATLAS 和 Intel MKL）提供了更好的性能。
正如您已经指出的那样，您没有在 Python 脚本中更改任何代码行，而仅通过将 NumPy 库与不同的加速器链接起来，便获得了巨大的性能提升。
如果您将更深入地研究这些低级库以了解提供了哪些特定的例程和函数，则将更好地了解程序的哪些部分将从这些实现中受益。
当然，起初您可能还不了解许多其他细节。 可能是您使用的函数未使用低级库或未并行化操作的情况。 在某些情况下，多线程会或不会有所帮助。 知识和经验最终取决于您的实验，并且您将从自己的经验中学到东西，因此您将更加精通各种应用。
许多研究人员发表了实验的设计和结果。 Google 的快速搜索将为您提供大量资源，以阅读和了解这些库在不同硬件和软件配置下的性能。
总结
在本章中，您探讨了执行计算密集型线性代数运算时不同配置的性能。
基准测试是一项严肃的工作，您至少现在已经具备运行基准测试的基本技能。 您在本章中学习的材料远未完成，但是它为您提供了从哪里开始的想法，并且您肯定可以在许多方面进行改进。
您可以看到的一件事是，逐渐增加向量和矩阵的大小时性能指标的行为。 理想情况下，您将需要功能更强大的硬件，但是t2.micro实例在大多数情况下是免费的，或者提供的价格非常便宜。
由于您将需要处理更多的计算密集型工作负载，因此重要的是要了解您的选择以及哪种选择将为您带来最佳性能。 您可以运行这些简单的实验，至少对性能有所了解，这将为您带来很多帮助，并节省时间和金钱。
如果您走了这么远，恭喜！ 我们认为，遍历所有章节并学习相关材料可以提高您在 Python 科学堆栈方面的技能。
我们希望您喜欢阅读这本书，并感谢您的宝贵时间。
浙公网安备 33010602011771号