python科学计算-数值计算-NumPy

1、NumPy 的基本对象

NumPy和SciPy的所有科学计算功能建立在NumPy两种基本类型的对象之上。第一种对象是N维数组对象，即ndarray；第二种对象是通用函数对象，即ufunc。除了这两种对象，还有其他一些对象建立在它们之上。

1.1 N 维数组对象

ndarray对象是一个同质元素的集合，这些元素用N个整型数索引（N是数组的维数）。 ndarray有两个重要的属性。第一个属性是数组中元素的数据类型，称作dtype；第二个属性是数组的维度。数组的数据类型可以是Python支持的任意数据类型。数组的维度是一个元组（N-tuple），即一个包含N维数组的N个元素的集合，元组中的每个元素定义了数组在该维度包含的元素个数。

1. 数组的属性

除了维度和dtype，数组还有其他属性，列举如下：
 大小（size）
 每项的大小（itemsize）
 数据（data）
 维度（ndim）
每项的大小（itemsize）是数组中每个元素的字节长度。数据（data）属性是一个Python缓存对象，该对象指向数组数据的起始位置。可通过以下的Python程序来理解维度、数据类型和其他属性的概念：

import numpy as npx2d = np.array(((100, 200, 300),                (111, 222, 333),                (123, 456, 789)))print(x2d.shape)print(x2d.dtype)print(x2d.size)print(x2d.itemsize)print(x2d.ndim)print(x2d.data)

输出：

(3, 3)int32942<memory at 0x000000000331AEA0>

2.数组的基本操作

使用方括号（[]）来索引数组的值称作数组索引（array indexing）。以上面的程序中定义和用到的二维数组x2d为例，该数组中的特定元素可以用x2d[row, column]表示。例如，x2d[1,1]表示第二行的第二个元素（即222），索引值的起始值是0。同样， x2d[2,1]表示第三行的第二个元素（即456）。

数组切片是从数组中选定某些元素来构成一个子数组的过程。对于一维数组，我们可以依次从数组中选择某些元素。进而可以通过切片获取二维数组的一整行或一整列。也就是说我们可以在任一维度选取数组元素。 Python的基本切片概念被扩展到了N维数组中。切片的基本语法是start:stop:step。第一个参数指切片的起始位置索引值，第二个参数是指切片结束位置的索引值，最后一个参数定义了相加于前一个选定元素索引值基础上的步长。如果我们省略前两个参数值中的任意一个，那么这个值会被默认当作零或者大于零。同样，步长的默认值是1。下面用几个例子来更清楚地理解切片。

以x2d这个6×3的数组为例。 x2d[1]等价于x2d[1, :]，都表示数组的第二行，该行包含三个元素。另一方面， x2d[:, 1]表示数组的第二列，该列包含六个元素。 x2d[::3,1]可以选中数组中第二列的每三个元素。省略号可以用来替换零个或多个冒号。一个省略号可以当作零个或多个全切片的对象，以匹配切片对象的所有维度，这个维度必须和原始数组的维度相等。例如，如果x4d是一个5×6×7×8的数组，那么x4d[2 :, …, 6]等价于x4d[2 :, :, :, 6]。同样， x4d[…, 4]等价于x4d[:, :, :, 4]。可用以下程序理清数组切片的概念。这个程序展示了一维数组和二维数组的切片：

import numpy as npx = np.array([1, 12, 25, 8, 15, 35, 50, 7, 2, 10])x[3:7]x[1:9:2]x[0:9:3]x2d = np.array(((100, 200, 300),                (111, 222, 333),                (123, 456, 789),                (125, 457, 791),                (127, 459, 793),                (129, 461, 795)))print(x2d[0:4, 0:3])print(x2d[0:4:2, 0:3:2])

out：
[ 8 15 35 50]

[12 8 35 7]

[ 1 8 50]

[[100 200 300]
[111 222 333]
[123 456 789]
[125 457 791]]

[[100 300]
[123 789]]

3.数组的特殊操作（变形及转换）

为了实现数组的变形，有ravel、 reshape、 resize和指定数组维度属性等方法。 ravel和reshape方法返回修改了维度的变量（被调用的对象），而resize和指定数组维度属性的办法则修改了实际数组的维度。 ravel方法将数组展开成C语言类型的数组，它返回的变量被当作一个一维数组，这个数组按行依次展开成一个一维数组。

我们用下面的程序来讨论这些方法的影响。这个程序实现了二维数组的变形操作。第一个print函数输出的是原始数组。 ravel函数将展示展开的数组。而ravel函数之后的print函数再次展示了原始数组，因为ravel函数并未改变原始数组。现在uresize函数将原始的6行3列的数组(6,3)变形为3行6列的数组(3,6)。 因此resize之后的print函数将输出变形之后的数组形状。现在我们对原始数组(6,3)应用了reshape函数。但是由于reshape不改变原始数组的形状，在reshape函数之后的print函数将会打印出(3,6)的数组。最后一种方法是指定数组的形状为(9,2)，将数组的形状变成(9,2)。最重要的是要记住，当做reshape变换时，新数组的总大小仍然不变：

import numpy as npx2d = np.array(((100, 200, 300),                (111, 222, 333),                (123, 456, 789),                (125, 457, 791),                (127, 459, 793),                (129, 461, 795)))print(x2d)x2d.ravel()print(x2d)x2d.resize((3, 6))print(x2d)x2d.reshape(6, 3)print(x2d)x2d.shape = (9, 2)print(x2d)

输出：
[[100 200 300]
[111 222 333]
[123 456 789]
[125 457 791]
[127 459 793]
[129 461 795]]

[[100 200 300]
[111 222 333]
[123 456 789]
[125 457 791]
[127 459 793]
[129 461 795]]

[[100 200 300 111 222 333]
[123 456 789 125 457 791]
[127 459 793 129 461 795]]

[[100 200 300 111 222 333]
[123 456 789 125 457 791]
[127 459 793 129 461 795]]

[[100 200]
[300 111]
[222 333]
[123 456]
[789 125]
[457 791]
[127 459]
[793 129]
[461 795]]

4. 与数组相关的类

有一些类和子类是和ndarray类相关的。这些类的目的是为了支持特定的增强功能。下面将讨论这些类和子类。

(1) 矩阵子类

matrix类是ndarray的Python子类。一个矩阵能够通过其他矩阵或字符串生成，或者通过其
他可以转化为ndarray的对象生成。 matrix子类拥有特殊的覆盖运算符，例如是矩阵的乘， *
是矩阵的幂运算。 matrix类提供的一些函数可以实现多种功能，例如元素排序、转置计算、矩
阵元素的求和、将矩阵转换为列表或者其他数据结构和数据类型。下面的程序定义了两个3行3
列的矩阵。最后程序输出了两个矩阵的乘积：

import numpy as npa = np.matrix('1 2 3; 4 5 6; 7 8 9')print(a)b = np.matrix('4 5 6; 7 8 9; 10 11 12')print(b)print(a * b)

out:

[[1 2 3] [4 5 6] [7 8 9]][[ 4  5  6] [ 7  8  9] [10 11 12]][[ 48  54  60] [111 126 141] [174 198 222]]

(2) 掩码数组

NumPy有一个生成掩码数组的模块numpy.ma。掩码数组是一个包含一些非法的、缺失的、
未预料的实体的正常数组。掩码数组有两个组成成分：原始ndarray和一个掩码。掩码是由布尔逻
辑值组成的数组，可以用来判定数组的值是有效还是无效。掩码中的一个true值反应了数组里
相应的值是无效的。在掩码数组后面的计算中，这些无效的实体将不会用到。下一个程序展示了
掩码数组的概念。假设原始数组x里包含不同人的脉搏率，并且其中有两个非法实体。为了掩盖
这两个非法实体，掩码中相应的值被设置成1（true）。最后我们计算原始数组和掩码数组的平
均值。没用掩码处理的平均值是61.1，因为包含两个负值，而用掩码进行处理后的平均值是94.5：

import numpy as npimport numpy.ma as max = np.array([72, 79, 85, 90, 150, -135, 120, -10, 60, 100])mx = ma.masked_array(x, mask=[0, 0, 0, 0, 0, 1, 0, 1, 0, 0])mx2 = ma.masked_array(x, mask=x < 0)print(x.mean())print(mx.mean())print(mx2.mean())

out：

61.194.594.5

注：mean()方法返回的数组元素的平均值

(3) 结构化的数组

NumPy的ndarray能包含记录类型的值。为了生成一个记录类型的数组，首先需要生成记录类型的数据类型，然后用这种数据类型的元素构建数组。这种数据类型可以通过dtype在数组定义中定义。下面的程序中生成了一个记录类型的数组，这个数组中包含城市的最低、最高和平均气温。 dtype函数由两部分组成：域的命名和格式。例子中用到的格式是32位整型（i4）、 32位浮点型（f4）和包含30或更少的字符的字符串（a30）：

import numpy as nprectype = np.dtype({'names': ['mintemp', 'maxtemp', 'avgtemp', 'city'],                    'formats': ['i4', 'i4', 'f4', 'a30']})a = np.array([(10, 44, 25.2, 'Indore'), (10, 42, 25.2, 'Mumbai'), (2, 48, 30,                                                                   'Delhi')], dtype=rectype)print(a[0])print(a['mintemp'])print(a['maxtemp'])print(a['avgtemp'])print(a['city'])

out:

(10, 44, 25.200000762939453, b'Indore')[10 10  2][44 42 48][ 25.20000076  25.20000076  30.        ][b'Indore' b'Mumbai' b'Delhi']

1.2 通用函数对象

通用函数（unfunc）是对ndarray进行一个元素一个元素操作的函数。它也支持广播、类型
转换和一些其他重要的特征。在NumPy中，广播是对不同维度数组的操作过程。特别是在数学运
算中，较小维度的数组要在大维度的数组中广播，以使其维度和大维度数组兼容。通用函数即是
NumPy中ufunc类的实例。

1. 属性

每个通用函数都有一些属性，但是用户无法设置这些属性的值。下面是通用函数的属性，这
是一些通用函数的信息属性。
 doc：它包含了ufunc函数的doc字符串。其第一部分是基于命名、输入的数量和输
出的数量动态生成的。第二部分是在函数生成时定义的，并存在函数中。
 name：这是ufunc的命名。
 ufunc.nin：这表示总的输入数量。
 ufunc.nout：这表示总的输出数量。
 ufunc.nargs：这表示总的变量数，包括输入和输出的数量。
 ufunc.ntypes：这表示这个函数定义的类型的种类。
 ufunc.types：这个属性返回函数定义的ntypes种具体类型的列表。
 ufunc.identity：这表示这个函数的值。

2. 方法

所有的通用函数ufunc有5种方法，如下所示。前4种方法的通用函数包含两个输入变量，返
回一个输出变量。这些方法在调用其他函数失败时会报出ValueError异常。第五种方法允许用
户利用索引做位置操作。下面是每个NumPy通用函数可用的方法。
 ufunc.reduce：在一个坐标轴应用通用函数时，它可以将数组降低一个维度。
 ufunc.accumulate：它可以在对所有元素使用同一个运算符时积累结果。
 ufunc.reduceat：它可以在单个坐标轴上减少指定的切片。
 ufunc.outer(A, B)：它对所有(a, b)应用通用函数操作符，其中a包含于A， b包含于B。
 ufunc.at：它对特定元素执行未缓存的位置操作。

3. 各种可用的通用函数

目前NumPy支持超过60种的通用函数。这些函数包括广泛的操作，如简单的数学运算（加、
减、求模、取绝对值）、开方、幂运算、指数运算、三角运算、比特位运算、比较和浮点运算。
通常最好选择利用这些函数而不选择循环，因为相比于循环这些函数更高效。
下面的程序展示了一些通用函数的运用：

import numpy as npx1 = np.array([72, 79, 85, 90, 150, -135, 120, -10, 60, 100])x2 = np.array([72, 79, 85, 90, 150, -135, 120, -10, 60, 100])x_angle = np.array([30, 60, 90, 120, 150, 180])x_sqr = np.array([9, 16, 25, 225, 400, 625])x_bit = np.array([2, 4, 8, 16, 32, 64])print(np.greater_equal(x1, x2))print(np.mod(x1, x2))print(np.exp(x1))print(np.reciprocal(x1))print(np.negative(x1))print(np.isreal(x1))# print(np.isnan(np.log10(x1)))print(np.sqrt(np.square(x_sqr)))print(np.sin(x_angle * np.pi / 180))print(np.tan(x_angle * np.pi / 180))print(np.right_shift(x_bit, 1))print(np.left_shift(x_bit, 1))

[ True  True  True  True  True  True  True  True  True  True][0 0 0 0 0 0 0 0 0 0][  1.85867175e+31   2.03828107e+34   8.22301271e+36   1.22040329e+39   1.39370958e+65   2.34555134e-59   1.30418088e+52   4.53999298e-05   1.14200739e+26   2.68811714e+43][0 0 0 0 0 0 0 0 0 0][ -72  -79  -85  -90 -150  135 -120   10  -60 -100][ True  True  True  True  True  True  True  True  True  True][   9.   16.   25.  225.  400.  625.][  5.00000000e-01   8.66025404e-01   1.00000000e+00   8.66025404e-01   5.00000000e-01   1.22464680e-16][  5.77350269e-01   1.73205081e+00   1.63312394e+16  -1.73205081e+00  -5.77350269e-01  -1.22464680e-16][ 1  2  4  8 16 32][  4   8  16  32  64 128]

在Python中，如果有一个值不能用数字表示，那么这个值就是空值nan。例如，如果对前一个程序中的x1数组做log10的通用函数运算，那么输出值就是nan。有一个通用函数isnan可用来检验一个值是不是nan。三角函数需要角度值度数的变量。十进制常数值是弧度值，可以通过乘以180/numpy.pi来实现度数的转换。比特位左移1位执行的是变量乘以2。类似地，比特位右移1位执行的是变量除以2。通常，这些通用函数都是对数组进行操作，如果有非数组的变量的话，那么这个变量通过广播机制当作一个数组，然后执行逐元素的操作。这就是前一个程序中最后四行的操作。

1.3 NumPy 的数学模块

NumPy加入了特定功能的模块，例如线性代数、离散傅里叶变换、随机采样和矩阵代数库。这些功能捆绑在单独的模块中，这些模块列举如下。

 numpy.linalg：这个模块支持线性代数的各种功能，如数组和向量的内积、外积和点积；向量和矩阵的范数；线性矩阵方程的解；矩阵转置的方法。
 numpy.fft：离散傅里叶变换在数字信号处理中有广泛的应用。这个模块中的函数可以计算各种类型的离散傅里叶变换，包括一维、二维、多维、转置和傅里叶变换。
 numpy.matlib：这个模块包括那些默认返回矩阵对象而不是ndarray对象的函数。这些函数包括empty、 zeros、 ones、 eye、 rapmat、 rand、 randn、 bmat、 mat和matrix。
 numpy.random：这个模块包括在特定人群或范围中执行随机抽样的函数。它也支持随机排列组合的生成。另外，它还包括一些支持各种基于统计分布生成的随机抽样数值的函数。

下面一个程序展示了linalg模块中的一些函数的应用。它计算了范数、转置、行列式、特征值和方阵右特征向量。它还通过求解两个方程2x+3y=4和3x+4y=5展示了线性方程的求解，并且是通过将它们当作数组来求解的。最后一行的allclose函数比较了传入的两个数组，并且当它们的每个元素都相等时返回true。 eig方法计算了方阵的特征值和特征向量。返回值如下： w是特征值， v是特征向量， v[:,i]列是w[i]的特征向量。

import numpy as npfrom numpy import linalg as LAarr2d = np.array(((100, 200, 300),                  (111, 222, 333),                  (129, 461, 795)))eig_val, eig_vec = LA.eig(arr2d)LA.norm(arr2d)LA.det(arr2d)LA.inv(arr2d)arr1 = np.array([[2, 3], [3, 4]])arr2 = np.array([4, 5])results = np.linalg.solve(arr1, arr2)print(results)print(np.allclose(np.dot(arr1, results), arr2))

out:

[-1.  2.]True

随机抽样是科学和商业计算中很重要的方面。下面的程序展示了numpy随机抽样模块支持的各类函数中的部分函数。除了样本维度与总体，有些分布还需要一些统计值，例如均值、众数、标准差。 permutation函数随机排列一个序列或者返回一个排列范围， randint函数随机返回最初的两个变量所指定的范围中的一些元素，而元素的总个数由第三个变量指定。剩余的方法返回特定分布中的抽样，如卡方检验、帕雷托分布、正态分布和对数正态分布。

import numpy as npnp.random.permutation(10)np.random.randint(20, 50, size=10)np.random.random_sample(10)np.random.chisquare(5, 10)  # 自由度alpha, location_param = 4., 2.s = np.random.pareto(alpha, 10) + location_params = np.random.standard_normal(20)mean, std_deviation = 4., 2.s = np.random.lognormal(mean, std_deviation, 10)