Python科学计算
NumPy库是用于科学计算的一个开源Python扩充程序库，是其他数据分析包的基础包，它为Python提供了高性能数组与矩阵运算处理能力。
- ndarray数组：学会数组的创建和属性。
- 数组选择：学会数组的切片和索引方法。
- 数组运算：学会数组的各类运算方法和使用。
- 数组存取：完成数组的存储和读取方法。
- 图像变换：了解图像的处理方法。
NumPy简介
NumPy是一个开源的Python科学计算的库，主要用在数据分析和科学计算，基本上可以把NumPy看做是Python数据计算的基础，因为很多非常优秀的数据分析和机器学习框架底层使用的都是NumPy。比如：Pandas, SciPy, Matplotlib, scikit-learn, scikit-image 等。
NumPy是一个Python包。它代表Numeric Python。NumPy库主要包含多维数组和矩阵数据结构。它为ndarray（一个n维数组对象）提供了对其进行有效操作的方法。NumPy可以用于对数组执行各种数学运算。并且提供了可在这些数组和矩阵上运行的庞大的高级数学函数库。
Numeric是NumPy的前身，是由Jim Hugunin开发的。他还开发了另一个包Numarray，拥有一些额外的功能。2005年，Travis Oliphant通过将Numarray的功能集成到Numeric包中来创建NumPy包。这个开源项目有很多贡献者。
使用NumPy，开发人员可以执行以下操作：
- 数组的算术和逻辑运算。
- 傅里叶变换和用于图形操作的例程。
- 与线性代数有关的操作，NumPy拥有线性代数和随机数生成的内置函数。 NumPy是SciPy家族的成员之一。SciPy家族是一个专门应用于数学、科学和工程领域的开源Python生态圈，SciPy家族的核心成员为Matplotlib、SciPy和NumPy，可以概括为MSN这三个字母。 NumPy通常与SciPy（Scientific Python）和Matplotlib（绘图库）一起使用。这种组合广泛用于替代Matlab，是一个流行的技术计算平台。Python作为Matlab的替代方案，现在被视为一种更加现代和完整的编程语言。
NumPy概述
NumPy是Python科学计算的基础软件包，提供了多维数组类（numpy.ndarray）及其派生类（掩码数组、矩阵等），以及用于快速操作数组的函数和API，涵盖基本线性代数、基本统计运算、随机模拟等诸多数学和工程领域。 和Python的列表相比，NumPy数组在运行速度上拥有明显的优势。NumPy底层使用C语言编写，内置并行运算功能，并且内部解除了GIL（Global Interpreter Lock，全局解释器锁），这意味着NumPy在运行速度和并行计算方面有着先天优势。
- NumPy对数组的操作速度不受Python解释器的限制。
- 当系统有多个CPU时，NumPy可以自动进行并行计算。 NumPy的数据组织结构，尤其是多维数组类（numpy.ndarray），几乎已经成为所有数据处理与机器学习的标准数据结构。越来越多基于Python的机器学习和数据处理软件包开始使用NumPy数组，虽然这些软件包通常都支持Python列表作为参数，但它们在处理之前还是会将输入的Python列表转换为NumPy数组，而且输出也通常为NumPy数组。 在Python的阵营里，NumPy的重要性和普遍性日趋增强。换句话说，为了高效地使用机器学习和数据处理等基于Python的工具包，只知道如何使用Python列表是不够的，还需要知道如何使用NumPy数组。
NumPy模块
NumPy的核心是多维数组类numpy.ndarray，矩阵类numpy.matrix是多维数组类的派生类。以多维数组类为数据组织结构，NumPy提供了众多的数学、科学和工程函数。此外，NumPy还提供了以下多个子模块。
- numpy.random：随机抽样子模块。
- numpy.linalg：线性代数子模块。
- numpy.fft：傅里叶变换子模块。
- numpy.ctypeslib：C-Types外部函数接口子模块。
- numpy.emath：具有自动域的数学函数子模块。
- numpy.testing：测试支持子模块。
- numpy.matlib：矩阵库子模块。
- numpy.dual：SciPy加速支持子模块。
- numpy.distutils：打包子模块。
NumPy环境安装配置
标准的Python发行版不会与NumPy模块捆绑在一起。NumPy的安装非常简单，可以使用pip命令直接安装。安装时还可以使用==选择版本号，使用-i参数选择下载速度更快的镜像源。
pip install numpy
导入NumPy时写成import numpy as np是程序员约定俗成的规则。
import numpy as np
列表VS数组
Python的列表操作非常灵活，特别是引入负数索引后，更是为列表操作注入了“神奇的力量”。每一位程序员第一次使用Python列表时，都会被它的便捷所打动。
- 首先，Python列表的元素类型不受限制。同一个列表内，列表元素可以是不同的数据类型，甚至可以是函数。
- 其次，Python列表的元素可以动态增减。例如，append( )方法可以向列表末尾追加元素，insert( )方法可以在指定位置插入元素，pop( )方法可以删除指定索引的元素，remove( )方法可以删除指定的元素。
- 最后，Python列表的索引非常简单、灵活、高效。例如，
[-1]可以取得最末位的元素，[1:-1]可以返回“掐头去尾”后的列表，[::-1]可以反转列表，[::2]可以隔一个取一个元素。
NumPy数组的操作便捷性，比Python列表更是有过之而无不及，并且还提供了大量的函数，但同时也增加了限制：一是数组元素必须具有相同的数据类型；二是数组一旦创建，其元素数量就不能再改变了。
说到这里，你也许会问：既然NumPy数组有这些限制条件，那我们为什么还要使用它呢？答案是，NumPy数组具有极高的、接近C语言的运行效率，同时又继承了Python列表操作便捷、灵活的特点，可以说NumPy数组是专为处理科学数据而生的。
数组的数据类型
NumPy支持的数据类型主要有整型（integer）、浮点型（float）、布尔型（bool）和复数型（complex），每一种数据类型根据占用内存的字节数又分为多个不同的子类型。 事实上，NumPy也支持字符串类型和自定义类型的数据，但绝大多数函数和方法不适用于非数值型数组。
创建数组
我们创建array，使用的是np.arange方法。我们还可以通过List来创建Array，List可以是一维列表，也可以是多维列表： 创建数组时，如果不指定数据类型，NumPy会根据输入的参数自动选择合适的数据类型。通常在指定数据类型的时候，可以省略类型后面的数字。如果省略数字，整型和无符号整型默认是32位，浮点型默认是64位，复数型默认是128位。
下面的代码生成数组并显示其数据类型。
import numpy as np
a = np.array([0, 1, 2])
print(a.dtype)
# int32
a = np.array([0, 1, 2.0])
print(a.dtype)
# float64
a = np.array([0, 1, 2], dtype=np.uint8)
print(a.dtype)
# uint8
使用正确的数据类型很重要。在图像处理、三维显示等领域，如果数据类型不正确就可能得不到正确的结果，并且很多情况下没有任何错误提示。例如，使用NumPy数组处理图像时，数据类型一般都必须指定为np.uint8类型。
数组的属性
我们用数组的dtype属性查看数组的数据类型。除了dtype属性，数组对象还有一些其他的属性，如shape属性表示数组的结构或叫作形状；size属性表示数组元素个数；itemsize属性表示数组元素字节数；flags属性表示数组的内存信息；real属性表示数组实部；imag属性表示数组虚部；data属性表示实际存储区域内存的起始地址，相当于指针。 数组的属性看起来有点多，但只需要记住dtype和shape这两个属性就足够了。这两个属性非常重要，重要到可以忽略其他属性。
数组的方法
数组对象的内置方法更多，重要的两个方法：ndarray.astype( )和ndarray.reshape( )。前者可以修改元素类型，后者可以重新定义数组的结构。 这两个方法的重要性和其对应的属性一样。记住ndarray.dtype和ndarray.shape这两个属性及其对应的修改数据类型和数组结构的两个方法，有助于在调试代码时快速定位问题。
下面的代演示了如何修改数组结构和数据类型。
import numpy as np
a = np.arange(6)
print(a)
# array([0, 1, 2, 3, 4, 5])
print(a.shape)
# (6,)
print(a.dtype)
# dtype('int32')
a = a.reshape((2,3)) # 改变数组结构
print(a)
#array([[0, 1, 2],[3, 4, 5]])
a = a.astype(np.float64) # 改变数据类型
print(a.dtype)
# dtype('float64')
接下来我们介绍几个常用的名词：
- vector — 表示的是一维数组
- matrix — 表示的是二维数组
- tensor — 表示的是三维或者更高维度的数组 在NumPy中维度也被称之为 axes 。
维、秩、轴
维，就是维度。通常说数组是几维的，就是指维度数，如三维数组的维度数就是3。维度数还有一个专用名字，即秩，也就是数组属性ndim。秩这个名字感觉有些多余，不如维度数更容易理解。 但是轴的概念大家一定要建立起来，并且要理解，因为轴的概念很重要。简单来说，我们可以把数组的轴和笛卡儿坐标系的轴来对应一下。 一维数组，类比于一维空间，只有一个轴，那就是0轴。二维数组，类比于二维空间，有两个轴，习惯表示成行和列，行的方向是0轴，列的方向是1轴。三维数组，类比于三维空间，有三个轴，习惯表示成层、行和列，层的方向是0轴，行的方向是1轴，列的方向是2轴。
下面用一个三维数组求总和与分层求和的例子来演示一下轴概念的重要性。 我们知道，列表求和需要使用Python内置的求和函数sum( )，且只能对数值型的一维列表求和。而数组则是自带求和方法，且支持按指定轴的方向求和，其代码如下。
import numpy as np
a = np.arange(18).reshape((3,2,3)) # 3层2行3列的结构
print(a)
#array([[[ 0, 1, 2],
# [ 3, 4, 5]],
#
# [[ 6, 7, 8],
# [ 9, 10, 11]],
#
# [[12, 13, 14],
# [15, 16, 17]]])
a.sum() # 全部数组元素之和
# 153
a.sum(axis=0) # 0轴方向求和：3层合并成1层，返回二维数组
# array([[18, 21, 24],
# [27, 30, 33]])
a.sum(axis=1) # 1轴方向求和：2行合并成1行，返回二维数组
# array([[ 3, 5, 7],
# [15, 17, 19],
# [27, 29, 31]])
a.sum(axis=2) # 2轴方向求和：3列合并成1列，返回二维数组
# array([[ 3, 12],
# [21, 30],
# [39, 48]])
a.sum(axis=1).sum(axis=1) # 分层求和方法1
# array([15, 51, 87])
a.sum(axis=2).sum(axis=1) # 分层求和方法2
# array([15, 51, 87])
广播和矢量化
NumPy数组具有极高的、接近C语言的运行效率，处理速度远比Python列表快得多。为什么数组比列表的处理速度快呢？原来，ndarray拥有区别于列表的两大“独门绝技”：广播（broadcast）和矢量化（vectorization）。广播可以理解为隐式地对每个元素实施操作；矢量化可以理解为代码中没有显式的循环、索引等。广播和矢量化对于初学者而言有点抽象，下面我们用两个例子来说明一下。
数值型数组的各个元素加1
使用Python列表实现列表的各个元素加1，似乎除了循环就没有更好的办法了。当然，如果你想到了用map( )函数来实现，说明你有扎实的基本功，但这个方法只是避免显式地使用循环，实际处理速度不会比循环更快。
a = [0, 1, 2, 3, 4]
for i in range(len(a)):
a[i] += 1
print(a)
# [1, 2, 3, 4, 5]
如果换成NumPy数组，利用其广播特性，无须循环就可以实现对数组每一个元素加1的操作，其代码如下。数组的广播特性，不仅省略了循环结构，更重要的是可以大幅度加快数据的处理速度。
import numpy as np
a = np.array([0, 1, 2, 3, 4])
a += 1
a
# array([1, 2, 3, 4, 5])
两个等长的数值型数组的对应元素相加
如果两个等长的Python列表对应元素相加，需要同时遍历两个列表，单是想把代码写“漂亮”就需要花一点心思。下面的代码使用zip( )函数，辅以列表推导式，实现两个等长Python列表的对应元素相加。虽然代码形式上接近“完美”，但是运行效率仍然很低。
a = [0, 1, 2, 3, 4]
b = [5, 6, 7, 8, 9]
[i+j for i, j in zip(a, b)]
# [5, 7, 9, 11, 13]
如果换成NumPy数组，利用其矢量化特性来实现两个数组对应元素相加，就像是进行两个整型变量相加，代码可读性强、处理速度快，其代码如下。
a = np.array([0, 1, 2, 3, 4])
b = np.array([5, 6, 7, 8, 9])
a + b
# array([ 5, 7, 9, 11, 13])
上面的两个例子分别用列表和数组的方式给出了答案。显然，用NumPy数组实现起来要比用Python列表更直观、更简洁。这正是得益于NumPy的两大特性：广播和矢量化。 广播和矢量化是NumPy最“精髓”的特性，是NumPy的灵魂。广播和矢量化体现在代码上则有以下几个特点。
- 矢量化代码更简洁，更易于阅读。
- 代码行越少意味着出错的概率越小。
- 代码更接近标准的数学符号。 矢量化代码更pythonic（意思是更有Python的“味道”）。
创建数组
NumPy为创建数组提供了非常丰富的手段，可以无中生有，可以移花接木，还可以举一反三。配合数据类型设置、结构设置，可以创建出任何形式的数组。
创建数组分成了创建简单数组和创建复杂数组两大类。其实简单数组和复杂数组并没有严格的分界线，大致上，“无中生有”创建出来的数组称为简单数组，如蛮力构造法、特殊数值法、随机数值法和定长分割法等；通过“移花接木”“举一反三”创建出来的数组称为复杂数组，如重复构造法、网格构造法等。
蛮力构造法
蛮力构造法使用np.array( )函数来创建数组，其原型如下。
np.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
该函数的参数看上去挺多，但固定参数（必要参数）只有一个object，也就是和我们要创建的数组相似的（array_like）数据结构，通常是Python列表或元组。
下面的例子使用一个列表创建数组，如果改成元组也完全符合规则。
import numpy as np
a = np.array([[1,2,3],[4,5,6]]) # 创建2行3列数组
a
# array([[1, 2, 3],
# [4, 5, 6]])
a.dtype
# dtype('int32')
在np.array( )函数的默认参数中，dtype参数用于指定数据类型。创建数组时，如果不指定数据类型，np.array( )函数会根据object参数自动选择合适的数据类型。
当然，也可以像下面代码中演示的这样，在创建数组时，指定元素的数据类型。
import numpy as np
a = np.array([[1,2,3],[4,5,6]], dtype=np.uint8) # 创建8位无符号整型数组
a
# array([[1, 2, 3],
# [4, 5, 6]], dtype=uint8)
蛮力构造法就是将想要创建数组的数据结构直接用Python列表或元组写出来，再用np.array( )函数转为数组。这个方法虽然看起来简单，但很容易出错，不适合创建体量较大的数组。
特殊数值法
这里的特殊数值指的是0、1、空值。特殊数值法适合构造全0、全1、空数组，或由0、1组成的类似单位矩阵（主对角线为1，其余为0）的数组。特殊数值法使用的4个函数原型如下。
np.zeros(shape, dtype=float, order='C')
np.ones(shape, dtype=float, order='C')
np.empty(shape, dtype=float, order='C')
np.eye(N, M=None, k=0, dtype=float, order='C')
固定参数shape表示生成的数组结构，默认参数dtype用于指定数据类型（默认浮点型）。虽然order参数几乎用不到，但作为常识，我们有必要了解一下。order参数指定的是数组在内存中的存储顺序，“C”表示C语言使用的行优先方式，“F”表示Fortran语言使用的列优先方式。
使用上面4个函数配合shape和dtype参数，可以很方便地创建出一些简单数组，其代码如下。
import numpy as np
np.zeros(6)
# array([0., 0., 0., 0., 0., 0.])
np.zeros((2,3))
# array([[0., 0., 0.],
# [0., 0., 0.]])
np.ones((2,3),dtype=np.int)
# array([[1, 1, 1],
# [1, 1, 1]])
np.empty((2,3))
# array([[1., 1., 1.],
# [1., 1., 1.]])
np.eye(3, dtype=np.uint8)
# array([[1, 0, 0],
# [0, 1, 0],
# [0, 0, 1]], dtype=uint8)
如果需要一个3行4列、初始值都是255的无符号整型数组，应该怎么做呢？全1数组乘以255，或全0数组加255，都是很好的解决方案。另外，使用填充函数fill( )也可以解决这个问题。fill( )函数不只可以填充空数组，任何数组都可以使用它来填充固定的值，其代码如下。
import numpy as np
a = np.empty((3,4), dtype=np.uint8)
a.fill(255)
a
# array([[255, 255, 255, 255],
# [255, 255, 255, 255],
# [255, 255, 255, 255]], dtype=uint8)
随机数值法
和Python的标准模块random类似，NumPy有一个random子模块，其功能更加强大。用随机数值法创建数组主要就是使用random子模块。random子模块的方法很多，这里只介绍3个最常用的函数。这3个最常用的函数原型如下。
np.random.random(size=None)
np.random.randint(low, high=None, size=None)
np.random.normal(loc=0.0, scale=1.0, size=None)
random( )函数用于生成 [0,1) 区间内的随机浮点型数组，randint()函数用于生成[low, high)区间内的随机整型数组。参数size是一个元组，用于指定生成数组的结构，其代码如下。请注意，这里描述的[0,1)区间和[low, high)区间都是左闭右开的。
import numpy as np
np.random.random(3)
# array([0.4129063 , 0.94242888, 0.10129428])
np.random.random((2,3))
# array([[0.80530845, 0.96161533, 0.89166972],
# [0.22920038, 0.84989557, 0.46865645]])
np.random.randint(5)
# 4
np.random.randint(1, 5, size=(2,3))
# array([[1, 4, 1],
# [2, 4, 3]])
normal()函数用于生成以loc为均值、以scale为标准差的正态分布数组。
下面用正态分布函数模拟生成1000位成年男性的身高数据（假定成年男性平均身高为170厘米，标准差为4厘米），并画出柱状图。
import numpy as np
import matplotlib.pyplot as plt # 导入绘图模块
tall = np.random.normal(170, 4, 1000) # 生成正态分布数据
bins = np.arange(156, 190, 2) # 从156厘米到190厘米，每2厘米一个分段
plt.hist(tall, bins) # 绘制柱状图
plt.show() # 显示图形
定长分割法
定长分割法最常用的函数是arange( )，它看起来和Python的range( )函数很像，只是前面多了一个字母a。另一个常用的定长分割函数是linspace( )，类似于arange( )函数，但功能更加强大，两个函数的原型如下。
np.arange(start, stop, step, dtype=None)
np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)
arange()函数和Python的range()函数用法相同，并且还可以接收浮点型参数，其代码如下。
import numpy as np
np.arange(5)
# array([0, 1, 2, 3, 4])
np.arange(5, 11)
# array([ 5, 6, 7, 8, 9, 10])
np.arange(5,11,2)
# array([5, 7, 9])
np.arange(5.5, 11, 1.5)
# array([ 5.5, 7. , 8.5, 10. ])
np.arange(3,15).reshape(3,4)
# array([[ 3, 4, 5, 6],
# [ 7, 8, 9, 10],
# [11, 12, 13, 14]])
linspace( )函数需要3个参数：一个起点、一个终点、一个返回元素的个数。linspace( )函数返回的元素包括起点和终点，我们可以通过endpoint参数选择是否包含终点，其代码如下。
import numpy as np
np.linspace(0, 5, 5) # 返回0到5之间的5个等距数值，包括0和5
# array([0. , 1.25, 2.5 , 3.75, 5. ])
np.linspace(0, 5, 5, endpoint=False) # 返回5个等距数值，包括0但不包括5
# array([0., 1., 2., 3., 4.])
重复构造法
重复构造法，顾名思义就是根据特定的规则对已有数组不断重复，从而生成新的数组。重复构造法主要使用repeat( )和tile( )这两个函数。 一般而言，repeat()函数用来重复数组元素。但如果被重复的数组是一个多维数组，且repeat( )函数指定了axis参数，情况就会变得有些复杂，其代码如下。
import numpy as np
a = np.arange(5)
a
# array([0, 1, 2, 3, 4])
np.repeat(a, 3) # 重复一维数组元素3次
# array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4])
a = np.arange(6).reshape((2,3))
a
# array([[0, 1, 2],
# [3, 4, 5]])
np.repeat(a, 3) # 重复二维数组元素3次，不指定轴
# array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5])
np.repeat(a, 3, axis=0) # 重复二维数组元素3次，指定0轴
# array([[0, 1, 2],
# [0, 1, 2],
# [0, 1, 2],
# [3, 4, 5],
# [3, 4, 5],
# [3, 4, 5]])
np.repeat(a, 3, axis=1) # 重复二维数组元素3次，指定1轴
# array([[0, 0, 0, 1, 1, 1, 2, 2, 2],
# [3, 3, 3, 4, 4, 4, 5, 5, 5]])
tile的原意是铺地砖或贴墙砖，总之是把一块一块的地砖或墙砖，一排排一列列地排列整齐。tile( )函数也是如此，它将整个数组而非数组元素水平和垂直重复指定的次数。因为没有axis参数，所以tile( )函数相对容易理解，其代码如下。
import numpy as np
a = np.arange(5)
a
# array([0, 1, 2, 3, 4])
np.tile(a, 3) # 重复一维数组3次
# array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4])
np.tile(a, (3,2)) # 重复一维数组3行2列
# array([[0, 1, 2, 3, 4, 0, 1, 2, 3, 4],
# [0, 1, 2, 3, 4, 0, 1, 2, 3, 4],
# [0, 1, 2, 3, 4, 0, 1, 2, 3, 4]])
a = np.arange(6).reshape((2,3))
a
# array([[0, 1, 2],
# [3, 4, 5]])
np.tile(a, 3) # 重复二维数组3次
# array([[0, 1, 2, 0, 1, 2, 0, 1, 2],
# [3, 4, 5, 3, 4, 5, 3, 4, 5]])
np.tile(a, (2,3)) # 重复二维数组2行3列
# array([[0, 1, 2, 0, 1, 2, 0, 1, 2],
# [3, 4, 5, 3, 4, 5, 3, 4, 5],
# [0, 1, 2, 0, 1, 2, 0, 1, 2],
# [3, 4, 5, 3, 4, 5, 3, 4, 5]])
网格构造法
众所周知，研究地球表面需要经纬度坐标，经度从西经180°（-180°）到东经180°（180°），纬度从北纬90°（90°）到南纬90°（-90°），把经纬度线画出来，就形成了一个经纬度网格。经纬度网格是科学数据中常用的概念。 通常，经度用longitude表示，简写为lon，纬度用latitude表示，简写为lat。那么，如何用数组表示经纬度网格呢？
用数组表示经纬度网格一般有两种方式。第一种方式，用两个一维数组表示。下面的代码使用定长分割函数linspace( )，将经度从-180°到180°分为间隔为10°的37个点，将纬度从90°到-90°分为间隔为10°的19个点，得到两个一维数组。
lon = np.linspace(-180,180,37) # 精度为10°，共计37个经度点
lat = np.linspace(90,-90,19) # 精度为10°，共计19个纬度点
经纬度网格的第二种表示方式是用两个二维数组分别表示经度网格和纬度网格。经度网格中每一列的元素都是相同的（同一个经度），纬度网格中每一行的元素都是相同的（同一个纬度）。
生成二维经纬度网格的常用函数是np.meshgrid( )，该函数以一维经度数组lon和一维纬度数组lat为参数，返回二维的经度数组和纬度数组，其代码如下。
lons,lats = np.meshgrid(lon,lat)
lons.shape
# (19, 37)
lats.shape
# (19, 37)
lons[:,0]
# array([-180., -180., -180., -180., -180., -180., -180., -180., -180.,
# -180., -180., -180., -180., -180., -180., -180., -180., -180.,
# -180.])
lats[0]
# array([90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90.,
# 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90.,
# 90., 90., 90., 90., 90., 90., 90., 90., 90., 90., 90.])
从上面的代码中可以看出，二维经度数组lons的第0列所有元素都是-180°，二维纬度数组lats的第0行所有元素都是90°。
构造经纬度网格，除了使用np.meshgrid( )函数外，还有一个更强大的方法。这个方法可以直接生成纬度网格和经度网格而无须借助于一维数组（请注意，纬度在前，经度在后）。
lats, lons = np.mgrid[90:-91:-5, -180:181:5] # 用实数指定网格精度为5°
lons.shape, lats.shape
# ((37, 73), (37, 73))
lats, lons = np.mgrid[90:-90:37j, -180:180:73j] # 也可以用虚数指定分割点数
lons.shape, lats.shape
# ((37, 73), (37, 73))
上面的例子中用到了虚数。相应构造复数的方法如下。
r, i = 2, 5
complex(r, i)
# (2+5j)
自定义数据类型
在讲解数据类型时已经说过，NumPy也支持字符串类型和自定义类型，但绝大多数函数和方法不适用于非数值型数组，因此，自定义数据类型将是最后的选择。
我们先来思考一个问题：同一个列表中，元素类型既有字符串，又有整型和浮点型，将该列表转成数组，会报错吗？如果不报错，数组的数据类型是什么呢？下面用一个例子演示一下，其代码如下。
np.array(['Anne', 1.70, 55])
# array([Anne, '1.70', '55'], dtype='<U4')
结果显示，数组会将所有元素的数据类型都转为’<U4’类型。这里的U表示Unicode字符串；<表示字节顺序，意为小端在前（低位字节存储在最小地址中）；4表示数组元素占用4字节，数组元素占用的字节数由所有元素中最长的那个元素决定。
接下来我们继续思考：怎样在数组中保留用以生成数组的列表中的元素类型呢？这就需要用到自定义数据类型了。自定义数据类型类似于C语言的结构体，其代码如下。
>>> mytype = np.dtype([('name','S32'), ('tall',np.float), ('bw',np.int)])
>>> np.array([('Anne', 1.70, 55)], dtype=mytype)
array([(b'Anne', 1.7, 55)],
dtype=[('name', 'S32'), ('tall', '<f8'), ('bw', '<i4')])
操作数组
数组的索引和切片、合并与拆分、复制和排序、查找和筛选，以及改变数组结构、数组I/O等，这些是数组操作的基本技术。其中最抽象的是查找和筛选，但这也是数组操作中最重要、最精髓的一部分。在数组操作中用好查找和筛选才能避免使用循环，这是数组操作的最高境界。
索引和切片
索引是定位一维或多维数组中的单个或多个元素的行为模式。切片是返回一维或多维数组中的单个或多个相邻元素的视图，目的是引用或赋值。NumPy数组对象的内容可以通过索引或切片来访问和修改。对于一维数组的索引和切片，NumPy数组使用起来和Python的列表一样灵活。
>>> a = np.arange(9)
>>> a[-1] # 最后一个元素
8
>>> a[2:5] # 返回第2到第5个元素
array([2, 3, 4])
>>> a[:7:3] # 返回第0到第7个元素，步长为3
array([0, 3, 6])
>>> a[::-1] # 返回逆序的数组
array([8, 7, 6, 5, 4, 3, 2, 1, 0])
对于多维数组操作，NumPy数组比Python的列表更加灵活、强大。假设有一栋楼，共2层，每层的房间都是3行4列，那我们可以用一个三维数组来保存每个房间的居住人数（也可以是房间面积等其他数值信息）。
>>> a = np.arange(24).reshape(2,3,4) # 2层3行4列
>>> a
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]]])
>>> a[1][2][3] # 虽然可以这样索引
23
>>> a[1,2,3] # 但这样才是规范的用法
23
>>> a[:,0,0] # 所有楼层的第0行第0列
array([ 0, 12])
>>> a[0,:,:] # 1层的所有房间，等价于a[0]或a[0,...]
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a[:,:,1:3] # 所有楼层所有排的第1到第3列
array([[[ 1, 2],
[ 5, 6],
[ 9, 10]],
[[13, 14],
[17, 18],
[21, 22]]])
>>> a[1,:,-1] # 2层每一行的最后一个房间
array([15, 19, 23])
从上面的代码中可以看出，对多维数组索引或切片得到的结果的维度不是确定的。另外还有一点需要特别提醒：切片返回的数组不是原始数据的副本，而是指向与原始数组相同的内存区域。数组切片不会复制内部数组数据，只是产生了原始数据的一个新视图。
>>> a = np.arange(12).reshape(3,4)
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> b = a[1:,2:] # 数组b是数组a的切片
>>> b
array([[ 6, 7],
[10, 11]])
>>> b[:,:] = 99 # 改变数组b的值，也会同时影响数组a
>>> b
array([[99, 99],
[99, 99]])
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 99, 99],
[ 8, 9, 99, 99]])
上面的代码中，数组b是数组a的切片，当改变数组b的元素时，数组a也同时发生了改变。这就证明了切片返回的数组不是一个独立数组，而是指向与原始数组相同的内存区域。
改变数组结构
NumPy之所以拥有极高的运算速度，除了并行、广播和矢量化等技术因素外，其数组存储顺序和数组视图相互独立也是一个很重要的原因。正因为如此，改变数组结构自然是非常便捷的操作。改变数组结构的操作通常不会改变所操作的数组本身的存储顺序，只是生成了一个新的视图。np. resize( )函数是个例外，它不返回新的视图，而是真正改变了数组的存储顺序。
ndarray自带多个改变数组结构的方法，在大部分情况下学会ndarray.reshape( )函数即可，我们在前面已经多次用到该函数。在某些情况下，翻滚轴函数numpy.rollaxis( )才是最佳的选择，需要多花一些时间去了解它。以下是改变数组结构的几个常用函数。
- ndarray.reshape( )：按照指定的结构（形状）返回数组的新视图，不改变原数组。
- ndarray.ravel( )：返回多维数组一维化的视图，不改变原数组。
- ndarray.transpose( )：返回行变列的视图，不改变原数组。
- ndarray.resize( )：按照指定的结构（形状）改变原数组，无返回值。
- numpy.rollaxis( )：翻滚轴，返回新的视图，不改变原数组。
如果你足够细心，读到这里也许会产生些许疑惑：为什么前面的几个函数都是ndarray的方法，而翻滚轴函数numpy.rollaxis( )却是numpy的呢？这是命名空间的问题，并且NumPy在命名空间问题上的确有些含糊不清。不过不用担心，在下一节会专门解释这个问题。下面继续演示这几个改变数组结构的函数的用法。
>>> a = np.arange(12) >>> b = a.reshape((3,4)) # reshape()函数返回数组a的一个新视图，但不会改变数组a >>>> a.shape (12,) >>> b.shape (3, 4) >>> b is a False >>> b.base is a True a.resize([4,3]) # resize()函数没有返回值，但真正改变了数组a的结构 >>> a.shape (4, 3) >>> a.ravel() # 返回多维数组一维化的视图，但不会改变原数组 array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]) >>> a.transpose() # 返回行变列的视图，但不会改变原数组 array([[ 0, 3, 6, 9], [ 1, 4, 7, 10], [ 2, 5, 8, 11]]) >>> a.T # 返回行变列的视图，等价于transpose()函数 array([[ 0, 3, 6, 9], [ 1, 4, 7, 10], [ 2, 5, 8, 11]]) >>> np.rollaxis(a, 1, 0) # 翻滚轴，1轴变0轴 array([[ 0, 3, 6, 9], [ 1, 4, 7, 10], [ 2, 5, 8, 11]])
翻滚轴函数有一个很容易理解的应用，就是用它来实现图像的通道分离。下面的代码生成了一个宽为800像素、高为600像素的彩色随机噪声图，使用翻滚轴函数可以将其分离成RGB三个颜色通道。最后两行代码导入pillow模块，从而可以直观地看到这幅噪声图。
>>> img = np.random.randint(0, 256, (600, 800, 3), dtype=np.uint8)
>>> img.shape
(600, 800, 3)
>>> r, g, b = np.rollaxis(img, 2, 0) # 将图像数据分离成RGB三个颜色通道
>>> r.shape, g.shape, b.shape
((600, 800), (600, 800), (600, 800))
>>> from PIL import Image # 导入pillow模块的Image
>>> Image.fromarray(img).show() # 显示随机生成的噪声图
合并
NumPy数组一旦创建就不能再改变其元素数量。如果要动态改变数组元素数量，只能通过合并或拆分的方法生成新的数组。对于刚上手NumPy的程序员来说，最大的困惑就是不能使用append( )函数向数组内添加元素，甚至都找不到append( )函数。其实，NumPy仍然保留了append( )函数，只不过这个函数不再是数组的函数，而是升级到最外层的NumPy命名空间了，并且该函数的功能不再是追加元素，而是合并数组，其代码如下。
>>> np.append([[1, 2, 3]], [[4, 5, 6]])
array([1, 2, 3, 4, 5, 6])
>>> np.append([[1, 2, 3]], [[4, 5, 6]], axis=0)
array([[1, 2, 3],
[4, 5, 6]])
>>> np.append([[1, 2, 3]], [[4, 5, 6]], axis=1)
array([[1, 2, 3, 4, 5, 6]])
不过，append( )函数还不够好用，推荐使用stack( )函数及其“同门小兄弟”：hstack( )水平合并函数、vstack( )垂直合并函数和dstack( )深度合并函数。下面演示这三个函数的用法。
>>> a = np.arange(4).reshape(2,2)
>>> b = np.arange(4,8).reshape(2,2)
>>> np.hstack((a,b)) # 水平合并
array([[0, 1, 4, 5],
[2, 3, 6, 7]])
>>> np.vstack((a,b)) # 垂直合并
array([[0, 1],
[2, 3],
[4, 5],
[6, 7]])
>>> np.dstack((a,b)) # 深度合并
array([[[0, 4],
[1, 5]],
[[2, 6],
[3, 7]]])
stack( )函数使用axis轴参数指定合并的规则，请仔细体会下面例子中axis轴参数的用法。
>>> a = np.arange(60).reshape(3,4,5)
>>> b = np.arange(60).reshape(3,4,5)
>>> a.shape, b.shape
>>> np.stack((a,b), axis=0).shape
(2, 3, 4, 5)
>>> np.stack((a,b), axis=1).shape
(3, 2, 4, 5)
>>> np.stack((a,b), axis=2).shape
(3, 4, 2, 5)
>>> np.stack((a,b), axis=3).shape
(3, 4, 5, 2)
拆分
因为数组切片非常简单，所以数组拆分应用较少。拆分是合并的逆过程，最常用的函数是split( )，其代码如下。
>>> a = np.arange(16).reshape(2,4,2)
>>> np.hsplit(a, 2) # 水平方向拆分成2部分
[array([[[ 0, 1], [ 2, 3]], [[ 8, 9], [10, 11]]]),
array([[[ 4, 5], [ 6, 7]], [[12, 13], [14, 15]]])]
>>> np.vsplit(a, 2) # 垂直方向拆分成2部分
[array([[[0, 1], [2, 3], [4, 5], [6, 7]]]),
array([[[ 8, 9], [10, 11], [12, 13], [14, 15]]])]
>>> np.dsplit(a, 2) # 深度方向拆分成2部分
[array([[[ 0], [ 2], [ 4], [ 6]], [[ 8], [10], [12], [14]]]),
array([[[ 1], [ 3], [ 5], [ 7]], [[ 9], [11], [13], [15]]])]
复制
改变数组结构返回的是原数组的一个新视图，而不是原数组的副本。浅复制（view）和深复制（copy）则是创建原数组的副本，但二者之间也有细微差别：浅复制（view）是共享内存，深复制（copy）是独享内存，其代码如下。
>>> a = np.arange(6).reshape((2,3))
>>> b = a.view()
>>> b is a
False
>>> b.base is a
False
>>> b.flags.owndata
False
>>> c = a.copy()
>>> c is a
False
>>> c.base is a
False
>>> c.flags.owndata
True
排序
NumPy数组有两个排序函数，一个是sort( )，另一个是argsort( )。sort( )函数返回输入数组的排序副本，argsort( )函数返回的是数组值从小到大的索引号。从函数原型看，这两个函数的参数完全一致。
numpy.sort(arr, axis=-1, kind='quicksort', order=None)
numpy.argsort(arr, axis=-1, kind='quicksort', order=None)
第1个参数arr，是要排序的数组；第2个参数axis，也就是轴，指定排序的轴，默认值-1表示没有指定排序轴，返回结果将沿着最后的轴排序；第3个参数kind，表示排序方法，默认为“quicksort”（快速排序），其他选项还有“mergesort”（归并排序）和“heapsort”（堆排序）；第4个参数order，指定用于排序的字段，前提是数组包含该字段。
>>> a = np.random.random((2,3))
>>> a
array([[0.79658569, 0.14507096, 0.63016223],
[0.24983103, 0.98368325, 0.71092079]])
>>> np.argsort(a) # 返回行内从小到大排序的索引号（列排序），相当于axis=1（最后的轴）
array([[1, 2, 0],
[0, 2, 1]], dtype=int64)
>>> np.sort(a) # 返回行内从小到大排序的一个新数组（列排序）
array([[0.14507096, 0.63016223, 0.79658569],
[0.24983103, 0.71092079, 0.98368325]])
>>> np.sort(a,axis=0) # 返回列内从小到大排序的一个新数组（行排序）
array([[0.24983103, 0.14507096, 0.63016223],
[0.79658569, 0.98368325, 0.71092079]])
下面演示的是排序字段的使用。先定义一个新的数据类型dt，它类似于一个字典，有姓名name和年龄age两个键值对，姓名的长度为10个字符，年龄的数据类型是整型。
>>> dt = np.dtype([('name', 'S10'),('age', int)])
>>> a = np.array([("zh",21),("wang",25),("li",17), ("zhao",27)], dtype = dt)
>>> np.sort(a, order='name') # 如果指定姓名排序，结果是李王张赵
array([(b'li', 17), (b'wang', 25), (b'zh', 21), (b'zhao', 27)],
dtype=[('name', 'S10'), ('age', '<i4')])
>>> np.sort(a, order='age') # 如果指定年龄排序，结果则是李张王赵
array([(b'li', 17), (b'zh', 21), (b'wang', 25), (b'zhao', 27)],
dtype=[('name', 'S10'), ('age', '<i4')])
查找
这里约定查找是返回数组中符合条件的元素的索引号，或返回和数组具有相同结构的布尔型数组，元素符合条件在布尔型数组中对应True，否则对应False。查找分为最大值和最小值查找、非零元素查找、使用逻辑表达式查找和使用where条件查找这4种方式。
- 最大值和最小值查找
下面的代码演示了返回数组中最大值和最小值的索引号，如果是多维数组，这个索引号是数组转成一维之后的索引号。
>>> a = np.random.random((2,3))
>>> a
array([[0.47881615, 0.55682904, 0.29173085],
[0.41107703, 0.91467593, 0.88852535]])
>>> np.argmax(a)
4
>>> np.argmin(a)
2
- 非零元素查找
下面的代码演示了返回数组中非零元素的索引号，返回的结果是一个元组。
>>> a = np.random.randint(0, 2, (2,3))
>>> a
array([[0, 0, 0],
[0, 1, 1]])
>>> np.nonzero(a)
(array([1, 1], dtype=int64), array([1, 2], dtype=int64))
- 使用逻辑表达式查找
下面的代码演示了使用逻辑表达式查找符合条件的元素，返回结果是一个和原数组结构相同的布尔型数组，元素符合条件在布尔型数组中对应True，否则对应False。
>>> a = np.arange(10).reshape((2,5))
>>> a
array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])
>>> (a>3)&(a<8)
array([[False, False, False, False, True],
[ True, True, True, False, False]])
- 使用where条件查找
np.where( )函数返回数组中满足给定条件的元素的索引号，其结构为元组，元组的第k个元素对应符合条件的元素在数组k轴上的索引号。这句话可以简单理解为，一维数组返回一个元素的元组，二维数组返回两个元素的元组，依此类推。np.where( )函数还可以用于替换符合条件的元素。
>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.where(a < 5)
(array([0, 1, 2, 3, 4], dtype=int64),)
>>> a = a.reshape((2, -1))
>>> a
array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])
>>> np.where(a < 5)
(array([0, 0, 0, 0, 0], dtype=int64), array([0, 1, 2, 3, 4], dtype=int64))
>>> np.where(a < 5, a, 10*a) # 满足条件的元素不变，其他元素乘10
array([[ 0, 1, 2, 3, 4],
[50, 60, 70, 80, 90]])
筛选
筛选是返回符合条件的元素。筛选条件有三种表示方式，一是使用np.where( )函数返回的Python元组，二是使用逻辑表达式返回的布尔型数组，三是使用整型数组，其代码如下。
>>> a = np.random.random((3,4))
>>> a
array([[0.41551063, 0.38984904, 0.01204226, 0.72323978],
[0.82425869, 0.64216573, 0.41475495, 0.21351508],
[0.30104819, 0.52046164, 0.58286043, 0.66749564]])
>>> a[np.where(a>0.5)] # 返回大于0.5的元素（使用np.where()函数返回的Python元组）
array([0.72323978, 0.82425869, 0.64216573, 0.52046164, 0.58286043, 0.66749564])
>>> a[(a>0.3)&(a<0.7)] # 返回大于0.3且小于0.7的元素（使用逻辑表达式返回的布尔型数组）
array([0.41551063, 0.38984904, 0.64216573, 0.41475495, 0.30104819,
0.52046164, 0.58286043, 0.66749564])
>>> a[np.array([2,1])] # 返回整型数组指定的项（使用整型数组）
array([[0.30104819, 0.52046164, 0.58286043, 0.66749564],
[0.82425869, 0.64216573, 0.41475495, 0.21351508]])
>>> a = a.ravel()
>>> a[np.array([3,5,7,11])] # 返回整型数组指定的项（使用整型数组）
array([0.72323978, 0.64216573, 0.21351508, 0.66749564])
>>> a[np.array([[3,5],[7,11]])] # 返回整型数组指定的项（使用整型数组）
array([[0.72323978, 0.64216573],
[0.21351508, 0.66749564]])
使用np.where( )函数或直接使用逻辑表达式来筛选数组元素，其目的是从数组中筛选出符合条件的元素。那么，使用整型数组来筛选数组元素的用途是什么呢？一个看似不起眼的功能，运用起来其实蕴含着无穷的想象空间。下面用一个例子来演示通过整型数组筛选数组元素的神奇魔法。图4-5所示的是用字符表现像素的灰度图。
一般而言，灰度图像每个像素的值域范围是[0, 255]。假如用于表现不同灰度的字符集是[’ ‘, ‘.’, ‘-‘, ‘+’, ‘=’, ‘*’, ‘#’, ‘@’]，从’ ‘到’@’表示从白到黑的8个灰度等级。我们需要将每个像素的灰度值分段转换成相应的字符。例如，灰度值小于32的像素用’@’表示，大于或等于32且小于64的像素用’#’表示，依次类推直至大于或等于224的像素用’ ‘表示。
如何实现图像数组从灰度值到对应字符的转换呢？乍一看，好像只有用循环的方式遍历所有像素才能实现。但是，下面的代码却用“整型数组筛选数组元素”的方法完成了这个转换，不但代码简洁，而且代码的执行速度也非常快。
>>> img = np.random.randint(0, 256, (5, 10), dtype=np.uint8)
>>> img
array([[145, 95, 60, 14, 66, 150, 221, 43, 184, 66],
[229, 138, 76, 90, 179, 217, 2, 20, 154, 191],
[165, 120, 77, 117, 42, 108, 156, 5, 208, 50],
[164, 196, 227, 111, 82, 84, 19, 208, 124, 16],
[146, 50, 107, 26, 34, 229, 137, 104, 93, 223]], dtype=uint8)
>>> img = (img/32).astype(np.uint8) # 将256级灰度值转为8级灰度值
>>> img
array([[0, 3, 3, 3, 1, 7, 7, 4, 2, 7],
[6, 3, 4, 5, 4, 5, 4, 7, 3, 4],
[6, 7, 1, 2, 2, 2, 2, 4, 7, 7],
[5, 7, 1, 2, 0, 2, 7, 0, 7, 5],
[3, 5, 0, 7, 0, 4, 6, 2, 5, 0]], dtype=uint8)
>>> chs = np.array([' ', '.', '-', '+', '=', '*', '#', '@']) # 灰度字符集
>>> chs[img] # 用整型数组筛选数组元素（我认为这是NumPy最精彩之处！）
array([[' ', '+', '+', '+', '.', '@', '@', '=', '-', '@'],
['#', '+', '=', '*', '=', '*', '=', '@', '+', '='],
['#', '@', '.', '-', '-', '-', '-', '=', '@', '@'],
['*', '@', '.', '-', ' ', '-', '@', ' ', '@', '*'],
['+', '*', ' ', '@', ' ', '=', '#', '-', '*', ' ']], dtype='<U1')
数组I/O
数组I/O就是讨论如何分发、交换数据。在机器学习算法模型的例子中，海量的训练数据通常都是从数据文件中读出来的，而数据文件一般是CSV格式。NumPy自带CSV格式文件读写函数，可以很方便地读写CSV格式的数据文件。除了支持通用的CSV格式的数据文件，NumPy还为数组对象引入了两个新的二进制文件格式，用于数据交换。扩展名为.npy的数据文件用于存储单个数组，扩展名为.npz的数据文件用于存储多个数组。NumPy支持的数据文件格式及读写函数的对应关系具体如图4-6所示。
CSV是一种通用的、相对简单的文件格式。CSV格式的数据文件以纯文本形式存储表格数据，由任意数目的记录组成，记录间以某种换行符分隔。每条记录由字段组成，字段间的分隔符是其他字符或字符串，最常见的是逗号或制表符。通常，所有记录都有完全相同的字段序列。下面的代码演示了NumPy读写CSV格式的数据文件的方法。实际操作下面的代码时，请注意结合实际情况替换对应的文件路径和文件名。
>>> a = np.random.random((15,5))
>>> np.savetxt('demo.csv', a, delimiter=',') # 将数组a保存成CSV格式的数据文件
>>> data = np.loadtxt('demo.csv', delimiter=',') # 打开CSV格式的数据文件
>>> data.shape, data.dtype
((15, 5), dtype('float64'))
NumPy自定义的数据交换格式也是一个非常好用的数据交换方式，使用它保存NumPy数组时不会丢失任何信息，特别是数据类型的信息。实际操作下面的代码时，请注意结合实际情况替换对应的文件路径和文件名。
>>> single_arr_fn = 'single_arr.npy' # 存储单个数组文件名
>>> multi_arr_fn = 'multi_arr.npz' # 存储多个数组文件名
>>> lon = np.linspace(10,90,9)
>>> lat = np.linspace(20,60,5)
>>> np.save(single_arr_fn, lon) # 用save()函数把经度数组保存成.npy文件
>>> lon = np.load(single_arr_fn) # 接着用load()函数读出来
>>> np.savez(multi_arr_fn, longitude=lon, latitude=lat) #保存两个数组到一个文件
>>> data = np.load(multi_arr_fn) # 用load()函数把这个.npz文件读成一个结构data
>>> data.files # 查看所有的数组名
>>> data['longitude'] # 使用data[数组名]，就可以取得想要的数据
>>> data['latitude'] # 使用data[数组名]，就可以取得想要的数据
常用函数
NumPy提供了大量的函数，其中和数组操作相关的函数有很大一部分是数组对象ndarray的方法在NumPy命名空间中的映射。也就是说，NumPy的很多函数其实是数组对象内置的方法。本节并不是对所有NumPy函数的总结，而是介绍和函数相关的几个概念，并对部分常用函数做示例性讲解。
常量
说起常量，我们首先会想到圆周率、自然常数、欧拉常数等。的确，NumPy的常量包括np.pi（圆周率）、np.e（自然常数）、np.euler_gamma（欧拉常数），此外还包括np.nan（非数字）和np.inf（无穷大）这两个特殊值。NumPy的特殊值还有正负无穷大、正负零等，但因为很少用到，这里就不进行重点介绍。
NumPy有两个很有趣的特殊值：np.nan和np.inf。nan是Not a Number的简写，意为非数字；inf是infinity的简写，意为无穷大。其代码如下。
>>> a = np.array([1, 2, np.nan, np.inf])
>>> a.dtype
dtype('float64')160
>>> a[0] = np.nan
>>> a[1] = np.inf
>>> a
array([nan, inf, nan, inf])
>>> a[0] == a[2] # 两个np.nan不相等
False
>>> a[1] == a[3] # 两个np.inf则相等
True
>>> np.isnan(a[0]) # 判断一个数组元素是否是np.nan
True
>>> np.isinf(a[1]) # 判断一个数组元素是否是np.inf
True
上面的代码中，两个np.nan居然不相等，但两个np.inf则是相等的。另外请注意，判断一个数组元素是否是np.nan或np.inf，需要使用np.isnan( )和np.isinf( )这两个相应的函数，而不是使用两个等号的逻辑表达式。
那么这两个特殊值有什么用途呢？原来，NumPy是用特殊值来表示缺值、空值和无效值的。想一想，Python语言和C语言如何表示缺值、空值和无效值呢？Python语言因为列表元素不受类型限制，可以用None或False等表示缺值、空值和无效值。而C语言只能在数据的值域范围之外，选一个特定值来表示。例如，假定数组存储的是学生的成绩，成绩一般都是正值，所以C语言可以用-1表示缺考。在NumPy数组中，因为有nan和inf这两个特殊值，所以不用在意数据的值域范围。说到这里，你可能会产生疑问，这两个特殊值，一个是非数字，一个是无穷大，数组运算的时候怎么处理呢？这就是NumPy神奇的地方，我们根本不用担心这个问题。
下面的代码演示了在数组相邻的两个元素之间插入它们的算术平均值。尽管数组元素包含np.nan，但这不影响数值计算。
>>> a = np.array([9, 3, np.nan, 5, 3])
>>> a = np.repeat(a,2)[:-1]
>>> a[1::2] += (a[2::2]-a[1::2])/2
>>> a
array([ 9., 6., 3., nan, nan, nan, 5., 4., 3.])
命名空间
刚开始使用NumPy函数时，对于函数的使用，你可能会有这样的困惑：实现同样的功能，一个函数却有两种写法；有时以为某个函数可以有两种写法，但用起来却会出错。归纳起来，这些困惑有以下三种类型。
- 都是求和、求极值，下面这两种写法有什么区别吗？
>>> a = np.random.random(10) >>> a.max(), np.max(a) (0.8975052328686041, 0.8975052328686041) >>> a.sum(), np.sum(a) (5.255303938070831, 5.255303938070831) - 同样是复制，为什么深复制copy( )两种写法都行，而浅复制view( )则只有数组的方法？
>>> a = np.random.random(5) >>> a.copy() array([0.59646094, 0.99280395, 0.1046394 , 0.11498018, 0.17936631]) >>> np.copy(a) array([0.59646094, 0.99280395, 0.1046394 , 0.11498018, 0.17936631]) >>> a.view() array([0.59646094, 0.99280395, 0.1046394 , 0.11498018, 0.17936631]) >>> np.view(a) Traceback (most recent call last): File "<pyshell#61>", line 1, in <module> np.view(a) AttributeError: module 'numpy' has no attribute 'view' - 为什么where( )不能作为数组ndarray的函数，必须作为NumPy的函数？
>>> np.where(a>0.5) (array([3, 4, 5, 8, 9], dtype=int64),) >>> a.where(a>0.5) Traceback (most recent call last): File "<pyshell#65>", line 1, in <module> a.where(a>0.5) AttributeError: 'numpy.ndarray' object has no attribute 'where'
以上这些差异取决于函数在不同的命名空间是否有映射。数组的大部分函数在顶层命名空间有映射，因此可以有两种写法。但数组的一小部分函数没有映射到顶层命名空间，所以只能有一种写法。而顶层命名空间的大部分函数，也都只有一种写法。表4-2所示的是常用函数和命名空间的关系，仅供参考。
数学函数
如果不熟悉NumPy，Python程序员一般都会选择使用math模块来解决数学问题。但实际上NumPy的数学函数比math模块更加方便，而且NumPy的数学函数可以广播到数组的每一个元素上，也就是说，如果用np.sqrt( )对数组arr开方，返回的是数组arr中每个元素的平方根组成的新数组。
下面把NumPy和math模块的数学函数罗列在一个表格中，分成了数学常数、舍入函数、快速转换函数、幂指数对数函数和三角函数这5类，如表4-3所示。其他如求和、求差、求积的函数被归类到下一小节的统计函数中。
下面的代码演示的是一些常用数学函数。
>>> import numpy as np
>>> import math
>>> math.e == np.e # 两个模块的自然常数相等
True
>>> math.pi == np.pi # 两个模块的圆周率相等
True
>>> np.ceil(5.3)
6.0
>>> np.ceil(-5.3)
-5.0
>>> np.floor(5.8)
5.0
>>> np.floor(-5.8)
-6.0
>>> np.around(5.87, 1)
5.9
>>> np.rint(5.87)
6.0
>>> np.degrees(np.pi/2)
90.0
>>> np.radians(180)
3.141592653589793
>>> np.hypot(3,4) # 求平面上任意两点的距离
5.0
>>> np.power(3,1/2)
1.7320508075688772
>>> np.log2(1024)
10.0
>>> np.exp(1)
2.718281828459045
>>> np.sin(np.radians(30)) #正弦、余弦函数的周期是2π
0.49999999999999994
>>> np.sin(np.radians(150))
0.49999999999999994
>>> np.degrees(np.arcsin(0.5)) # 反正弦、反余弦函数的周期则是π
30.000000000000004
上述代码中使用的函数参数都是单一的数值，实际上，这些函数也都可以用到NumPy数组上。例如，平面直角坐标系中有1000万个点，它们的x坐标和y坐标都分布在[0,1)区间，哪一个点距离点(0.5,0.5)最近呢？使用NumPy数组计算的代码如下。
>>> p = np.random.random((10000000,2))
>>> x, y = np.hsplit(p,2) # 分离每一个点的x和y坐标
>>> d = np.hypot(x-0.5,y-0.5) # 计算每一个点距离点(0.5,0.5)的距离
>>> i = np.argmin(d) # 返回最短距离的点的索引号
>>> print('距离点(0.5,0.5)最近的点的坐标是(%f, %f)，距离为%f'%(*p[i], d[i]))
距离点(0.5,0.5)最近的点的坐标是(0.499855, 0.499877)，距离为0.000190
统计函数
NumPy的统计函数有很多，并且很多函数还同时提供了忽略nan（缺值或无效值）的形式。常用的统计函数大致上可以分成查找特殊值、求和差积、均值和方差以及相关系数这4类，详细说明如表4-4所示。
在实际编程实践中，我们所获得的数据远没有想象中的那么理想，存在缺值或无效值是常态。假定用np.nan表示无效值，一旦数据中存在无效值，对一个函数而言，是否忽略无效值将会得到完全不同的结果。下面先以求最大值和最小值为例，演示忽略np.nan的必要性。
>>> a = np.random.random(10)
>>> np.max(a), np.min(a)
(0.9690291560970601, 0.19240165472992765)
>>> a[1::2] = np.nan # 将索引号为1、3、5、7、9的元素设置为np.nan
>>> a
array([0.80138474, nan, 0.8615121 , nan, 0.19240165,
nan, 0.61915229, nan, 0.96902916, nan])
>>> np.max(a), np.min(a) # 此时，min()函数和max()函数失效了
(nan, nan)
>>> np.nanmax(a), np.nanmin(a) # 必须使用nanmax()函数和nanmin()函数
(0.9690291560970601, 0.19240165472992765)
方差和标准差是衡量数据离散程度最重要且最常用的指标，也是统计学上最重要的分析工具和手段。方差是各个数据与其算术平均值的离差平方和的平均值。方差的算术平方根，即为标准差。统计学上，方差和标准差使用比较频繁，下面来演示一下。
>>> a = np.random.randint(0,50,(3,4))
>>> np.sum(np.square(a-a.mean()))/a.size # 用方差定义求方差
238.25
>>> np.var(a) # 直接用方差函数求方差，与用方差定义求方差的结果相同
238.25
>>> np.sqrt(a.var()) # 对方差求算术平方根就是标准差
15.435349040433131
>>> a.std() # 直接用标准差函数求标准差，与对方差求算术平方根的结果相同
15.435349040433131
下面的例子综合运用统计函数，来分析两只股票的关联关系和收益率。pa和pb是两只股票连续30个交易日的股价数组。每日股价收益率定义为当日股价与前一个交易日股价之差再除以最后一个交易日的股价。
>>> pa = np.array([79.66, 81.29, 80.37, 79.31, 79.84, 78.53, 78.29, 78.51, 77.99,
79.82, 80.41, 79.27, 80.26, 81.61, 81.39, 80.29, 80.18, 78.38, 75.06, 76.15, 75.66,
73.90, 72.14, 74.27, 75.27, 76.15, 75.40, 76.51, 77.57, 77.06])
>>> pb = np.array([30.93, 31.61, 31.62, 31.77, 32.01, 31.52, 30.09, 30.54, 30.78,
30.84, 30.80, 30.38, 30.88, 31.38, 31.05, 29.90, 29.96, 29.59, 28.71, 28.95, 29.19,
28.71, 27.93, 28.35, 28.92, 29.17, 29.02, 29.43, 29.12, 29.11])
>>> np.corrcoef(pa, pb) # 两只股票的相关系数为0.867，关联比较密切
array([[1. , 0.86746674],
[0.86746674, 1. ]])
>>> pa_re = np.diff(pa)/pa[:-1] # 股价收益率
>>> pb_re = np.diff(pb)/pb[:-1] # 股价与前一个交易日股价之差再除以最后一个交易日的股价
>>> import matplotlib.pyplot as plt
>>> plt.plot(pa_re)
[<matplotlib.lines.Line2D object at 0x000002262AEBD9C8>]
>>> plt.plot(pb_re)
[<matplotlib.lines.Line2D object at 0x000002262BB96408>]
>>> plt.show()
插值函数
数据插值是数据处理过程中经常用到的技术，常用的插值有一维插值、二维插值、高阶插值等，常见的算法有线性插值、B样条插值、临近插值等。不过，NumPy只提供了一个简单的一维线性插值函数np.interp( )，其他更加复杂的插值功能放到了SciPy中，具体讲解详见本书7.2节。
下面用一个实例来演示NumPy一维线性插值函数的使用方法。假定_x和_y是原始的样本点x坐标和y坐标构成的数组，总数只有11个点。如果想在_x的值域范围内插值更多的点，如增加到33个点，就需要在_x的值域范围内生成33个点的x坐标构成的数组x，再利用插值函数np.interp( )得到对应的33个点的y坐标构成的数组y。为了更直观地理解线性插值，这里提前使用了第5章介绍的数据可视化方法。
>>> import matplotlib.pyplot as plt
>>> _x = np.linspace(0, 2*np.pi, 11)
>>> _y = np.sin(_x)
>>> x = np.linspace(0, 2*np.pi, 33)
>>> y = np.interp(x, _x, _y)
>>> plt.plot(x, y, 'o')
[<matplotlib.lines.Line2D object at 0x0000020A2A8D1048>]
>>> plt.plot(_x, _y, 'o')
[<matplotlib.lines.Line2D object at 0x0000020A2A5ED148>]
>>> plt.show()
橘黄色（浅色）的点是原始样本点，蓝色（深色）的点是进行一维线性插值后的点。
多项式拟合函数
拟合与插值看起来有一些相似，所以初学者比较容易混淆，实际上二者是完全不同的概念。拟合又称回归，是指已知某函数的若干离散函数值，通过调整该函数中若干待定系数，使得该函数与已知离散函数值的误差达到最小。
多项式拟合是最常见的拟合方法。对函数f(x)，我们可以使用一个k阶多项式去近似。
f(x)≈g(x)=a0+a1x+a2x2+a3x3+…+axxk
通过选择合适的系数（最佳系数），可以使函数f(x)和g(x)之间的误差达到最小。最小二乘法被用于寻找多项式的最佳系数。NumPy提供了一个非常简单易用的多项式拟合函数np.polyfit( )，只需要输入一组自变量的离散值，和一组对应的函数f(x)，并指定多项式次数k，就可以返回一组最佳系数。函数np.poly1d( )，则可以将一组最佳系数转成函数g(x)。
下面的例子首先生成了原始数据点_x和_y，然后分别用4次、5次、6次和7次多项式去拟合原始数据，并计算出每次拟合的误差。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> plt.rcParams['font.sans-serif'] = ['FangSong'] # 指定字体以保证中文正常显示
>>> plt.rcParams['axes.unicode_minus'] = False # 正确显示连字符
>>> _x = np.linspace(-1, 1, 201)
>>> _y = ((_x**2-1)**3 + 0.5)*np.sin(2*_x) + np.random.random(201)/10 - 0.1
>>> plt.plot(_x, _y, ls='', marker='o', label="原始数据")
[<matplotlib.lines.Line2D object at 0x0000011D87DFFC08>]
>>> for k in range(4, 8):
g = np.poly1d(np.polyfit(_x, _y, k)) # g是k次多项式
loss = np.sum(np.square(g(_x)-_y)) # g(x)和f(x)的误差
plt.plot(_x, g(_x), label="%d次多项式，误差：%0.3f"%(k,loss))
[<matplotlib.lines.Line2D object at 0x0000011D87FDAA08>]
[<matplotlib.lines.Line2D object at 0x0000011D8A2A9A88>]
[<matplotlib.lines.Line2D object at 0x0000011D87FE0AC8>]
[<matplotlib.lines.Line2D object at 0x0000011D87FE5348>]
>>> plt.legend()
<matplotlib.legend.Legend object at 0x0000011D87FE00C8>
>>> plt.show()
分别用4次、5次、6次和7次多项式拟合的效果，可以看出7次多项式拟合的误差最小。
自定义广播函数
前面讲过，广播是NumPy最具特色的特性之一，几乎所有的NumPy函数都可以通过广播特性将操作映射到数组的每一个元素上。然而NumPy函数并不能完成所有的工作，有些工作还需要我们自己来定义函数。如何让我们自己定义的函数也可以广播到数组的每一个元素上，就是自定义广播函数要做的事情。
假定a和b是两个结构相同的数组，数据类型为8位无符号整型。我们希望用a和b生成一个具有相同结构的新数组c，生成规则是：同一位置的元素，若a或b任何一方等于0，则c等于0；若a等于b，则c等于0；若a和b均等于2的整数次幂，则c取a和b中的较大者；若a或b只有一方等于2的整数幂，则c等于满足条件的一方；以上条件都不满足时，c取a和b中的较大者。
很显然，现有的NumPy函数都无法实现这个计算功能，因此需要自定义函数，其代码如下。
>>> def func_demo(x, y):
if x == 0 or y == 0 or x == y:
return 0
elif x&(x-1) == 0 and y&(y-1) == 0: # x和y都是2的整数次幂
return max(x, y)
elif x&(x-1) == 0: # 仅有x等于2的整数次幂
return x
elif y&(y-1) == 0: # 仅有y等于2的整数次幂
return y
else:
return max(x, y)
将自定义函数变成广播函数的方法有两个，下面分别进行详细讲解。
- 使用np.frompyfunc( )定义广播函数
使用np.frompyfunc( )将数值函数转换成数组函数需要提供三个参数：数值函数、输入参数的个数和返回值的个数。另外，np.frompyfunc( )返回的广播函数，其返回值是object类型，最终需要根据实际情况显式地转换数据类型，其代码如下。
>>> uf = np.frompyfunc(func_demo, 2, 1)
>>> a = np.random.randint(0, 256, (2,5), dtype=np.uint8)
>>> b = np.random.randint(0, 256, (2,5), dtype=np.uint8)
>>> a
array([[118, 33, 164, 187, 48],
[ 41, 128, 242, 225, 34]], dtype=uint8)
>>> b
array([[170, 207, 35, 61, 251],
[251, 206, 70, 208, 85]], dtype=uint8)
>>> c = uf(a, b)
>>> c
array([[170, 207, 164, 187, 251],
[251, 128, 242, 225, 85]], dtype=object)
>>> c = c.astype(np.uint8)
>>> c
array([[170, 207, 164, 187, 251],
[251, 128, 242, 225, 85]], dtype=uint8)
- 使用np.vectorize( )定义广播函数
np.frompyfunc( )适用于多个返回值的函数。如果返回值只有一个，使用np.vectorize( )定义广播函数更方便，并且还可以通过otypes参数指定返回数组的元素类型，其代码如下。
>>> uf = np.vectorize(func_demo, otypes=[np.uint8])
>>> a = np.array([[118, 33, 164, 187, 48],
[ 41, 128, 242, 225, 34]], dtype=np.uint8)
>>> b = np.array([[170, 207, 35, 61, 251],
[251, 206, 70, 208, 85]], dtype=np.uint8)
>>> c = uf(a, b)
>>> c
array([[170, 207, 164, 187, 251],
[251, 128, 242, 225, 85]], dtype=uint8)
自定义广播函数并不是真正的广播函数，其运行效率和循环遍历几乎没有差别，因此除非确实必要，否则不应该滥用自定义广播函数。事实上，总有一些技巧可以不用遍历数组也能实现对数组元素的操作，如对数组元素分组操作等。
掩码数组
在科研活动和实际工作中，我们获得的数据集往往是有缺失或被污染的，如卫星上各种载荷的传感器在某一瞬间甚至某一段时间内可能无法记录数据或记录值被干扰。上一节简单介绍了NumPy处理数据缺值或无效值的思路，本节则是针对这个问题的完整的实现方案。
numpy.ma子模块通过引入掩码数组提供了一种解决数据缺失或无效问题的安全、便捷的方法。numpy.ma子模块的主体是MaskedArray类，它是numpy.ndarray的派生类，可以把numpy.ma子模块当作ndarray来用，且无须考虑数组的无效值是否会给操作带来无法预知的意外。
创建掩码数组
这里约定导入掩码数组子模块的方法和导入NumPy模块的风格保持一致。
>>> import numpy as np
>>> import numpy.ma as ma
- 由列表生成掩码数组
掩码数组子模块的ma.array( )函数和NumPy的np.array( )函数类似，可以直接将列表生成掩码数组，默认mask参数为False，生成的数组类型是MaskedArray类。数组掩码处理后，无论是查找最大值、最小值，还是计算均值、方差，都不用再担心数据是否无效的问题了。
>>> a = ma.array([0, 1, 2, 3], mask=[0, 0, 1, 0]) # 指定第三个元素无效
>>> a
masked_array(data=[0, 1, --, 3],
mask=[False, False, True, False],
fill_value=999999)
>>> type(a)
<class 'numpy.ma.core.MaskedArray'>
>>> a.min(), a.max(), a.mean(), a.var()
(0, 3, 1.3333333333333333, 1.5555555555555556)
- 由数组生成掩码数组
ma.asarray( )函数可以将普通的NumPy数组转成掩码数组。新生成的掩码数组不会对原数组中的np.nan或np.inf做掩码处理，但是会相应调整填充值（fill_value）。
>>> a = np.arange(5)
>>> ma.asarray(a)
masked_array(data=[0, 1, 2, 3, 4],
mask=False,
fill_value=999999)
>>> a = np.array([1, np.nan, 2, np.inf, 3]) # 包含特殊值的数组
>>> ma.asarray(a)
masked_array(data=[ 1., nan, 2., inf, 3.],
mask=False,
fill_value=1e+20) # 填充值会相应变化
- 对数组中的无效值做掩码处理
ma.asarray( )函数不会对原数组中的np.nan或np.inf做掩码处理，ma.masked_invalid( )函数则可以实现这个功能。
>>> a = np.array([1, np.nan, 2, np.inf, 3])
>>> ma.masked_invalid(a)
masked_array(data=[1.0, --, 2.0, --, 3.0],
mask=[False, True, False, True, False],
fill_value=1e+20)
- 对数组中的给定值做掩码处理
有时需要将数组中的某个给定值设置为无效（掩码），ma.masked_equal( )函数可以实现这个功能。
>>> a = np.arange(3).repeat(2)
>>> ma.masked_equal(a, 1) # 对数组元素1做掩码处理
masked_array(data=[0, 0, --, --, 2, 2],
mask=[False, False, True, True, False, False],
fill_value=1)
- 对数组中符合条件的特定值做掩码处理
有时需要将数组中符合条件的某些特定值设置为无效（掩码），掩码数组子模块提供了若干函数实现条件掩码。这些可能的筛选条件包括大于、大于等于、小于、小于等于、区间内、区间外等6种。下面的代码演示了6种筛选条件对应的6个掩码函数的使用方法。
>>> a = np.arange(8)
>>> ma.masked_greater(a, 4) # 掩码大于4的元素
masked_array(data=[0, 1, 2, 3, 4, --, --, --],
mask=[False, False, False, False, False, True, True, True],
fill_value=999999)
>>> ma.masked_greater_equal(a, 4) # 掩码大于等于4的元素
masked_array(data=[0, 1, 2, 3, --, --, --, --],
mask=[False, False, False, False, True, True, True, True],
fill_value=999999)
>>> ma.masked_less(a, 4) # 掩码小于4的元素
masked_array(data=[--, --, --, --, 4, 5, 6, 7],
mask=[ True, True, True, True, False, False, False, False],
fill_value=999999)
>>> ma.masked_less_equal(a, 4) # 掩码小于等于4的元素
masked_array(data=[--, --, --, --, --, 5, 6, 7],
mask=[ True, True, True, True, True, False, False, False],
fill_value=999999)
>>> ma.masked_inside(a, 2, 5) # 掩码在 [2,5]内的元素
masked_array(data=[0, 1, --, --, --, --, 6, 7],
mask=[False, False, True, True, True, True, False, False],
fill_value=999999)
>>> ma.masked_outside(a, 2, 5) # 掩码在 [2,5]之外的元素
masked_array(data=[--, --, 2, 3, 4, 5, --, --],
mask=[ True, True, False, False, False, False, True, True],
fill_value=999999)
- 用一个数组的条件筛选结果对另一个数组做掩码处理
a和b是两个结构相同的数组，如果用a>5的条件对数组b掩码，上面那些函数就失效了。这种情况正是ma.masked_where( )函数可以大显身手的时候。当然，该函数也可以对数组自身掩码。
>>> a = np.arange(8)
>>> b = np.random.random(8)
>>> ma.masked_where(a>5, b) # 用a>5的条件对数组b掩码
masked_array(data=[0.08445100592764732, 0.3502664957826195,
0.38403008762851243, 0.13025516166581663,
0.393489225584158, 0.6623482125806512, --, --],
mask=[False, False, False, False, False, False, True, True],
fill_value=1e+20)
访问掩码数组
- 索引和切片
因为掩码数组MaskedArray类是numpy.ndarray的派生类，所以那些用在普通NumPy数组上的索引和切片操作也依然有效。
>>> a = np.array([1, np.nan, 2, np.inf, 3])
>>> a = ma.masked_invalid(a)
>>> a[0], a[1], a[-1]
(1.0, masked, 3.0)
>>> a[1:-1]
masked_array(data=[--, 2.0, --],
mask=[ True, False, True],
fill_value=1e+20)
- 函数应用
掩码数组内置方法的使用和普通数组没有区别，下面的代码演示了最大值、最小值、均值和方差的使用。
使用NumPy命名空间的函数则要慎重，如果掩码数组子模块有对应函数，应优先使用掩码数组子模块的对应函数。例如，对掩码数组求正弦，如果使用np.sin( )函数，会发出警告信息；如果使用ma.sin( )函数，则无任何问题。
>>> a = np.array([1, np.nan, 2, np.inf, 3])
>>> a = ma.masked_invalid(a)
>>> a.min(), a.max(), a.mean(), a.var()
(1.0, 3.0, 2.0, 0.6666666666666666)
>>> np.sin(a)
Warning (from warnings module):
File "__main__", line 1
RuntimeWarning: invalid value encountered in sin
masked_array(data=[0.8414709848078965, --, 0.9092974268256817, --,
0.1411200080598672],
mask=[False, True, False, True, False],
fill_value=1e+20)
>>> ma.sin(a)
masked_array(data=[0.8414709848078965, --, 0.9092974268256817, --,
0.1411200080598672],
mask=[False, True, False, True, False],
fill_value=1e+20)
- 掩码数组转为普通数组
任何情况下，我们都可以通过掩码数组的data属性来获得掩码数组的数据视图，其类型就是np.ndarray数组。另外，还可以使用掩码数组的__array__( )函数或ma.getdata( )函数来获取掩码数组的数据视图。上述三种方法获得数据视图的操作，本质上都是操作掩码数组本身。如果需要数据视图副本，需使用copy( )函数。
>>> a = ma.array([1, np.nan, 2, np.inf, 3])
>>> a
masked_array(data=[ 1., nan, 2., inf, 9.],
mask=False,
fill_value=1e+20)
>>> x = a.data
>>> y = a.__array__()
>>> z = ma.getdata(a)
>>> w = np.copy(a.__array__()) # 复制数据视图
>>> x
array([ 1., nan, 2., inf, 3.])
>>> y
array([ 1., nan, 2., inf, 3.])
>>> z
array([ 1., nan, 2., inf, 3.])
>>> w
array([ 1., nan, 2., inf, 3.])
>>> a[-1] = 9
>>> x
array([ 1., nan, 2., inf, 9.])
>>> y
array([ 1., nan, 2., inf, 9.])
>>> z
array([ 1., nan, 2., inf, 9.])
>>> w
array([ 1., nan, 2., inf, 3.])
- 修改掩码
通过掩码数组的mask属性可以查看当前数组的掩码情况，其代码如下。通常，数组的掩码是一个布尔型数组，或是一个布尔值。
>>> a = ma.masked_invalid(np.array([1, np.nan, 2, np.inf, 3]))
>>> a.mask
array([False, True, False, True, False])
如果要对数组切片掩码或对数组的某个元素掩码，直接令该切片或该元素等于ma.masked常量即可，其代码如下。
>>> a = ma.masked_invalid(np.array([1, np.nan, 2, np.inf, 3]))
>>> a
masked_array(data=[1.0, --, 2.0, --, 3.0],
mask=[False, True, False, True, False],
fill_value=1e+20)
>>> a[:2] = ma.masked
>>> a
masked_array(data=[--, --, 2.0, --, 3.0],
mask=[ True, True, False, True, False],
fill_value=1e+20)
如果要撤销对数组切片或数组中的某个元素的掩码，只需要对该切片或该元素做赋值操作即可，其代码如下。
>>> a = ma.masked_invalid(np.array([1, np.nan, 2, np.inf, 3]))
>>> a
masked_array(data=[1.0, --, 2.0, --, 3.0],
mask=[False, True, False, True, False],
fill_value=1e+20)
>>> a[1] = 1.5
>>> a[2:4] = 5
>>> a
masked_array(data=[1.0, 1.5, 5.0, 5.0, 3.0],
mask=[False, False, False, False, False],
fill_value=1e+20)
矩阵对象
在数学上，矩阵（Matrix）是一个按照矩形阵列排列的复数或实数集合，但在NumPy中，矩阵np.matrix是数组np.ndarray的派生类。这意味着矩阵本质上是一个数组，拥有数组的所有属性和方法；同时，矩阵又有一些不同于数组的特性和方法。
首先，矩阵是二维的，不能像数组一样幻化成任意维度，即使展开或切片，返回也是二维的；其次，矩阵和矩阵、矩阵和数组都可以做加减乘除运算，运算结果总是返回矩阵；最后，矩阵的乘法不同于数组乘法。
创建矩阵
np.mat( )函数用于创建矩阵，它可以接受列表、数组甚至是字符串等形式的参数，还可以使用dtype参数指定数据类型，其代码如下。
>>> np.mat([[1,2,3],[4,5,6]], dtype=np.int) # 使用列表创建矩阵
matrix([[1, 2, 3],
[4, 5, 6]])
>>> np.mat(np.arange(6).reshape((2,3))) # 使用数组创建矩阵
matrix([[0, 1, 2],
[3, 4, 5]])
>>> np.mat('1 4 7; 2 5 8; 3 6 9') # 使用MATLAB风格的字符串创建矩阵
matrix([[1, 4, 7],
[2, 5, 8],
[3, 6, 9]])
此外，和生成特殊值数组类似，numpy.matlib子模块也提供了多个函数用于生成特殊值矩阵和随机数矩阵，其代码如下。
>>> import numpy.matlib as mat
>>> mat.zeros((2,3)) # 全0矩阵
matrix([[0., 0., 0.],
[0., 0., 0.]])
>>> mat.ones((2,3)) # 全1矩阵
matrix([[1., 1., 1.],
[1., 1., 1.]])
>>> mat.eye(3) # 单位矩阵
matrix([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
>>> mat.empty((2,3)) # 空矩阵
matrix([[1., 1., 1.],
[1., 1., 1.]])
>>> mat.rand((2,3)) # [0,1)区间的随机数矩阵
matrix([[0.99553367, 0.53978781, 0.96430715],
[0.02990035, 0.19574645, 0.2558724 ]])
>>> mat.randn((2,3)) # 均值为0、方差为1的高斯（正态）分布矩阵
matrix([[ 0.44248199, -1.51660328, 0.75149411],
[-0.8000728 , -0.4941788 , -0.33822774]])
矩阵特有属性
矩阵有几个特有的属性，如转置矩阵、逆矩阵、共轭矩阵、共轭转置矩阵等。熟悉这些属性对矩阵计算会有很大的帮助。
>>> m = np.mat(np.arange(6).reshape((2,3)))
>>> m
matrix([[0, 1, 2],
[3, 4, 5]])
>>> m.T # 返回自身的转置矩阵
matrix([[0, 3],
[1, 4],
[2, 5]])
>>> m.H # 返回自身的共轭转置矩阵
matrix([[0, 3],
[1, 4],
[2, 5]])
>>> m.I # 返回自身的逆矩阵
matrix([[-0.77777778, 0.27777778],
[-0.11111111, 0.11111111],
[ 0.55555556, -0.05555556]])
>>> m.A # 返回自身数据的视图（ndarray类）
array([[0, 1, 2],
[3, 4, 5]])
矩阵乘法
矩阵运算和数组运算大致相同，只有乘法运算有较大差别。在讲广播和矢量化时，我们已经知道，两个数组相乘就是对应元素相乘，条件是两个数组的结构相同。事实上，即使两个数组的结构不同，只要满足特定条件，也能做乘法运算。
>>> a = np.random.randint(0,10,(2,3))
>>> a
array([[7, 2, 3],
[1, 2, 7]])
>>> b = np.random.randint(0,10,3)
>>> b
array([9, 6, 2])
>>> a*b # 结构不同的两个数组也可以相乘
array([[63, 12, 6],
[ 9, 12, 14]])
>>> b*a
array([[63, 12, 6],
[ 9, 12, 14]])
除了对应元素相乘，数组还可以使用np.dot( )函数相乘，其代码如下。
>>> a = np.random.randint(0,10,(2,3))
>>> b = np.random.randint(0,10,3)
>>> c = np.random.randint(0,10,(3,2))
>>> np.dot(a,b)
array([38, 46])
>>> np.dot(a,c)
array([[44, 35],
[45, 40]])
对于数组而言，使用星号相乘和使用np.dot( )函数相乘是完全不同的两种乘法；对于矩阵来说，不管是使用星号相乘还是使用np.dot( )函数相乘，结果都是np.dot( )函数相乘的结果，因为矩阵没有对应元素相乘这个概念。np.dot( ) 函数实现的乘法就是矩阵乘法。那么矩阵乘法究竟是怎么运算的呢？图4-10显示了矩阵相乘的具体算法。
不是所有的矩阵都能相乘。图4-10中矩阵A乘矩阵B，二者可以相乘的条件是矩阵A的列数必须等于矩阵B的行数。例如，矩阵A是4行2列，矩阵B是2行3列，A×B没问题，但是反过来，B×A就无法运算了。可见，矩阵乘法不满足交换律。再来看一看乘法规则。概括来说，就是矩阵A的各行逐一去乘矩阵B的各列。例如，矩阵A的第1行和矩阵B的第1列，它们的元素个数一定相等，对应元素相乘后求和的值作为结果矩阵第1行第1列的值。又如，矩阵A的第3行和矩阵B的第3列，对应元素相乘后求和的值作为结果矩阵第3行第3列的值。依此类推，最终得到矩阵A乘矩阵B的结果矩阵。
那么，这个令人眼花缭乱的矩阵乘法有什么使用价值吗？答案是有，不但有，而且有非常大的使用价值。对于编程来说，矩阵乘法最常见的应用是图像的旋转、位移、透视等操作。图4-11显示了一个平面直角坐标系旋转矩阵的推导过程。
下面应用这个推导结果定义一个函数，来返回平面上的点围绕原点旋转给定角度后的坐标。
>>> def rotate(p,d):
a = np.radians(d)
m = np.array([[np.cos(a), np.sin(a)],[-np.sin(a), np.cos(a)]])
return np.dot(np.array(p), m)
>>> rotate((5.7,2.8), 35) # 旋转35°
array([3.06315263, 5.56301141])
>>> rotate((5.7,2.8), 90) # 旋转90°
array([-2.8, 5.7])
>>> rotate((5.7,2.8), 180) # 旋转180°
array([-5.7, -2.8])
>>> rotate((5.7,2.8), 360) # 旋转360°
array([5.7, 2.8])
随机抽样子模块
在科研活动和日常工作中，经常会用到随机数，特别是仿真、人工智能等领域更离不开随机数。计算机系统生成的随机数都是伪随机数，NumPy自然也不例外，但它的random随机抽样子模块提供了很多便捷的函数，可以满足绝大多数的应用需求。
随机数
np.random.random( )是最常用的随机数生成函数，该函数生成的随机数均匀分布于[0, 1)区间（左闭右开）。如果不提供参数，np.random.random( )函数返回一个浮点型随机数。np.random.random( )函数还可以接受一个整型或元组参数，用于指定返回的浮点型随机数数组的结构（shape）。也有很多人习惯使用np.random.rand( )函数生成随机数，其功能和np.random.random( )函数一样，只是np.random.rand( )函数不接受元组参数，必须要写成两个整型参数。
>>> np.random.random()
0.5325439713619637
>>> np.random.random(2)
array([0.58230722, 0.59662556])
>>> np.random.random((2,3))
array([[0.30103372, 0.16186198, 0.01281106],
[0.43598428, 0.58364464, 0.37358535]])
np.random.randint( )是另一个常用的随机数生成函数，该函数生成的随机整数均匀分布于[low, high)区间（左闭右开）。如果省略low参数，则默认low的值等于0。np.random.randint( )函数还有一个默认参数size，用于指定返回的整型随机数数组的结构（shape），其代码如下。
>>> np.random.randint(10)
3
>>> np.random.randint(10, size=5)
array([4, 5, 3, 0, 0])
>>> np.random.randint(10, size=(2,5))
array([[7, 0, 7, 8, 7],
[3, 9, 6, 9, 6]])
>>> np.random.randint(10, 100, size=(2,5))
array([[74, 84, 35, 78, 95],
[36, 28, 61, 25, 21]])
随机抽样
随机抽样是从指定的有序列表中随机抽取指定数量的元素。随机抽样的应用比较广泛，如产品抽检、抽签排序等。NumPy的随机抽样函数是np.random.choice( )，其原型如下。
np.random.choice(a, size=None, replace=True, p=None)
参数a表示待抽样的全体样本，它只接受整数或一维的数组（列表）。参数a如果是整数，相当于将数组np.arange(a)作为全体样本。参数size用于指定返回抽样结果数组的结构（shape）。参数replace用于指定是否允许多次抽取同一个样本，默认为允许。参数p是和全体样本集合等长的权重数组，用于指定对应样本被抽中的概率。
>>> np.random.choice(1,5) # 抽签样本只有1个元素0
array([0, 0, 0, 0, 0])
>>> np.random.choice(['a','b','c'], size=(3,5), p=[0.5,0.25,0.25]) # 指定权重
array([['a', 'a', 'a', 'c', 'b'],
['a', 'b', 'c', 'a', 'a'],
['a', 'a', 'a', 'a', 'b']], dtype='<U1')
>>> np.random.choice(np.arange(100), size=(2,5), replace=False) # 不允许重复
array([[39, 83, 57, 1, 45],
[ 6, 42, 81, 96, 80]])
正态分布
使用np.random.randn( )函数是最简单的生成标准正态分布随机数的方法。np.random.randn( )函数用于生成均值为0、标准差为1的正态分布（标准正态分布）随机数。该函数可以接受一个或两个整型参数，用来指定返回的符合标准正态分布的随机数数组的结构（shape）。
>>> np.random.randn()
-0.6810877665131689
>>> np.random.randn(5)
array([-0.66705291, -0.36616713, -1.16756153, 0.87119457, 0.30739537])
>>> np.random.randn(2,5)
array([[ 0.13724214, -0.71052821, 0.72988774, 0.16853684, -1.89836678],
[-0.94679367, -0.94165807, 2.31610439, -0.46910794, 0.02178768]])
如果需要生成非标准正态分布随机数，则应该使用np.random.normal( )函数。np.random.normal( )函数默认生成均值为0、标准差为1的正态分布随机数。参数loc用于指定均值，参数scale用于指定标准差，参数size用于指定返回的符合正态分布的随机数数组的结构（shape）。从下面的代码可以看出，和使用默认标准差相比，指定标准差为0.2时，数据分布更加靠近均值。
>>> np.random.normal()
0.49623559879844314
>>> np.random.normal(loc=2, size=5)
array([1.53367679, 1.64747689, 2.98585801, 0.41111733, 3.31062258])
>>> np.random.normal(loc=2, scale=0.2, size=(2,5))
array([[2.16220351, 1.73378492, 2.3342871 , 1.60388221, 1.85757531],
[2.05661231, 1.76707957, 1.90711622, 1.9695985 , 1.90088622]])
伪随机数的深度思考
计算机程序或编程语言中的随机数都是伪随机数。因为计算机硬件是确定的，代码是固定的，算法是准确的，通过这些确定的、固定的、准确的东西不会产生真正的随机数，除非引入这个封闭系统以外的因素。计算机系统的随机数算法一般使用线性同余或平方取中的算法，通过一个种子（通常用时钟代替）产生。这意味着，如果知道了种子和已经产生的随机数，就可能获得接下来随机数序列的信息，这就是伪随机数的可预测性。
NumPy随机数函数内部使用了一个伪随机数生成器，这个生成器每次实例化时都需要一个种子（整数）完成初始化。如果两次初始化的种子相同，则每次初始化后生成的随机数序列就完全一致。np.random.seed( )函数可以指定伪随机数生成器的初始化种子，其代码如下。
>>> np.random.seed(12345)
>>> np.random.random(5)
array([0.92961609, 0.31637555, 0.18391881, 0.20456028, 0.56772503])
>>> np.random.random((2,3))
array([[0.5955447 , 0.96451452, 0.6531771 ],
[0.74890664, 0.65356987, 0.74771481]])
>>> np.random.seed(12345)
>>> np.random.random(5) # 和上面完全一致
array([0.92961609, 0.31637555, 0.18391881, 0.20456028, 0.56772503])
>>> np.random.random((2,3)) # 和上面完全一致
array([[0.5955447 , 0.96451452, 0.6531771 ],
[0.74890664, 0.65356987, 0.74771481]])
从上述代码中可以看出，只要指定相同的种子，接下来的随机序列就完全一致。这意味着，只有从外部引入真正的随机因子（如天空云朵的形状、邻居家无线网络信号的强度等）作为种子，才可以得到真正的随机数。
此外，NumPy还提供了随机数生成器，可以直接操作这个生成器来生成随机数，其代码如下。
>>> r = np.random.RandomState(12345)
>>> r.random(5)
array([0.92961609, 0.31637555, 0.18391881, 0.20456028, 0.56772503])
>>> r.random((2,3))
array([[0.5955447 , 0.96451452, 0.6531771 ],
[0.74890664, 0.65356987, 0.74771481]])
