Python 数据分析——NumPy 函数库

除了前面介绍的ndarray数组对象和ufunc函数之外，NumPy还提供了大量对数组进行处理的函数。充分利用这些函数，能够简化程序的逻辑，提高运算速度。我们通过一些较常用的例子，说明它们的一些使用技巧和注意事项。

一、随机数

本节介绍的函数如表1所示。

表1 本节要介绍的函数

函数名	功能	函数名	功能
rand	0到1之间的随机数	randn	标准正态分布的随机数
randint	指定范围内的随机整数	normal	正态分布
uniform	均匀分布	poisson	泊松分布
permutation	随机排列	shuffle	随机打乱顺序
choice	随机抽取样本	seed	设置随机数种子

numpy.random模块中提供了大量的随机数相关的函数，为了方便后面用随机数测试各种运算函数，让我们首先来看看如何产生随机数：

• rand()产生0到1之间的随机浮点数，它的所有参数用于指定所产生的数组的形状。
• randn()产生标准正态分布的随机数，参数的含义与rand()相同。
• randint()产生指定范围的随机整数，包括起始值，但是不包括终值，在下面的例子中，产生0到9的随机数，它的第三个参数用于指定数组的形状：

random模块提供了许多产生符合特定随机分布的随机数的函数，它们的最后一个参数size都用于指定输出数组的形状，而其他参数都是分布函数的参数。例如：

• normal()：正态分布，前两个参数分别为期望值和标准差。
• uniform()：均匀分布，前两个参数分别为区间的起始值和终值。
• poisson()：泊松分布，第一个参数指定λ系数，它表示单位时间(或单位面积)内随机事件的平均发生率。由于泊松分布是一个离散分布，因此它输出的数组是一个整数数组。

permutation()可以用于产生一个乱序数组，当参数为整数n时，它返回[0,n)这n个整数的随机排列；当参数为一个序列时，它返回一个随机排列之后的序列：

a = np.array([1, 10, 20, 30, 40])print nr.permutation(10)print nr.permutation(a)[2 4 3 5 6 8 0 1 9 7][40 1 10 20 30]

permutation()返回一个新数组，而shuffle()则直接将参数数组的顺序打乱：

nr.shuffle(a)aarray([ 1, 20, 30, 10, 40])

choice()从指定的样本中随机进行抽取：

• size参数用于指定输出数组的形状。
• replace参数为True时，进行可重复抽取，而为False时进行不重复抽取，默认值为True。所以在下面的例子中，c1中可能有重复数值，而c2中的每个数值都是不同的。
• p参数指定每个元素对应的抽取概率，如果不指定，所有的元素被抽取到的概率相同。在下面的例子中，值越大的元素被抽到的概率越大，因此c3中数值较大的元素比较多。

为了保证每次运行时能重现相同的随机数，可以通过seed()函数指定随机数的种子。在下面的例子中，计算r3和r4之前，都使用42作为种子，因此得到的随机数组是相同的：

二、求和、平均值、方差

本节介绍的函数如表2所示。

表2　本节要介绍的函数

函数名	功能	函数名	功能
sum	求和	mean	求期望
average	加权平均数	std	标准差
var	方差	product	连乘积

sum()计算数组元素之和，也可以对列表、元组等与数组类似的序列进行求和。当数组是多维时，它计算数组中所有元素的和。这里我们使用random.randint()模块中的函数创建一个随机整数数组。

np.random.seed(42)a = np.random.randint(0,10,size=(4,5))a np.sum(a)----------------- ---------[[6, 3, 7, 4, 6], 96[9, 2, 6, 7, 4],[3, 7, 7, 2, 5],[4, 1, 7, 5, 1]]

如果指定axis参数，则求和运算沿着指定的轴进行。在上面的例子中，数组a的第0轴的长度为4，第1轴的长度为5。如果axis参数为1，则对每行上的5个数求和，所得的结果是长度为4的一维数组。如果参数axis为0，则对每列上的4个数求和，结果是长度为5的一维数组。即结果数组的形状是原始数组的形状除去其第axis个元素：

np.sum(a, axis=1) np.sum(a, axis=0)----------------- --------------------[26, 28, 24, 18] [22, 13, 27, 18, 16]

当axis参数是一个轴的序列时，对指定的所有轴进行求和运算。例如下面的程序对一个形状为(2,3,4)的三维数组的第0和第2轴求和，得到的结果为一个形状为(3,)的数组。由于数组的所有元素都为1，因此求和的结果都是8：

np.sum(np.ones((2, 3, 4)), axis=(0, 2))array([ 8., 8., 8.])

有时我们希望能够保持原数组的维数，这时可以设置keepdims参数为True：

np.sum(a, 1, keepdims=True) np.sum(a, 0, keepdims=True)--------------------------- ---------------------------[[26], [[22, 13, 27, 18, 16]][28],[24],[18]]

sum()默认使用和数组的元素类型相同的累加变量进行计算，如果元素类型为整数，则使用系统的默认整数类型作为累加变量，例如在32位系统中使用32位整数作为累加变量。因此对整数数组进行累加时可能会出现溢出问题，即数组元素的总和超过了累加变量的取值范围。下面的程序计算数组a中每个元素占其所在行总和的百分比。在调用sum()函数时：

• 设置dtype参数为float，这样得到的结果是浮点数组，能避免整数的整除运算。
• 设置keepdims参数为True，这样sum()得到的结果的形状为(4,1)，能够和原始数组进行广播运算。

pa = a / np.sum(a, 1, dtype=float, keepdims=True) * 100pa pa.sum(1, keepdims=True)------------------------------------------ ------------------------[[ 23.08, 11.54, 26.92, 15.38, 23.08],[[ 100.],[ 32.14, 7.14, 21.43, 25. , 14.29],[ 100.],[ 12.5 , 29.17, 29.17, 8.33, 20.83],[ 100.],[ 22.22, 5.56, 38.89, 27.78, 5.56]][ 100.]]

对很大的单精度浮点数类型的数组进行计算时，也可能出现精度不够的现象，这时也可以通过dtype参数指定累加变量的类型。在下面的例子中，我们对一个元素都为1.1的单精度数组进行求和，比较单精度累加变量和双精度累加变量的计算结果：

np.set_printoptions(precision=8)b = np.full(1000000, 1.1, dtype=np.float32) # 创建一个很大的单精度浮点数数组b # 1.1无法使用浮点数精确表示，存在一些误差array([ 1.10000002, 1.10000002, 1.10000002, ..., 1.10000002,1.10000002, 1.10000002], dtype=float32)

使用单精度累加变量进行累加计算，误差将越来越大，而使用双精度浮点数则能够得到较精确的结果：

np.sum(b) np.sum(b, dtype=np.double)--------- --------------------------1099999.3 1100000.0238418579

上面的例子将产生一个新的数组来保存求和的结果，如果希望将结果直接保存到另一个数组中，可以和ufunc函数一样使用out参数指定输出数组，它的形状必须和结果数组的形状相同。

mean()求数组的平均值，它的参数与sum()相同。和sum()不同的是：对于整数数组它使用双精度浮点数进行计算，而对于其他类型的数组，则使用和数组元素类型相同的累加变量进行计算：

np.mean(a, axis=1) # 整数数组使用双精度浮点数进行计算array([ 5.2, 5.6, 4.8, 3.6])np.mean(b) np.mean(b, dtype=np.double)---------- ---------------------------1.0999993 1.1000000238418579

此外，average()也可以对数组进行平均计算。它没有out和dtype参数，但有一个指定每个元素权值的weights参数，可以用于计算加权平均数。例如有三个班级，number数组中保存每个班级的人数，score数组中保存每个班级的平均分，下面计算所有班级的加权平均分，得到整个年级的平均分：

score = np.array([83, 72, 79])number = np.array([20, 15, 30])print np.average(score, weights=number)78.6153846154

相当于进行如下计算：

print np.sum(score * number) / np.sum(number, dtype=float)78.6153846154

std()和var()分别计算数组的标准差和方差，有axis、out、dtype以及keepdims等参数。方差有两种定义：偏样本方差(biased sample variance)和无偏样本方差(unbiased sample variance)。偏样本方差的计算公式为：

而无偏样本方差的公式为：

当ddof参数为0时，计算偏样本方差；当ddof为1时，计算无偏样本方差，默认值为0。下面我们用程序演示这两种方差的差别。

首先产生一个标准差为2.0、方差为4.0的正态分布的随机数组。我们可以认为总体样本的方差为4.0。假设从总体样本中随机抽取10个样本，我们分别计算这10个样本的两种方差，这里我们用一个二维数组重复上述操作100000次，然后计算所有这些方差的期望值：

a = nr.normal(0, 2.0, (100000, 10))v1 = np.var(a, axis=1, ddof=0) #可以省略ddof=0v2 = np.var(a, axis=1, ddof=1)np.mean(v1) np.mean(v2)------------------ ------------------3.6008566906846693 4.0009518785385216

可以看到无偏样本方差的期望值接近于总体方差4.0，而偏样本方差比4.0小一些。

偏样本方差是正态分布随机变量的最大似然估计。如果有一个样本包含n个随机数，并且知道它们符合正态分布，通过该样本可以估算出正态分布的概率密度函数的参数。所估算的那组正态分布参数最符合给定的样本，就称为最大似然估计。

正态分布的概率密度函数的定义如下，其中μ表示期望，σ2表示方差：

所谓最大似然估计，就是找到一组参数使得下面的乘积最大，其中x1为样本中的值：

f(x1)f(x2)…f(xn)

专业术语总是很难理解，下面我们还是用程序来验证：

def normal_pdf(mean, var, x):return 1 / np.sqrt(2 * np.pi * var) * np.exp(-(x - mean) ** 2 / (2 * var))nr.seed(42)data = nr.normal(0, 2.0, size=10) ❶mean, var = np.mean(data), np.var(data) ❷var_range = np.linspace(max(var - 4, 0.1), var + 4, 100) ❸p = normal_pdf(mean, var_range[:, None], data) ❹p = np.product(p, axis=1) ❺

normal_pdf()为计算正态分布的概率密度的函数。❶产生10个正态分布的随机数。❷计算其最大似然估计的参数。❸以最大似然估计的方差为中心，产生一组方差值。❹用正态分布的概率密度函数计算每个样本、每个方差所对应的概率密度。由于使用了广播运算，得到的结果p是一个二维数组，它的第0轴对应var_range中的各个方差，第1轴对应data中的每个元素。❺沿着p的第1轴求所有概率密度的乘积。product()和sum()的用法类似，用于计算数组所有元素的乘积。

下面绘制var_range中各个方差对应的似然估计值，并用一条竖线表示偏样本方差。由图1可以看到偏样本方差位于似然估计曲线的最大值处。

import pylab as plpl.plot(var_range, p)pl.axvline(var, 0, 1, c='r')pl.show()

图1　偏样本方差位于似然估计曲线的最大值处

三、大小与排序

本节介绍的函数如表3所示。

表3 本节要介绍的函数

函数名	功能	函数名	功能
min	最小值	max	最大值
minimum	二元最小值	maximum	二元最大值
ptp	最大值与最小值的差	argmin	最小值的下标
argmax	最大值的下标	unravel_index	一维下标转换成多维下标
sort	数组排序	argsort	计算数组排序的下标
lexsort	多列排序	partition	快速计算前k位
argpartition	前k位的下标	median	中位数
percentile	百分位数	searchsorted	二分查找

用min()和max()可以计算数组的最小值和最大值，它们都有axis、out、keepdims等参数。这些参数的用法和sum()基本相同，但是axis参数不支持序列。此外，ptp()计算最大值和最小值之间的差，有axis和out参数。这里就不再多举例了，请读者自行查看函数的文档。minimum()和maximum()用于比较两个数组对应下标的元素，返回数组的形状为两参数数组广播之后的形状。

a = np.array([1, 3, 5, 7])b = np.array([2, 4, 6])np.maximum(a[None, :], b[:, None])array([[2, 3, 5, 7],[4, 4, 5, 7],[6, 6, 6, 7]])

用argmax()和argmin()可以求最大值和最小值的下标。如果不指定axis参数，则返回平坦化之后的数组下标，例如下面的程序找到a中最大值的下标，有多个最值时得到第一个最值的下标：

np.random.seed(42)a = np.random.randint(0, 10, size=(4, 5))max_pos = np.argmax(a)max_pos5

下面查看a平坦化之后的数组中下标为max_pos的元素：

a.ravel()[max_pos] np.max(a)------------------ ---------9 9

可以通过unravel_index()将一维数组的下标转换为多维数组的下标，它的第一个参数为一维数组的下标，第二个参数是多维数组的形状：

idx = np.unravel_index(max_pos, a.shape)idx a[idx]------ ------(1, 0) 9

当使用axis参数时，可以沿着指定轴计算最大值的下标。例如下面的结果表示，在数组a中第0行的最大值的下标为2，第1行的最大值的下标为0：

idx = np.argmax(a, axis=1)idxarray([2, 0, 1, 2])

使用下面的语句可以通过idx选出每行的最大值：

a[np.arange(a.shape[0]), idx]array([7, 9, 7, 7])

数组的sort()方法对数组进行排序，它会改变数组的内容；而sort()函数则返回一个新数组，不改变原始数组。它们的axis默认值都为-1，即沿着数组的最终轴进行排序。sort()函数的axis参数可以设置为None，此时它将得到平坦化之后进行排序的新数组。在下面的例子中，np.sort(a)对a中每行的数值进行排序，而np.sort(a, axis=0)对数组a每列上的数值进行排序。

np.sort(a) np.sort(a, axis=0)----------------- ------------------[[3, 4, 6, 6, 7], [[3, 1, 6, 2, 1],[2, 4, 6, 7, 9], [4, 2, 7, 4, 4],[2, 3, 5, 7, 7], [6, 3, 7, 5, 5],[1, 1, 4, 5, 7]] [9, 7, 7, 7, 6]]

argsort()返回数组的排序下标，参数axis的默认值为-1：

sort_axis1 = np.argsort(a)sort_axis0 = np.argsort(a, axis=0)sort_axis1 sort_axis0----------------- -----------------[[1, 3, 0, 4, 2], [[2, 3, 1, 2, 3],[1, 4, 2, 3, 0], [3, 1, 0, 0, 1],[3, 0, 4, 1, 2], [0, 0, 2, 3, 2],[1, 4, 0, 3, 2]] [1, 2, 3, 1, 0]]

为了使用sort_axis0和sort_axis1计算排序之后的数组，即np.sort(a)的结果，需要产生非排序轴的下标。下面使用ogrid对象产生第0轴和第1轴的下标axis0和axis1：

axis0, axis1 = np.ogrid[:a.shape[0], :a.shape[1]]

然后使用这些下标数组得到排序之后的数组：

a[axis0, sort_axis1] a[sort_axis0, axis1]-------------------- --------------------[[3, 4, 6, 6, 7], [[3, 1, 6, 2, 1],[2, 4, 6, 7, 9], [4, 2, 7, 4, 4],[2, 3, 5, 7, 7], [6, 3, 7, 5, 5],[1, 1, 4, 5, 7]] [9, 7, 7, 7, 6]]

使用这种方法可以对两个相关联的数组进行排序，即从数组a产生排序下标数组，然后使用它对数组b进行排序。排序相关的函数或方法还可以通过kind参数指定排序算法，对于结构数组可以通过order参数指定排序所使用的字段。

lexsort()类似于Excel中的多列排序。它的参数是一个形状为(k, N)的数组，或者包含k个长度为N的数组的序列，可以把它理解为Excel中N行k列的表格。lexsort()返回排序下标，注意数组中最后的列为排序的主键。在下面的例子中，按照“姓名-年龄”的顺序对数据排序：

names = ['zhang', 'wang', 'li', 'wang', 'zhang']ages = [37, 33, 32, 31, 36]idx = np.lexsort([ages, names])sorted_data = np.array(zip(names, ages), 'O')[idx]idx sorted_data--------------- ---------------[2, 3, 1, 4, 0][['li', 32],['wang', 31],['wang', 33],['zhang', 36],['zhang', 37]]

如果需要对一个N行k列的数组以第一列为主键进行排序，可以先通过切片下标::-1反转数组的第1轴，然后对其转置进行lexsort()排序：

b = np.random.randint(0, 10, (5, 3))b b[np.lexsort(b[:, ::-1].T)]----------- ---------------------------[[4, 0, 9], [[3, 8, 2],[5, 8, 0], [4, 0, 9],[9, 2, 6], [4, 2, 6],[3, 8, 2], [5, 8, 0],[4, 2, 6]] [9, 2, 6]]

partition()和argpartition()对数组进行分割，可以很快地找出排序之后的前k个元素，由于它不需要对整个数组进行完整排序，因此速度比调用sort()之后再取前k个元素要快许多。下面从10万个随机数中找出前5个最小的数，注意partition()得到的前5个数值没有按照从小到大排序，如果需要，可以再调用sort()对这5个数进行排序即可：

r = np.random.randint(10, 1000000, 100000)np.sort(r)[:5] np.partition(r, 5)[:5]-------------------- ----------------------[15, 23, 25, 37, 47] [15, 47, 25, 37, 23]

下面用%timeit测试sort()和partition()的运行速度：

%timeit np.sort(r)[:5]%timeit np.sort(np.partition(r, 5)[:5])100 loops, best of 3: 6.02 ms per loop1000 loops, best of 3: 348 µs per loop

用median()可以获得数组的中值，即对数组进行排序之后，位于数组中间位置的值。当长度是偶数时，则得到正中间两个数的平均值。它也可以指定axis和out参数：

np.median(a, axis=1)array([ 6., 6., 5., 4.])

percentile()用于计算百分位数，即将数值从小到大排列，计算处于p%位置上的值。下面的程序计算标准正态分布随机数的绝对值在68.3%、95.4%以及99.7%处的百分位数，它们应该约等于1倍、2倍和3倍的标准差：

r = np.abs(np.random.randn(100000))np.percentile(r, [68.3, 95.4, 99.7])array([ 1.00029686, 1.99473003, 2.9614485 ])

当数组中的元素按照从小到大的顺序排列时，可以使用searchsorted()在数组中进行二分搜索。在下面的例子中，a是一个已经排好序的列表，v是需要搜索的数值列表。searchsorted()返回一个下标数组，将v中对应的元素插入到a中的位置，能够保持数据的升序排列。当v中的元素在a中出现时，通过side参数指定返回最左端的下标还是最右端的下标。在下面的例子中，16放到a的下标为3、4、5的位置都能保持升序排列，side参数为默认值'left'时返回3，而为'right'时返回5。

a = [2, 4, 8, 16, 16, 32]v = [1, 5, 33, 16]np.searchsorted(a, v) np.searchsorted(a, v, side='right')--------------------- -----------------------------------[0, 2, 6, 3] [0, 2, 6, 5]

searchsorted()可以用于在两个数组中查找相同的元素。下面看一个比较复杂的例子：有x和y两个一维数组，找到y中每个元素在x中的下标。若不存在，将下标设置为-1。

x = np.array([3, 5, 7, 1, 9, 8, 6, 10])y = np.array([2, 1, 5, 10, 100, 6])def get_index_searchsorted(x, y):index = np.argsort(x) ❶sorted_x = x[index] ❷sorted_index = np.searchsorted(sorted_x, y) ❸yindex = np.take(index, sorted_index, mode='clip') ❹mask = x[yindex] != y ❺yindex[mask] = -1return yindexget_index_searchsorted(x, y)array([-1, 3, 1, 7, -1, 6])

❶由于x并不是按照升序排列，因此先调用argsort()获得升序排序的下标index。❷使用index获得将x排序之后的sorted_x。❸使用searchsorted()在sorte_x中搜索y中每个元素对应的下标sorted_index。

❹如果搜索的值大于x的最大值，那么下标会越界，因此这里调用take()函数，take(index, sorted_index)与index[sorted_index]的含义相同，但是能处理下标越界的情况。通过设置mode参数为'clip'，将下标限定在0到len(x)-1之间。

❺使用yindex获取x中的元素并和y比较，若值相同则表示该元素确实存在于x之中，否则表示不存在。

这段算法有些复杂，但由于利用了NumPy提供的数组操作函数，它的运算速度比使用字典的纯Python程序要快。下面我们用两个较大的数组测试运算速度。为了比较的公平性，我们调用tolist()方法将数组转换成列表：

x = np.random.permutation(1000)[:100]y = np.random.randint(0, 1000, 2000)xl, yl = x.tolist(), y.tolist()def get_index_dict(x, y):idx_map = {v:i for i,v in enumerate(x)}yindex = [idx_map.get(v, -1) for v in y]return yindexyindex1 = get_index_searchsorted(x, y)yindex2 = get_index_dict(xl, yl)print np.all(yindex1 == yindex2)%timeit get_index_searchsorted(x, y)%timeit get_index_dict(xl, yl)True10000 loops, best of 3: 122 µs per loop1000 loops, best of 3: 368 µs per loop

四、统计函数

本节介绍的函数如表4所示。

表4 本节要介绍的函数

函数名	功能	函数名	功能
unique	去除重复元素	bincount	对整数数组的元素计数
histogram	一维直方图统计	digitze	离散化

unique()返回其参数数组中所有不同的值，并且按照从小到大的顺序排列。它有两个可选参数：

• return_index：Ture表示同时返回原始数组中的下标。
• return_inverse：True表示返回重建原始数组用的下标数组。

下面通过几个例子介绍unique()的用法。首先用randint()创建有10个元素、值在0到9范围之内的随机整数数组，通过unique(a)可以找到数组a中所有的整数，并按照升序排列：

np.random.seed(42)a = np.random.randint(0, 8, 10)a np.unique(a)------------------------------ ------------------[6, 3, 4, 6, 2, 7, 4, 4, 6, 1] [1, 2, 3, 4, 6, 7]

如果参数return_index为True，则返回两个数组，第二个数组是第一个数组在原始数组中的下标。在下面的例子中，数组index保存的是数组x中每个元素在数组a中的下标：

x, index = np.unique(a, return_index=True)x index a[index]------------------ ------------------ ------------------[1, 2, 3, 4, 6, 7] [9, 4, 1, 2, 0, 5] [1, 2, 3, 4, 6, 7]

如果参数return_inverse为True，则返回的第二个数组是原始数组a的每个元素在数组x中的下标：

x, rindex = np.unique(a, return_inverse=True)rindex x[rindex]------------------------------ ------------------------------[4, 2, 3, 4, 1, 5, 3, 3, 4, 0] [6, 3, 4, 6, 2, 7, 4, 4, 6, 1]

bincount()对整数数组中各个元素所出现的次数进行统计，它要求数组中的所有元素都是非负的。其返回数组中第i个元素的值表示整数i出现的次数。

np.bincount(a)array([0, 1, 1, 1, 3, 0, 3, 1])

由上面的结果可知，在数组a中有1个1、1个2、1个3、3个4、3个6和1个7，而0、5等数没有在数组a中出现。

通过weights参数可以指定每个数所对应的权值。当指定weights参数时，bincount(x, weights=w)返回数组x中的每个整数所对应的w中的权值之和。用文字解释比较难以理解，下面我们看一个实例：

x = np.array([0 , 1, 2, 2, 1, 1, 0])w = np.array([0.1, 0.3, 0.2, 0.4, 0.5, 0.8, 1.2])np.bincount(x, w)array([ 1.3, 1.6, 0.6])

在上面的结果中，1.3是数组x中0所对应的w中的元素(0.1和1.2)之和，1.6是1所对应的w中的元素(0.3、0.5和0.8)之和，而0.6是2所对应的w中的元素(0.2和0.4)之和。如果要求平均值，可以用求和的结果与次数相除：

np.bincount(x, w) / np.bincount(x)array([ 0.65 , 0.53333333, 0.3 ])

histogram()对一维数组进行直方图统计。其参数列表如下：

histogram(a, bins=10, range=None, weights=None, density=False)

其中a是保存待统计数据的数组，bins指定统计的区间个数，即对统计范围的等分数。range是一个长度为2的元组，表示统计范围的最小值和最大值，默认值为None，表示范围由数据的范围决定，即(a.min(), a.max())。当density参数为False时，函数返回a中的数据在每个区间的个数，参数为True则返回每个区间的概率密度。weights参数和bincount()的类似。

histogram()返回两个一维数组—— hist和bin_edges，第一个数组是每个区间的统计结果，第二个数组的长度为len(hist) + 1，每两个相邻的数值构成一个统计区间。下面我们看一个例子：

a = np.random.rand(100)np.histogram(a, bins=5, range=(0, 1))(array([28, 18, 17, 19, 18]), array([ 0. , 0.2, 0.4, 0.6, 0.8, 1. ]))

首先创建了一个有100个元素的一维随机数组a，取值范围在0到1之间。然后用histogram()对数组a中的数据进行直方图统计。结果显示有28个元素的值在0到0.2之间，18个元素的值在0.2到0.4之间。读者可以尝试用rand()创建更大的随机数组，由统计结果可知每个区间出现的次数近似相等，因此rand()所创建的随机数在0到1范围之间是平均分布的。

如果需要统计的区间的长度不等，可以将表示区间分隔位置的数组传递给bins参数，例如：

np.histogram(a, bins=[0, 0.4, 0.8, 1.0])(array([46, 36, 18]), array([ 0. , 0.4, 0.8, 1. ]))

结果表示0到0.4之间有46个值，0.4到0.8之间有36个值。

如果用weights参数指定了数组a中每个元素所对应的权值，则histogram( )对区间中数值所对应的权值进行求和。下面看一个使用histogram( )统计男性青少年年龄和身高的例子。“height.csv”文件是100名年龄在7到20岁之间的男性青少年的身高统计数据。

首先用loadtxt( )从数据文件载入数据。在数组d中，第0列是年龄，第1列是身高。可以看到年龄的范围在7到20之间：

下面对数据进行统计，sums是每个年龄段的身高总和，cnts是每个年龄段的数据个数，因此很容易计算出每个年龄段的平均身高：

五、分段函数

本节介绍的函数如表5所示。

表5　本节要介绍的函数

函数名	功能
where	矢量化判断表达式
piecewise	分段函数
select	多分支判断选择

在前面的小节中介绍过如何使用frompyfunc( )函数计算三角波形。由于三角波形是分段函数，需要根据自变量的取值范围决定计算函数值的公式，因此无法直接通过ufunc函数计算。NumPy提供了一些计算分段函数的方法。

在Python 2.6中新增加了如下判断表达式语法，当condition条件为True时，表达式的值为y，否则为z：

x = y if condition else z

在NumPy中，where( )函数可以看作判断表达式的数组版本：

x = where(condition, y, z)

其中condition、y和z都是数组，它的返回值是一个形状与condition相同的数组。当condition中的某个元素为True时，x中对应下标的值从数组y获取，否则从数组z获取：

x = np.arange(10)np.where(x < 5, 9 - x, x)array([9, 8, 7, 6, 5, 5, 6, 7, 8, 9])

如果y和z是单个数值或者它们的形状与condition的不同，将先通过广播运算使其形状一致：

np.where(x > 6, 2 * x, 0)array([ 0, 0, 0, 0, 0, 0, 0, 14, 16, 18])

使用where( )很容易计算前面介绍过的三角波形：

def triangle_wave1(x, c, c0, hc):x = x - x.astype(np.int) # 三角波的周期为1，因此只取x坐标的小数部分进行计算return np.where(x >= c,0,np.where(x < c0,x / c0 * hc,(c - x) / (c - c0) * hc))

由于三角波形分为三段，因此需要两个嵌套的where( )进行计算。由于所有的运算和循环都在C语言级别完成，因此它的计算效率比frompyfunc( )高。

随着分段函数的分段数量的增加，需要嵌套更多层where( )。这样不便于程序的编写和阅读。可以用select( )解决这个问题，它的调用形式如下：

select(condlist, choicelist, default=0)

其中condlist是一个长度为N的布尔数组列表，choicelist是一个长度为N的存储候选值的数组列表，所有数组的长度都为M。如果列表元素不是数组而是单个数值，那么它相当于元素值都相同、长度为M的数组。

对于从0到M-1的数组下标i，从布尔数组列表中找出满足条件condlist[j][i]==True的j的最小值，则out[i]=choicelist[j][i]，其中out是select( )的返回数组。我们可以使用select( )计算三角波形：

def triangle_wave2(x, c, c0, hc):x = x - x.astype(np.int)return np.select([x >= c, x < c0 , True ],[0 , x/c0*hc, (c-x)/(c-c0)*hc])

由于分段函数分为三段，因此每个列表都有三个元素。choicelist的最后一个元素为True，表示前面的所有条件都不满足时，将使用choicelist的最后一个数组中的值。也可以用default参数指定条件都不满足时的候选值数组：

return np.select([x>= c, x < c0 ],[0 , x / c0 * hc],default=(c-x)/(c-c0)*hc)

但是where( )和select( )的所有参数都需要在调用它们之前完成计算，因此NumPy会计算下面4个数组：

x >= c, x < c0, x / c0 * hc, (c - x) / (c -c0 ) * hc

在计算时还会产生许多保存中间结果的数组，因此如果输入的数组x很大，将会发生大量内存分配和释放。

为了解决这个问题，NumPy提供了piecewise( )专门用于计算分段函数，它的调用参数如下：

piecewise(x, condlist, funclist)

参数x是一个保存自变量值的数组，condlist是一个长度为M的布尔数组列表，其中的每个布尔数组的长度都和数组x相同。funclist是一个长度为M或M+1的函数列表，这些函数的输入和输出都是数组。它们计算分段函数中的每个片段。如果不是函数而是数值，则相当于返回此数值的函数。每个函数与condlist中下标相同的布尔数组对应，如果funclist的长度为M+1，则最后一个函数计算所有条件都为False时的值。下面是使用piecewise( )计算三角波形的程序：

def triangle_wave3(x, c, c0, hc):x = x - x.astype(np.int)return np.piecewise(x,[x >= c, x < c0],[0, # x>=clambda x: x / c0 * hc, # x<c0lambda x: (c - x) / (c - c0) * hc]) # else

使用piecewise( )的好处在于它只计算需要计算的值。因此在上面的例子中，表达式x/c0*hc和(c-x)/(c-c0)*hc只对输入数组x中满足条件的部分进行计算。下面运行前面定义的三个分段函数，并使用%timeit命令比较这三个函数的运行时间：

x = np.linspace(0, 2, 10000)y1 = triangle_wave1(x, 0.6, 0.4, 1.0)y2 = triangle_wave2(x, 0.6, 0.4, 1.0)y3 = triangle_wave3(x, 0.6, 0.4, 1.0)np.all(y1 == y2), np.all(y1 == y3)(True, True)%timeit triangle_wave1(x, 0.6, 0.4, 1.0)%timeit triangle_wave2(x, 0.6, 0.4, 1.0)%timeit triangle_wave3(x, 0.6, 0.4, 1.0)1000 loops, best of 3: 614 µs per loop1000 loops, best of 3: 736 µs per loop1000 loops, best of 3: 311 µs per loop

六、操作多维数组

本节介绍的函数如表6所示。

表6　本节要介绍的函数

函数名	功能	函数名	功能
concatenate	连接多个数组	vstack	沿第0轴连接数组
hstack	沿第1轴连接数组	column_stack	按列连接多个一维数组
split、array_split	将数组分为多段	transpose	重新设置轴的顺序
swapaxes	交换两个轴的顺序

concatenate( )是连接多个数组的最基本的函数，其他函数都是它的快捷版本。它的第一个参数是包含多个数组的序列，它将沿着axis参数指定的轴（默认为第0轴）连接数组。所有这些数组的形状除了第axis轴之外都相同。

vstack( )沿着第0轴连接数组，当被连接的数组是长度为N的一维数组时，将其形状改为(1, N)。

hstack( )沿着第1轴连接数组。当所有数组都是一维时，沿着第0轴连接数组，因此结果数组仍然为一维的。

column_stack( )和hstack( )类似，沿着第1轴连接数组，但是当数组为一维时，将其形状改为(N, 1)，经常用于按列连接多个一维数组。

此外，c_[ ]对象也可以用于按列连接数组：

np.c_[a, b, a+b]array([[ 0, 10, 10],[ 1, 11, 12],[ 2, 12, 14]])

split( )和array_split( )的用法基本相同，将一个数组沿着指定轴分成多个数组，可以直接指定切分轴上的切分点下标。下面的代码把随机数组a切分为多个数组，保证每个数组中的元素都是升序排列的。注意通过diff( )和nonzero( )获得的下标是每个升序片段中最后一个元素的下标，而切分点为每个片段第一个元素的下标，因此需要+1。

当第二个参数为整数时，表示分组个数。split( )只能平均分组，而array_split( )能尽量平均分组：

np.split(a, 6) np.array_split(a, 5)--------------- --------------------[array([6, 3]), [array([6, 3, 7]),array([7, 4]), array([4, 6, 9]),array([6, 9]), array([2, 6]),array([2, 6]), array([7, 4]),array([7, 4]), array([3, 7])]array([3, 7])]

transpose( )和swapaxes( )用于修改轴的顺序，它们得到的是原数组的视图。transpose( )通过其第二个参数axes指定轴的顺序，默认时表示将整个形状翻转。而swapaxes( )通过两个整数指定调换顺序的轴。在下面的例子中：

• transpose( )的结果数组的形状为(3, 4, 2, 5)，它们分别位于原数组形状(2, 3, 4, 5)的(1, 2, 0, 3)下标位置处。
• swapaxes( )的结果数组的形状为(2, 4, 3, 5)，它是通过将原数组形状的中间两个轴对调得到的。

a = np.random.randint(0, 10, (2, 3, 4, 5))print u'原数组形状:', a.shapeprint u'transpose:', np.transpose(a, (1, 2, 0, 3)).shapeprint u'swapaxes:', np.swapaxes(a, 1, 2).shape原数组形状: (2, 3, 4, 5)transpose: (3, 4, 2, 5)swapaxes: (2, 4, 3, 5)

下面以将多个缩略图拼成一幅大图为例，帮助读者理解多维数组中变换轴的顺序。在data/thumbnails目录之下有30个160×90像素的PNG图标图像，需要将这些图像拼成一幅6行5列的大图像。首先调用glob和cv2模块中的函数，获得一个数组列表imgs。

import globimport numpy as npfrom cv2 import imread, imwriteimgs = [ ]for fn in glob.glob('thumbnails/*.png'):imgs.append(imread(fn, -1))print imgs[0].shape(90, 160, 3)

imgs中每个元素都是一个多维数组，它的形状为(90, 160, 3)，其中第0轴的长度为图像的高度，第1轴的长度为图像的宽度，第2轴为图像的通道数，彩色图像包含红、绿、蓝三个通道，所以第2轴的长度为3。

调用concatenate( )将这些数组沿第0轴拼成一个大数组，结果img是一个宽为160像素、高为2700像素的图像：

img = np.concatenate(imgs, 0)img.shape(2700, 160, 3)

由于我们的最终目标是把它们拼成一幅如图3（左）所示的6行5列的缩略图，因此需要将img的第0轴分解为3个轴，长度分别为(6, 5, 90)。下面使用reshape( )完成这个工作。使用img1[i, j]可以获取第i行、第j列上的图像：

img1 = img.reshape(6, 5, 90, 160, 3)img1[0, 1].shape(90, 160, 3)

图3　使用操作多维数组的函数拼接多幅缩略图

根据目标图像的大小，可以算出目标数组的形状为(540, 800, 3)，即(6*90, 5*160, 3)，也可以把它看作形状为(6, 90, 5, 160, 3)的多维数组。与img1的形状相比，可以看出需要交换img1的第1轴和第2轴。这个操作可以通过img1.swapaxes( )或img1.transpose( )完成。然后再通过reshape( )将数组的形状改为(540, 800, 3)。

img2 = img1.swapaxes(1, 2).reshape(540, 800, 3)

请读者思考下面的img3会得到怎样的图像：

img = np.concatenate(imgs, 0)img3 = img.reshape(5, 6, 90, 160, 3) \.transpose(1, 2, 0, 3, 4) \.reshape(540, 800, 3)

下面的程序将每幅缩略图的边沿上的两个像素填充为白色，效果如图3（右）所示。❶这里使用一个形状与img1的前4个轴相同的mask布尔数组，该数组的初始值为True。❷通过切片将mask中除去边框的部分设置为False。❸将img1中与mask为True的对应像素填充为白色。

img = np.concatenate(imgs, 0)img1 = img.reshape(6, 5, 90, 160, 3)mask = np.ones(img1.shape[:-1], dtype=bool) ❶mask[:, :, 2:-2, 2:-2] = False ❷img1[mask] = 230 ❸img4 = img1.swapaxes(1, 2).reshape(540, 800, 3)

七、多项式函数

多项式函数是变量的整数次幂与系数的乘积之和，可以用下面的数学公式表示：

由于多项式函数只包含加法和乘法运算，因此它很容易计算，可用于计算其他数学函数的近似值。多项式函数的应用非常广泛，例如在嵌入式系统中经常会用它计算正弦、余弦等函数。在NumPy中，多项式函数的系数可以用一维数组表示，例如可以用下面的数组表示，其中a[0]是最高次的系数，a[-1]是常数项，注意x2的系数为0。

a = np.array([1.0, 0, -2, 1])

我们可以用poly1d( )将系数转换为poly1d（一元多项式）对象，此对象可以像函数一样调用，它返回多项式函数的值：

p = np.poly1d(a)print type(p)p(np.linspace(0, 1, 5))<class 'numpy.lib.polynomial.poly1d'>array([ 1. , 0.515625, 0.125 , -0.078125, 0. ])

对poly1d对象进行加减乘除运算相当于对相应的多项式函数进行计算。例如：

p + [-2, 1] # 和 p + np.poly1d([-2, 1]) 相同poly1d([ 1., 0., -4., 2.])p * p # 两个3次多项式相乘得到一个6次多项式poly1d([ 1., 0., -4., 2., 4., -4., 1.])p / [1, 1] # 除法返回两个多项式，分别为商式和余式(poly1d([ 1., -1., -1.]), poly1d([ 2.]))

由于多项式的除法不一定能正好整除，因此它返回除法所得到的商式和余式。在上面的例子中，商式为<im><??-??>，余式为2。因此将商式和被除式相乘，再加上余式就等于原来的p：

p == np.poly1d([ 1., -1., -1.]) * [1,1] + 2True

多项式对象的deriv( )和integ( )方法分别计算多项式函数的微分和积分：

p.deriv( )poly1d([ 3., 0., -2.])p.integ( )poly1d([ 0.25, 0. , -1. , 1. , 0. ])p.integ( ).deriv( ) == pTrue

多项式函数的根可以使用roots( )函数计算：

r = np.roots(p)rarray([-1.61803399, 1. , 0.61803399])p(r) # 将根带入多项式计算，得到的值近似为0array([ 2.33146835e-15, 1.33226763e-15, 1.11022302e-16])

而poly( )函数可以将根转换回多项式的系数：

np.poly(r)array([ 1.00000000e+00, -1.66533454e-15, -2.00000000e+00,1.00000000e+00])

除了使用多项式对象之外，也可以直接使用NumPy提供的多项式函数对表示多项式系数的数组进行运算。可以在IPython中使用自动补全查看函数名：

>>> np.poly # 按Tab键np.poly np.polyadd np.polydiv np.polyint np.polysubnp.poly1d np.polyder np.polyfit np.polymul np.polyval

其中的polyfit( )函数可以对一组数据使用多项式函数进行拟合，找到与这组数据的误差平方和最小的多项式的系数。下面的程序用它计算-π/2～π/2区间与sin(x)函数最接近的多项式的系数：

np.set_printoptions(suppress=True, precision=4)x = np.linspace(-np.pi / 2, np.pi / 2, 1000) ❶y = np.sin(x) ❷for deg in [3, 5, 7]:a = np.polyfit(x, y, deg) ❸error = np.abs(np.polyval(a, x) - y) ❹print 'degree {}: {}'.format(deg, a)print 'max error of order %d:' % deg, np.max(error)degree 3: [-0.145 -0. 0.9887 0. ]max error of order 3: 0.00894699376707degree 5: [ 0.0076 -0. -0.1658 -0. 0.9998 -0. ]max error of order 5: 0.000157408614187degree 7: [-0.0002 -0. 0.0083 0. -0.1667 -0. 1. 0. ]max error of order 7: 1.52682558063e-06

❶首先通过linspace( )将-π/2～π/2区间等分为(1000-1)等份。❷计算拟合目标函数sin(x)的值。❸将表示目标函数的数组传递给polyfit( )进行拟合，第三个参数deg为多项式函数的最高阶数。polyfit( )所得到的多项式和目标函数在给定的1000个点之间的误差最小，polyfit( )返回多项式的系数数组。❹使用polyval( )计算多项式函数的值，并计算与目标函数的差的绝对值。

从程序的输出可以看到，由于正弦函数是一个奇函数，因此拟合的多项式系数中偶数次项的系数接近于0。图4显示了各阶多项式与正弦函数之间的误差，请注意图中Y轴为对数坐标。

图4　各阶多项式近似正弦函数的误差

八、多项式函数类

numpy.polynomial模块中提供了更丰富的多项式函数类，例如Polynomial、Chebyshev、Legendre等。它们和前面介绍的numpy.poly1d相反，多项式各项的系数按照幂从小到大的顺序排列，下面使用Polynomial类表示多项式x3-2x＋1，并计算x=2处的值：

from numpy.polynomial import Polynomial, Chebyshevp = Polynomial([1, -2, 0, 1])print p(2.0)5.0

Polynomial对象提供了众多的方法对多项式进行操作，例如deriv( )计算导函数：

p.deriv( )Polynomial([-2., 0., 3.], [-1., 1.], [-1., 1.])

切比雪夫多项式是一个正交多项式序列Ti(x)，一个n次多项式可以表示为多个切比雪夫多项式的加权和。在NumPy中，使用Chebyshev类表示由切比雪夫多项式组成的多项式p(x)：

Ti(x)多项式可以通过Chebyshev.basis(i)获得，图5显示了0到4次切比雪夫多项式。通过多项式类的convert( )方法可以在不同类型的多项式之间相互转换，转换的目标类型由kind参数指定。例如下面将T4(x)转换成Polynomial类。由结果可知：

Chebyshev.basis(4).convert(kind=Polynomial)Polynomial([ 1., 0., -8., 0., 8.], [-1., 1.], [-1., 1.])

图5　0到4次切比雪夫多项式

切比雪夫多项式的根被称为切比雪夫节点，可以用于多项式插值。相应的插值多项式能最大限度地降低龙格现象，并且提供多项式在连续函数的最佳一致逼近。下面以

函数插值为例演示切比雪夫节点与龙格现象。

❶在[-1,1]区间上等距离取n个取样点。❷使用n阶切比雪夫多项式的根作为取样点。❸使用两种取样点分别对f(x)进行多项式插值，即计算一个多项式经过所有的插值点。图6显示了两种插值点所得到的插值多项式，由左图可知等距离插值多项式在两端有非常大的振荡，这种现象被称为龙格现象，n越大振荡也越大；而右图采用切比雪夫节点作为插值点，插值多项式的振荡明显减小，并且n越大振荡越小。

插值与拟合

所谓多项式插值就是找到一个多项式经过所有的插值点。一个n阶多项式有n+1个系数，因此可以通过解方程求解经过n+1个插值点的n阶多项式的系数。fit( )方法虽然计算与目标点拟合的多项式系数，但是当使用n阶多项式拟合n+1的目标点时，多项式将经过所有目标点，因此其结果与多项式插值相同。

def f(x):return 1.0 / ( 1 + 25 * x**2)n = 11x1 = np.linspace(-1, 1, n) ❶x2 = Chebyshev.basis(n).roots( ) ❷xd = np.linspace(-1, 1, 200)c1 = Chebyshev.fit(x1, f(x1), n - 1, domain=[-1, 1]) ❸c2 = Chebyshev.fit(x2, f(x2), n - 1, domain=[-1, 1])print u'插值多项式的最大误差：',print u'等距离取样点：', abs(c1(xd) - f(xd)).max( ),print u'切比雪夫节点：', abs(c2(xd) - f(xd)).max( )插值多项式的最大误差：等距离取样点： 1.91556933029 切比雪夫节点： 0.109149825014

图6　等距离插值点（左）、切比雪夫插值点（右）

在使用多项式逼近函数时，使用切比雪夫多项式进行插值的误差比一般多项式要小许多。在下面的例子中，对g(x)在100个切比雪夫节点之上分别使用Polynomial和Chebyshev进行插值，结果如图7所示。在使用Polynomial.fit( )插值时，产生了RankWarning: The fit may be poorly conditioned警告，因此其结果多项式未能经过所有插值点。

def g(x):x = (x - 1) * 5return np.sin(x**2) + np.sin(x)**2n = 100x = Chebyshev.basis(n).roots( )xd = np.linspace(-1, 1, 1000)p_g = Polynomial.fit(x, g(x), n - 1, domain=[-1, 1])c_g = Chebyshev.fit(x, g(x), n - 1, domain=[-1, 1])print 'Max Polynomial Error:', abs(g(xd) - p_g(xd)).max( )print 'Max Chebyshev Error:', abs(g(xd) - c_g(xd)).max( )Max Polynomial Error: 1.19560558744Max Chebyshev Error: 6.47575726376e-09

图7　Chebyshev插值与Polynomial插值比较

trim( )方法可以降低多项式的次数，将尾部绝对值小于参数tol的高次系数截断。下面使用trim( )方法获取c中前68个系数，得到一个新的Chebyshev对象c_trimed，其最大误差上升到0.09左右。

c_trimed = c_g.trim(tol=0.05)print 'degree:', c_trimed.degree( )print 'error:', abs(g(xd) - c_trimed(xd)).max( )degree: 68error: 0.0912094835458

下面用同样的方法对函数h(x)进行19阶的切比雪夫多项式插值，得到插值多项式c_h：

def h(x):x = 5 * xreturn np.exp(-x**2 / 10)n = 20x = Chebyshev.basis(n).roots( )c_h = Chebyshev.fit(x, h(x), n - 1, domain=[-1, 1])print 'Max Chebyshev Error:', abs(h(xd) - c_h(xd)).max( )Max Chebyshev Error: 1.66544267266e-09

多项式类支持四则运算，下面将c_g和c_h相减得到c_diff，并调用其roots( )计算其所有根。然后找出其中所有的实数根real_roots，它们就是g(x)与h(x)交点的横坐标。图8显示了这两条函数曲线以及通过插值多项式计算的交点：

c_diff = c_g - c_hroots = c_diff.roots( )real_roots = roots[roots.imag == 0].realprint np.allclose(c_diff(real_roots), 0)True

图8　使用Chebyshev插值计算两条曲线在[-1，1]之间的所有交点

切比雪夫多项式在区间[-1,1]上为正交多项式，因此只有在该区间才能对目标函数正确插值。为了对任何区域的目标函数进行插值，需要对自变量的区间进行缩放和平移变换。可以通过domain参数指定拟合点的区间。在下面的例子中，对g2(x)在区间[-10,0]之内使用切比雪夫多项式进行插值。❶为了产生目标区间的切比雪夫节点，在通过basis( )方法创建Tn(x)时，通过domain参数指定目标区间。❷在调用fit( )方法进行拟合时，通过domain参数指定同样的区间。❸最后输出拟合得到的c_g2多项式在[-10,0]中与目标函数的最大误差。

def g2(x):return np.sin(x**2) + np.sin(x)**2n = 100x = Chebyshev.basis(n, domain=[-10, 0]).roots( ) ❶xd = np.linspace(-10, 0, 1000)c_g2 = Chebyshev.fit(x, g2(x), n - 1, domain=[-10, 0]) ❷print 'Max Chebyshev Error:', abs(g2(xd) - c_g2(xd)).max( ) ❸Max Chebyshev Error: 6.47574571744e-09

九、各种乘积运算

本节介绍的函数如表9所示。

表9　本节要介绍的函数

函数名	功能	函数名	功能
dot	矩阵乘积	inner	内积
outter	外积	tensordot	张量乘积

矩阵的乘积可以使用dot( )计算。对于二维数组，它计算的是矩阵乘积；对于一维数组，它计算的是内积。当需要将一维数组当作列矢量或行矢量进行矩阵运算时，先将一维数组转换为二维数组：

a = np.array([1, 2, 3])a[:, None] a[None, :]---------- -----------[[1], [[1, 2, 3]][2],[3]]

对于多维数组，dot( )的通用计算公式如下，即结果数组中的每个元素都是：数组a的最后轴上的所有元素与数组b的倒数第二轴上的所有元素的乘积和：

dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])

下面以两个三维数组的乘积演示dot( )的计算结果。首先创建两个三维数组，这两个数组的最后两轴满足矩阵乘积的条件：

a = np.arange(12).reshape(2, 3, 2)b = np.arange(12, 24).reshape(2, 2, 3)c = np.dot(a, b)c.shape(2, 3, 2, 3)

c是数组a和b的多个子矩阵的乘积。我们可以把数组a看作两个形状为(3,2)的矩阵，而把数组b看作两个形状为(2,3)的矩阵。a中的两个矩阵分别与b中的两个矩阵进行矩阵乘积，就得到数组c，c[i, :, j, :]是a中第i个矩阵与b中第j个矩阵的乘积。

for i, j in np.ndindex(2, 2):assert np.alltrue( c[i, :, j, :] == np.dot(a[i], b[j]) )

对于两个一维数组，inner( )和dot( )一样，计算两个数组对应下标元素的乘积和。而对于多维数组，它计算的结果数组中的每个元素都是：数组a和b的最后轴的内积。因此数组a和b的最后轴的长度必须相同：

inner(a, b)[i,j,k,m] = sum(a[i,j,:]*b[k,m,:])

下面是对inner( )的演示：

a = np.arange(12).reshape(2, 3, 2)b = np.arange(12, 24).reshape(2, 3, 2)c = np.inner(a, b)c.shape(2, 3, 2, 3)for i, j, k, l in np.ndindex(2, 3, 2, 3):assert c[i, j, k, l] == np.inner(a[i, j], b[k, l])

outer( )只对一维数组进行计算，如果传入的是多维数组，则先将此数组展平为一维数组之后再进行运算。它计算列向量和行向量的矩阵乘积：

a = np.array([1, 2, 3])b = np.array([4, 5, 6, 7])np.outer(a, b) np.dot(a[:, None], b[None, :])------------------ ------------------------------[[ 4, 5, 6, 7], [[ 4, 5, 6, 7],[ 8, 10, 12, 14],[ 8, 10, 12, 14],[12, 15, 18, 21]][12, 15, 18, 21]]

tensordot( )将两个多维数组a和b指定轴上的对应元素相乘并求和，它是最一般化的乘积运算函数。下面通过一些例子逐步介绍其用法。下面计算两个矩阵的乘积：❶axes参数有两个元素，第一个元素表示a中的轴，第二个元素表示b中的轴，这两个轴上对应的元素相乘之后求和。❷axes也可以是一个整数，它表示把a中的后axes个轴和b中的前axes个轴进行乘积和运算，而对于乘积和之外的轴则保持不变。

a = np.random.rand(3, 4)b = np.random.rand(4, 5)c1 = np.tensordot(a, b, axes=[[1], [0]]) ❶c2 = np.tensordot(a, b, axes=1) ❷c3 = np.dot(a, b)assert np.allclose(c1, c3)assert np.allclose(c2, c3)

对于多维数组的dot( )乘积，可以用tensordot(a, b, axes=[[-1], [-2]])表示，即将a的最后轴和b中的倒数第二轴求乘积和：

a = np.arange(12).reshape(2, 3, 2)b = np.arange(12, 24).reshape(2, 2, 3)c1 = np.tensordot(a, b, axes=[[-1], [-2]])c2 = np.dot(a, b)assert np.alltrue(c1 == c2)

在下面的例子中，将a的第1、第2轴与b的第1、第0轴求乘积和，因此c中的每个元素都是按照如下表达式计算的：

c[i, j, k, l] = np.sum(a[i, :, :, j] * b[:, :, k, l].T)

注意由于b对应的axes中的轴是倒序的，因此需要做转置操作。

a = np.random.rand(4, 5, 6, 7)b = np.random.rand(6, 5, 2, 3)c = np.tensordot(a, b, axes=[[1, 2], [1, 0]])for i, j, k, l in np.ndindex(4, 7, 2, 3):assert np.allclose(c[i, j, k, l], np.sum(a[i, :, :, j] * b[:, :, k, l].T))c.shape(4, 7, 2, 3)

十、广义ufunc函

从NumPy 1.8开始正式支持广义ufunc函数(generalized ufunc，以下简称gufunc)。gufunc是对ufunc的推广，所谓ufunc就是将对单个数值的运算通过广播运用到整个数组中的所有元素之上，而gufunc则是将对单个矩阵的运算通过广播运用到整个数组之上。例如numpy.linalg.inv( )是求逆矩阵的gufunc函数。在其文档中描述其输入输出数组的形状如下：

ainv = inv(a)

a : (..., M, M)

ainv : (..., M, M)

输入数组a的形状中带有“...”，它表示0到任意多个轴。当它为空时，就是对单个矩阵求逆，gufunc函数将对这些轴进行广播运算。最后两个轴的长度为M，表示任意大小的方形矩阵。

NumPy中的线性代数模块linalg中提供的函数大都为广义ufunc函数。在SciPy中也提供了线性代数模块linalg，但其中的函数都是一般函数，只能对单个矩阵进行计算。关于线性代数函数库的用法将在下一章进行详细介绍。

在输出数组ainv中，由于逆矩阵的形状与原矩阵相同，因此ainv的最后两轴的形状也是(M,M)。“...”表示广播运算之后的形状，而由于矩阵求逆只对一个矩阵进行运算，因此“...”的形状和a中的“...”的形状相同。

在下面的例子中，a的形状为(10, 20, 3, 3)，其中(10, 20)与“...”对应，3与M对应。而inv( )通过广播运算对10×20个形状为(3, 3)的矩阵求逆。得到的结果ainv的形状与a相同，也是(10, 20, 3, 3)。

a = np.random.rand(10, 20, 3, 3)ainv = np.linalg.inv(a)ainv.shape(10, 20, 3, 3)

下面的程序验证第i行、第j列的矩阵及其逆矩阵的乘积，应该近似等于3阶单位矩阵

i, j = 3, 4np.allclose(np.dot(a[i, j], ainv[i, j]), np.eye(3))True

numpy.linalg.det( )计算矩阵的行列式，它也是一个gufunc函数。它的输入输出的形状为：

adet = det(a)

a : (..., M, M)

adet : (...)

由于矩阵的行列式是将一个M×M的矩阵映射到一个标量，因此输出adet的形状中只包含“...”。

adet = np.linalg.det(a)adet.shape(10, 20)

下面以多个二次函数的数据拟合为例，介绍如何使用gufunc函数提高运算效率。首先通过随机数函数创建测试用的数据x和y，这两个数组的形状都为(n,10)。其中的每行数据(x[i]和y[i])构成一个曲线拟合的数据集，它们的关系为：y=β2+β1x+β0x2。现在需要计算每对数据所对应的系数β。

n = 10000np.random.seed(0)beta = np.random.rand(n, 3)x = np.random.rand(n, 10)y = beta[:,2, None] + x*beta[:, 1, None] + x**2*beta[:, 0, None]

显然使用前面介绍过的numpy.polyfit( )可以很方便地完成这个任务，下面的程序输出第42组的实际系数以及拟合的结果：

print beta[42]print np.polyfit(x[42], y[42], 2)[ 0.0191932 0.30157482 0.66017354][ 0.0191932 0.30157482 0.66017354]

只需要循环调用n次numpy.polyfit( )即可得到所需的结果，但是它的运算速度有些慢：

%time beta2 = np.vstack([np.polyfit(x[i], y[i], 2) for i in range(n)])Wall time: 1.52 snp.allclose(beta, beta2)True

在numpy.polyfit( )内部实际上是通过调用最小二乘法函数numpy.linalg.lstsq( )来实现多项式拟合的，我们也可以直接调用lstsq( )计算系数：

xx = np.column_stack(([x[42]**2, x[42], np.ones_like(x[42])]))print np.linalg.lstsq(xx, y[42])[0][ 0.0191932 0.30157482 0.66017354]

但遗憾的是，目前numpy.linalg.lstsq( )还不是gufunc函数，因此无法直接使用它计算所有数据组的拟合系数。然而numpy.linalg中对线性方程组求解的函数solve( )是一个gufunc函数。并且根据最小二乘法的公式：

只需要求出XTX和XTy，就可以使用numpy.linalg.solve( )计算出

。为了实现这个运算，还需要计算矩阵乘积的gufunc函数。然而dot( )并不是一个gufunc函数，因为它不遵循广播规则。NumPy中目前还没有正式提供计算矩阵乘积的gufunc函数，不过在umath_tests模块中提供了一个测试用的函数：matrix_multiply( )。下面的程序使用它和solve( )实现高速多项式拟合运算，它所需的时间约为polyfit( )版本的五十分之一。

%%timeX = np.dstack([x**2, x, np.ones_like(x)])Xt = X.swapaxes(-1, -2)import numpy.core.umath_tests as umathA = umath.matrix_multiply(Xt, X)b = umath.matrix_multiply(Xt, y[..., None]).squeeze( )beta3 = np.linalg.solve(A, b)print np.allclose(beta3, beta2)TrueWall time: 30 ms

在上面的运算中，X的形状为(10000, 10, 3)，Xt的形状为(10000, 3, 10)。matrix_multiply( )的各个参数和返回值的形状如下

c = matrix_multiply(a, b)

a : (..., M, N)

b : (..., N, K)

c : (..., M, K)

调用matrix_multiply( )对Xt和X中的每对矩阵进行乘积运算，得到的结果A的形状为(10000, 3, 3)。而为了计算XTy，需要通过y[..., None]将y变成形状为(10000,10,1)的数组。matrix_multiply(Xt, y[..., None])所得到的形状为(10000, 3, 1)。调用其squeeze( )，删除长度为1的轴。这样b的形状为(10000, 3)。

solve( )的参数b支持两种形状，其中第一种情况的形状如下：

x = solve(a, b)

a : (..., M, M)

b : (..., M)

x : (..., M)

因此solve( )的返回值beta3的形状为(10000, 3)。

前面的例子中，使用的都是最简单的广播规则。实际上gufunc函数支持所有的ufunc函数的广播规则。因此形状分别为(a, m, n)和(b, 1, n, k)的两个数组通过matrix_multiply( )乘积之后得到的数组的形状为(b, a, m, k)。下面看一个使用gufunc函数广播运算的例子：

在二维平面上的旋转矩阵为：

它能将平面上的某点的坐标围绕原点旋转θ。对于形状为(N, 2)的矩阵P，可以表示平面上N个点的坐标。而矩阵乘积得到的则是将这N个点绕坐标原点旋转θ之后的坐标。下面的程序使用matrix_multiply( )将3条曲线上的坐标点分别旋转4个角度，得到12条曲线。

调用matrix_multiply( )时两个参数数组的形状分别为(3, 100, 2)和(4, 1, 2, 2)，其中广播轴的形状分别为(3,)和(4, 1)，运算轴的形状分别为(100, 2)和(2, 2)。广播轴进行广播之后的形状为(4, 3)，而运算轴进行矩阵乘积之后的形状为(100, 2)，因此结果rpoints的形状为(4, 3, 100, 2)。

M = np.array([[[np.cos(t), -np.sin(t)],[np.sin(t), np.cos(t)]]for t in np.linspace(0, np.pi, 4, endpoint=False)])x = np.linspace(-1, 1, 100)points = np.array((np.c_[x, x], np.c_[x, x**3], np.c_[x**3, x]))rpoints = umath.matrix_multiply(points, M[:, None, ...])print points.shape, M.shape, rpoints.shape(3, 100, 2) (4, 2, 2) (4, 3, 100, 2)

将这12条曲线绘制成图表之后的效果如图9所示。

图9　使用矩阵乘积的广播运算将3条曲线分别旋转4个角度