本人以简书作者 SeanCheney 系列专题文章并结合原书为学习资源,记录个人笔记,仅作为知识记录及后期复习所用,原作者地址查看 简书 SeanCheney,如有错误,还望批评指教。——ZJ
原作者:SeanCheney | 《利用Python进行数据分析·第2版》第4章 NumPy基础:数组和矢量计算 | 來源:简书
环境:Python 3.6
第 4 章 NumPy 基础:数组和矢量计算
NumPy(Numerical Python 的简称)是 Python数 值计算最重要的基础包。大多数提供科学计算的包都是用 NumPy 的数组作为构建基础。
NumPy 的部分功能如下:
- ndarray,一个具有矢量算术运算和复杂广播能力的快速且节省空间的多维数组。
- 用于对整组数据进行快速运算的标准数学函数(无需编写循环)。
- 用于读写磁盘数据的工具以及用于操作内存映射文件的工具。
- 线性代数、随机数生成以及傅里叶变换功能。
- 用于集成由 C、C++、Fortran 等语言编写的代码的 A C API。
对于大部分数据分析应用而言,我最关注的功能主要集中在:
- 用于数据整理和清理、子集构造和过滤、转换等快速的矢量化数组运算。
- 常用的数组算法,如排序、唯一化、集合运算等。
- 高效的描述统计和数据聚合/摘要运算。
- 用于异构数据集的合并/连接运算的数据对齐和关系型数据运算。
- 将条件逻辑表述为数组表达式(而不是带有if-elif-else分支的循环)。
- 数据的分组运算(聚合、转换、函数应用等)
NumPy 之于数值计算特别重要的原因之一,是因为它可以高效处理大数组的数据。这是因为:
- NumPy 是在一个连续的内存块中存储数据,独立于其他 Python 内置对象。NumPy 的 C 语言编写的算法库可以操作内存,而不必进行类型检查或其它前期工作。比起 Python 的内置序列,NumPy 数组使用的内存更少。
- NumPy 可以在整个数组上执行复杂的计算,而不需要 Python 的 for 循环。
要搞明白具体的性能差距,考察一个包含一百万整数的数组,和一个等价的Python 列表:
In [1]: import numpy as np
In [2]: my_arr = np.arange(1000000)
In [3]: my_list = list(range(1000000))
# 各个序列分别乘以 2:
In [4]: %time for _ in range(10):my_arr2 = my_arr * 2
Wall time: 46.9 ms
In [5]: %time for _ in range(10): my_list = [x * 2 for x in my_list]
Wall time: 913 ms
基于 NumPy 的算法要比纯 Python 快 10 到 100 倍(甚至更快),并且使用的内存更少。
4.1 NumPy 的 ndarray:一种多维数组对象
NumPy 最重要的一个特点就是其 N 维数组对象(即 ndarray),该对象是一个快速而灵活的大数据集容器。
Python 是如何利用与标量值类似的语法进行批次计算:
In [6]: data = np.random.randn(2,3)
In [7]: data
Out[7]:
array([[-0.7949888 , 0.59269599, 0.6047145 ],
[-0.05579778, -1.97623901, -0.52438662]])
In [8]: data * 10
Out[8]:
array([[ -7.94988799, 5.92695992, 6.04714496],
[ -0.55797776, -19.76239014, -5.24386616]])
In [9]: data + data
Out[9]:
array([[-1.5899776 , 1.18539198, 1.20942899],
[-0.11159555, -3.95247803, -1.04877323]])
ndarray 是一个通用的同构数据多维容器,也就是说,其中的所有元素必须是相同类型的。每个数组都有一个 shape(一个表示各维度大小的元组)和一个dtype(一个用于说明数组数据类型的对象):
In [10]: data.shape
Out[10]: (2, 3)
In [11]: data.dtype
Out[11]: dtype('float64')
- 笔记:当你在本书中看到“数组”、“NumPy 数组”、”ndarray”时,基本上都指的是同一样东西,即 ndarray 对象。
创建 ndarray
创建数组最简单的办法就是使用 array 函数。它接受一切序列型的对象(包括其他数组),然后产生一个新的含有传入数据的 NumPy 数组。以一个列表的转换为例:
In [12]: data1 = [6, 7.5, 8, 0, 1]
In [13]: arr1 = np.array(data1)
In [14]: arr1
Out[14]: array([6. , 7.5, 8. , 0. , 1. ])
# 嵌套序列(比如由一组等长列表组成的列表)将会被转换为一个多维数组:
In [15]: data2 = [[1,2,3,4],[5,6,7,8]]
In [17]: arr2 = np.array(data2)
In [18]: arr2
Out[18]:
array([[1, 2, 3, 4],
[5, 6, 7, 8]])
因为 data2 是列表的列表,NumPy 数组 arr2 的两个维度的 shape 是从 data2引入的。可以用属性 ndim 和 shape 验证:
In [19]: arr2.ndim
Out[19]: 2
In [20]: arr2.shape
Out[20]: (2, 4)
除非特别说明(稍后将会详细介绍),np.array 会尝试为新建的这个数组推断出一个较为合适的数据类型。数据类型保存在一个特殊的dtype对象中。比如说,在上面的两个例子中,我们有:
In [21]: arr1.dtype
Out[21]: dtype('float64')
In [22]: arr2.dtype
Out[22]: dtype('int32')
除 np.array 之外,还有一些函数也可以新建数组。比如,zeros 和 ones 分别可以创建指定长度或形状的全 0 或全 1 数组。empty 可以创建一个没有任何具体值的数组。要用这些方法创建多维数组,只需传入一个表示形状的元组即可:
In [26]: np.zeros(10)
Out[26]: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [27]: np.zeros((4,5))
Out[27]:
array([[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.]])
In [28]: np.empty((2,3,2))
Out[28]:
array([[[0., 0.],
[0., 0.],
[0., 0.]],
[[0., 0.],
[0., 0.],
[0., 0.]]])
注意:认为np.empty
会返回全 0 数组的想法是不安全的。很多情况下(如前所示),它返回的都是一些未初始化的垃圾值。
- arange 是 Python 内置函数 range 的数组版:
In [29]: np.arange(12)
Out[29]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
表 4-1 列出了一些数组创建函数。由于 NumPy 关注的是数值计算,因此,如果没有特别指定,数据类型基本都是 float64(浮点数)。
ndarray 的数据类型
- dtype(数据类型)是一个特殊的对象,它含有 ndarray 将一块内存解释为特定数据类型所需的信息:
In [31]: arr1 = np.array([1,2,3], dtype=np.float64)
In [32]: arr2 = np.array([1,2,3], dtype=np.int32)
In [33]: arr1
Out[33]: array([1., 2., 3.])
In [34]: arr2
Out[34]: array([1, 2, 3])
In [35]: arr1.dtype
Out[35]: dtype('float64')
In [36]: arr2.dtype
Out[36]: dtype('int32')
dtype 是 NumPy 灵活交互其它系统的源泉之一。多数情况下,它们直接映射到相应的机器表示,这使得“读写磁盘上的二进制数据流”以及“集成低级语言代码(如 C、Fortran)”等工作变得更加简单。数值型 dtype 的命名方式相同:一个类型名(如 float 或 int),后面跟一个用于表示各元素位长的数字。标准的双精度浮点值(即 Python 中的 float对象)需要占用 8 字节(即64 位)。因此,该类型在 NumPy 中就记作 float64。表 4-2 列出了 NumPy 所支持的全部数据类型。
笔记:记不住这些 NumPy 的 dtype 也没关系,新手更是如此。通常只需要知道你所处理的数据的大致类型是浮点数、复数、整数、布尔值、字符串,还是普通的 Python 对象即可。当你需要控制数据在内存和磁盘中的存储方式时(尤其是对大数据集),那就得了解如何控制存储类型。
- 你可以通过 ndarray 的 astype 方法明确地将一个数组从一个 dtype 转换成另一个 dtype:
In [37]: arr = np.array([1,2,3,4,5])
In [38]: arr.dtype
Out[38]: dtype('int32')
In [39]: float_arr = arr.astype(np.float64)
In [40]: float_arr.dtype
Out[40]: dtype('float64')
- 在本例中,整数被转换成了浮点数。如果将浮点数转换成整数,则小数部分将会被截取删除:
In [41]: arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
In [42]: arr
Out[42]: array([ 3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
In [43]: arr.astype(np.int32)
Out[43]: array([ 3, -1, -2, 0, 12, 10])
- 如果某字符串数组表示的全是数字,也可以用 astype 将其转换为数值形式:
In [44]: numeric_strings = np.array(['1.25', '9.6', '45'], dtype=np.s
...: tring_ )
In [45]: numeric_strings.astype(float)
Out[45]: array([ 1.25, 9.6 , 45. ])
注意:使用
numpy.string_
类型时,一定要小心,因为 NumPy 的字符串数据是大小固定的,发生截取时,不会发出警告。pandas 提供了更多非数值数据的便利的处理方法。数组的 dtype 还有另一个属性:
In [46]: int_array = np.arange(10)
In [47]: calibers = np.array([.22, .270, .357, .380, .44, .50], dtype=np.float64)
In [48]: int_array.astype(calibers.dtype)
Out[48]: array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
# 你还可以用简洁的类型代码来表示 dtype:
In [49]: empty_uint32 = np.empty(8, dtype='u4')
In [50]: empty_uint32
Out[50]:
array([ 0, 1075314688, 0, 1075707904, 0,
1075838976, 0, 1072693248], dtype=uint32)
- 笔记:调用 astype 总会创建一个新的数组(一个数据的备份),即使新的dtype 与旧的 dtype 相同。
NumPy 数组的运算
数组很重要,因为它使你不用编写循环即可对数据执行批量运算。NumPy 用户称其为矢量化(vectorization)。大小相等的数组之间的任何算术运算都会将运算应用到元素级:
In [57]: arr = np.array([[1., 2., 3.], [4., 5., 6.]])
In [58]: arr
Out[58]:
array([[1., 2., 3.],
[4., 5., 6.]])
In [59]: arr * arr
Out[59]:
array([[ 1., 4., 9.],
[16., 25., 36.]])
In [60]: %time arr - arr
Wall time: 15.6 ms
Out[60]:
array([[0., 0., 0.],
[0., 0., 0.]])
# 数组与标量的算术运算会将标量值传播到各个元素:
In [61]: 1/arr
Out[61]:
array([[1. , 0.5 , 0.33333333],
[0.25 , 0.2 , 0.16666667]])
In [62]: arr ** 0.5
Out[62]:
array([[1. , 1.41421356, 1.73205081],
[2. , 2.23606798, 2.44948974]])
# 大小相同的数组之间的比较会生成布尔值数组:
In [63]: arr2 = np.array([[0., 4., 1.], [7., 2., 12.]])
In [64]: arr2
Out[64]:
array([[ 0., 4., 1.],
[ 7., 2., 12.]])
In [65]: arr2 > arr
Out[65]:
array([[False, True, False],
[ True, False, True]])
不同大小的数组之间的运算叫做广播(broadcasting),将在附录A中对其进行详细讨论。本书的内容不需要对广播机制有多深的理解。
基本的索引和切片
NumPy 数组的索引是一个内容丰富的主题,因为选取数据子集或单个元素的方式有很多。一维数组很简单。从表面上看,它们跟 Python 列表的功能差不多:
In [60]: arr = np.arange(10)
In [61]: arr
Out[61]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [62]: arr[5]
Out[62]: 5
In [63]: arr[5:8]
Out[63]: array([5, 6, 7])
In [64]: arr[5:8] = 12
In [65]: arr
Out[65]: array([ 0, 1, 2, 3, 4, 12, 12, 12, 8, 9])
如上所示,当你将一个标量值赋值给一个切片时(如 arr[5:8]=12
),该值会自动传播(也就说后面将会讲到的“广播”)到整个选区。跟列表最重要的区别在于,数组切片是原始数组的视图。这意味着数据不会被复制,视图上的任何修改都会直接反映到源数组上。
# 作为例子,先创建一个 arr 的切片:
In [73]: arr_slice = arr[5:8]
In [74]: arr_slice
Out[74]: array([12, 12, 12])
# 现在,当我修稿 arr_slice 中的值,变动也会体现在原始数组 arr 中:
In [75]: arr_slice[1] = 142345
In [76]: arr
Out[76]:
array([ 0, 1, 2, 3, 4, 12, 142345, 12,
8, 9])
# 切片[ : ]会给数组中的所有值赋值:
In [77]: arr_slice[:] = 64
In [78]: arr
Out[78]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
注意:如果你想要得到的是 ndarray 切片的一份副本而非视图,就需要明确地进行复制操作,例如
arr[5:8].copy()
。对于高维度数组,能做的事情更多。在一个二维数组中,各索引位置上的元素不再是标量而是一维数组:
In [79]: arr2d = np.array([[1,2,3],[4,5,6],[7,8,9]])
In [80]: arr2d[2]
Out[80]: array([7, 8, 9])
# 因此,可以对各个元素进行递归访问,但这样需要做的事情有点多。你可以传入一个以逗号隔开的索引列表来选取单个元素。也就是说,下面两种方式是等价的:
In [81]: arr2d[0][2]
Out[81]: 3
In [82]: arr2d[0,2]
Out[82]: 3
图4-1 说明了二维数组的索引方式。轴 0 作为行,轴 1 作为列。
在多维数组中,如果省略了后面的索引,则返回对象会是一个维度低一点的ndarray(它含有高一级维度上的所有数据)。因此,在 2×2×3 数组 arr3d中:
In [83]: arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 1
...: 1, 12]]])
# arr3d[0] 是一个 2×3 数组:
In [84]: arr3d[0]
Out[84]:
array([[1, 2, 3],
[4, 5, 6]])
# 标量值和数组都可以被赋值给arr3d[0]:
In [88]: old_values = arr3d[0].copy()
In [89]: arr3d[0] = 42
In [90]: arr3d
Out[90]:
array([[[42, 42, 42],
[42, 42, 42]],
[[ 7, 8, 9],
[10, 11, 12]]])
In [91]: arr3d[0] = old_values
In [92]: arr3d
Out[92]:
array([[[ 1, 2, 3],
[ 4, 5, 6]],
[[ 7, 8, 9],
[10, 11, 12]]])
# 相似的,arr3d[1,0] 可以访问索引以 (1,0) 开头的那些值(以一维数组的形式返回):
In [100]: arr3d[1,0]
Out[100]: array([7, 8, 9])
# 虽然是用两步进行索引的,表达式是相同的:
In [101]: x = arr3d[1]
In [102]: x
Out[102]:
array([[ 7, 8, 9],
[10, 11, 12]])
In [103]: x[0]
Out[103]: array([7, 8, 9])
注意,在上面所有这些选取数组子集的例子中,返回的数组都是视图。
切片索引
# ndarray 的切片语法跟 Python 列表这样的一维对象差不多:
In [4]: arr = np.arange(10)
In [5]: arr[1:6]
Out[5]: array([1, 2, 3, 4, 5])
# 对于之前的二维数组 arr2d,其切片方式稍显不同:
In [6]: arr2d = np.arr ay([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
In [7]: arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11
...: , 12]]])
In [8]: arr2d[:2]
Out[8]:
array([[1, 2, 3],
[4, 5, 6]])
# 可以看出,它是沿着第 0 轴(即第一个轴)切片的。也就是说,切片是沿着一
# 个轴向选取元素的。表达式 arr2d[:2] 可以被认为是“选取 arr2d 的前两行”。
# 你可以一次传入多个切片,就像传入多个索引那样
In [9]: arr2d[:2,1:]
Out[9]:
array([[2, 3],
[5, 6]])
# 像这样进行切片时,只能得到相同维数的数组视图。
# 通过将整数索引和切片混合,可以得到低维度的切片。
# 例如,我可以选取第二行的前两列:
In [10]: arr2d[1,:2]
Out[10]: array([4, 5])
# 相似的,还可以选择第三列的前两行:
In [11]: arr2d[:2,2]
Out[11]: array([3, 6])
In [12]: arr2d[:,:1]
Out[12]:
array([[1],
[4],
[7]])
自然,对切片表达式的赋值操作也会被扩散到整个选区:
In [96]: arr2d[:2, 1:] = 0
In [97]: arr2d
Out[97]:
array([[1, 0, 0],
[4, 0, 0],
[7, 8, 9]])
布尔型索引
来看这样一个例子,假设我们有一个用于存储数据的数组以及一个存储姓名的数组(含有重复项)。在这里,我将使用numpy.random
中的randn
函数生成一些正态分布的随机数据:
In [13]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe'
...: , 'Joe'])
In [14]: data = np.random.randn(7,4)
In [15]: names
Out[15]: array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dt
ype='<U4')
In [16]: data
Out[16]:
array([[-0.40297344, -0.04858888, 0.94186535, -0.43067234],
[ 0.94293832, 0.42824721, -0.47260321, -0.82525473],
[-1.15078413, 0.12698473, 1.22078788, -0.21719696],
[-3.39869444, -1.66794053, -0.93024575, -0.77330647],
[ 0.93072027, -0.53585592, -0.11522692, -1.38992789],
[-0.37161027, -0.91641934, -0.5625268 , -0.85193452],
[ 0.42138626, 0.80474264, -0.93333829, 0.54982325]])
# 假设每个名字都对应data数组中的一行,而我们想要选出对应于名字"Bob"的所有行。
# 跟算术运算一样,数组的比较运算(如==)也是矢量化的。
# 因此,对names和字符串"Bob"的比较运算将会产生一个布尔型数组:
In [17]: names == 'Bob'
Out[17]: array([ True, False, False, True, False, False, False])
# 这个布尔型数组可用于数组索引:
In [18]: data[names=='Bob']
Out[18]:
array([[-0.40297344, -0.04858888, 0.94186535, -0.43067234],
[-3.39869444, -1.66794053, -0.93024575, -0.77330647]])
布尔型数组的长度必须跟被索引的轴长度一致。此外,还可以将布尔型数组跟切片、整数(或整数序列,稍后将对此进行详细讲解)混合使用:
In [103]: data[names == 'Bob']
Out[103]:
array([[ 0.0929, 0.2817, 0.769 , 1.2464],
[ 1.669 , -0.4386, -0.5397, 0.477 ]])
注意:如果布尔型数组的长度不对,布尔型选择就会出错,因此一定要小心。
# 下面的例子,我选取了names == 'Bob'的行,并索引了列:
In [19]: data[names == 'Bob', 2:]
Out[19]:
array([[ 0.94186535, -0.43067234],
[-0.93024575, -0.77330647]])
In [20]: data[names == 'Bob', 3:]
Out[20]:
array([[-0.43067234],
[-0.77330647]])
In [21]: data
Out[21]:
array([[-0.40297344, -0.04858888, 0.94186535, -0.43067234],
[ 0.94293832, 0.42824721, -0.47260321, -0.82525473],
[-1.15078413, 0.12698473, 1.22078788, -0.21719696],
[-3.39869444, -1.66794053, -0.93024575, -0.77330647],
[ 0.93072027, -0.53585592, -0.11522692, -1.38992789],
[-0.37161027, -0.91641934, -0.5625268 , -0.85193452],
[ 0.42138626, 0.80474264, -0.93333829, 0.54982325]])
# 要选择除"Bob"以外的其他值,既可以使用不等于符号(!=),也可以通过~对条件进行否定:
In [22]: names != 'Bob'
Out[22]: array([False, True, True, False, True, True, True])
In [23]: data[~(names == 'Bob')]
Out[23]:
array([[ 0.94293832, 0.42824721, -0.47260321, -0.82525473],
[-1.15078413, 0.12698473, 1.22078788, -0.21719696],
[ 0.93072027, -0.53585592, -0.11522692, -1.38992789],
[-0.37161027, -0.91641934, -0.5625268 , -0.85193452],
[ 0.42138626, 0.80474264, -0.93333829, 0.54982325]])
# ~操作符用来反转条件很好用:
In [24]: cond = names == 'Bob'
In [25]: data[~cond]
Out[25]:
array([[ 0.94293832, 0.42824721, -0.47260321, -0.82525473],
[-1.15078413, 0.12698473, 1.22078788, -0.21719696],
[ 0.93072027, -0.53585592, -0.11522692, -1.38992789],
[-0.37161027, -0.91641934, -0.5625268 , -0.85193452],
[ 0.42138626, 0.80474264, -0.93333829, 0.54982325]])
In [26]: cond
Out[26]: array([ True, False, False, True, False, False, False])
In [27]: ~cond
Out[27]: array([False, True, True, False, True, True, True])
选取这三个名字中的两个需要组合应用多个布尔条件,使用 &(和)、|(或)之类的布尔算术运算符即可:
In [30]: mask = (names == 'Bob') | (names == 'Will')
In [31]: mask
Out[31]: array([ True, False, True, True, True, False, False])
In [32]: data[mask]
Out[32]:
array([[-0.40297344, -0.04858888, 0.94186535, -0.43067234],
[-1.15078413, 0.12698473, 1.22078788, -0.21719696],
[-3.39869444, -1.66794053, -0.93024575, -0.77330647],
[ 0.93072027, -0.53585592, -0.11522692, -1.38992789]])
通过布尔型索引选取数组中的数据,将总是创建数据的副本,即使返回一模一样的数组也是如此。
注意:Python 关键字 and 和 or 在布尔型数组中无效。要使用 & 与 |。
通过布尔型数组设置值是一种经常用到的手段。为了将 data 中的所有负值都设置为 0,我们只需:
In [33]: data[data < 0] =0
In [34]: data
Out[34]:
array([[0. , 0. , 0.94186535, 0. ],
[0.94293832, 0.42824721, 0. , 0. ],
[0. , 0.12698473, 1.22078788, 0. ],
[0. , 0. , 0. , 0. ],
[0.93072027, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. ],
[0.42138626, 0.80474264, 0. , 0.54982325]])
通过一维布尔数组设置整行或列的值也很简单:
In [35]: data[names != 'Joe'] = 7
In [36]: data
Out[36]:
array([[7. , 7. , 7. , 7. ],
[0.94293832, 0.42824721, 0. , 0. ],
[7. , 7. , 7. , 7. ],
[7. , 7. , 7. , 7. ],
[7. , 7. , 7. , 7. ],
[0. , 0. , 0. , 0. ],
[0.42138626, 0.80474264, 0. , 0.54982325]])
后面会看到,这类二维数据的操作也可以用 pandas
方便的来做。
花式索引
花式索引(Fancy indexing)是一个 NumPy 术语,它指的是利用整数数组进行索引。假设我们有一个8×4数组:
In [37]: arr = np.empty((8,4))
In [38]: for i in range(8):
...: arr[i] = i
...:
In [39]: arr
Out[39]:
array([[0., 0., 0., 0.],
[1., 1., 1., 1.],
[2., 2., 2., 2.],
[3., 3., 3., 3.],
[4., 4., 4., 4.],
[5., 5., 5., 5.],
[6., 6., 6., 6.],
[7., 7., 7., 7.]])
# 为了以特定顺序选取行子集,只需传入一个用于指定顺序的整数列表或ndarray 即可:
In [40]: arr[[4,3,0,6]]
Out[40]:
array([[4., 4., 4., 4.],
[3., 3., 3., 3.],
[0., 0., 0., 0.],
[6., 6., 6., 6.]])
# 这段代码确实达到我们的要求了!使用负数索引将会从末尾开始选取行:
In [41]: arr[[-3,-5,-7]]
Out[41]:
array([[5., 5., 5., 5.],
[3., 3., 3., 3.],
[1., 1., 1., 1.]])
# 一次传入多个索引数组会有一点特别。它返回的是一个一维数组,其中的元素对应各个索引元组:
In [42]: arr = np.arange(32).reshape((8,4))
In [43]: arr
Out[43]:
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]])
In [44]: arr[[1,5,7,2],[0,3,1,2]]
Out[44]: array([ 4, 23, 29, 10])
附录A中会详细介绍 reshape 方法。
最终选出的是元素 (1,0)、(5,3)、(7,1)和(2,2)。无论数组是多少维的,花式索引总是一维的。
这个花式索引的行为可能会跟某些用户的预期不一样(包括我在内),选取矩阵的行列子集应该是矩形区域的形式才对。下面是得到该结果的一个办法:
In [45]: arr[[1,5,7,2]][:,[0,3,1,2]]
Out[45]:
array([[ 4, 7, 5, 6],
[20, 23, 21, 22],
[28, 31, 29, 30],
[ 8, 11, 9, 10]])
- 记住,花式索引跟切片不一样,它总是将数据复制到新数组中。
数组转置和轴对换
转置是重塑的一种特殊形式,它返回的是源数据的视图(不会进行任何复制操作)。数组不仅有transpose
方法,还有一个特殊的 T 属性:
In [46]: arr = np.arange(15).reshape((3,5))
In [47]: arr
Out[47]:
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])
In [48]: arr.T
Out[48]:
array([[ 0, 5, 10],
[ 1, 6, 11],
[ 2, 7, 12],
[ 3, 8, 13],
[ 4, 9, 14]])
在进行矩阵计算时,经常需要用到该操作,比如利用np.dot
计算矩阵内积:
(3,6) * (6,3)–>(3, 3)
In [49]: arr = np.random.randn(6,3)
In [50]: arr
Out[50]:
array([[ 2.38335214, 0.09150081, -1.10938566],
[ 0.32481852, 1.5597611 , -1.3504366 ],
[-1.45644045, -0.54867809, 0.76936126],
[ 0.14373714, 0.96910097, 0.58995405],
[-0.43944108, 0.8263271 , -1.1597216 ],
[ 0.25505929, -0.90902731, 0.68217187]])
In [51]: np.dot(arr.T, arr)
Out[51]:
array([[ 8.18591734, 1.06815277, -3.43481047],
[ 1.06815277, 5.19057854, -3.63669704],
[-3.43481047, -3.63669704, 5.80469074]])
对于高维数组,transpose
需要得到一个由轴编号组成的元组才能对这些轴进行转置(比较费脑子):
In [52]: arr = np.arange(16).reshape((2,2,4))
In [53]: arr
Out[53]:
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7]],
[[ 8, 9, 10, 11],
[12, 13, 14, 15]]])
In [54]: arr.transpose((1,0,2))
Out[54]:
array([[[ 0, 1, 2, 3],
[ 8, 9, 10, 11]],
[[ 4, 5, 6, 7],
[12, 13, 14, 15]]])
这里,第一个轴被换成了第二个,第二个轴被换成了第一个,最后一个轴不变。
arr.transpose((1,0,2))
1,0,2 是索引位置
简单的转置可以使用.T,它其实就是进行轴对换而已。ndarray
还有一个swapaxes
方法,它需要接受一对轴编号:
In [55]: arr
Out[55]:
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7]],
[[ 8, 9, 10, 11],
[12, 13, 14, 15]]])
In [56]: arr.swapaxes(1,2)
Out[56]:
array([[[ 0, 4],
[ 1, 5],
[ 2, 6],
[ 3, 7]],
[[ 8, 12],
[ 9, 13],
[10, 14],
[11, 15]]])
swapaxes
也是返回源数据的视图(不会进行任何复制操作)。
4.2 通用函数:快速的元素级数组函数
通用函数(即 ufunc
)是一种对ndarray
中的数据执行元素级运算的函数。你可以将其看做简单函数(接受一个或多个标量值,并产生一个或多个标量值)的矢量化包装器。
许多ufunc
都是简单的元素级变体,如sqrt
和exp
:
In [57]: arr = np.arange(10)
In [58]: arr
Out[58]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [59]: np.sqrt(arr)
Out[59]:
array([0. , 1. , 1.41421356, 1.73205081, 2. ,
2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
In [60]: np.exp(arr)
Out[60]:
array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01,
5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03,
2.98095799e+03, 8.10308393e+03])
这些都是一元 (unary)ufunc
。另外一些(如add
或maximum
)接受 2 个数组(因此也叫二元(binary)ufunc
),并返回一个结果数组:
In [61]: x = np.random.randn(8)
In [62]: y = np.random.randn(8)
In [63]: x
Out[63]:
array([-0.49658376, 0.34930792, -0.98828359, 1.33465412, -0.81113981
,
1.50500909, -0.2623108 , 0.45881796])
In [64]: y
Out[64]:
array([ 1.90256952, -1.87612173, -0.32216431, 2.31260537, -0.97106547
,
0.94585092, 0.62071541, 0.27382797])
In [65]: np.maximum(x, y)
Out[65]:
array([ 1.90256952, 0.34930792, -0.32216431, 2.31260537, -0.81113981
,
1.50500909, 0.62071541, 0.45881796])
这里,numpy.maximum
计算了 x 和 y 中元素级别最大的元素。
虽然并不常见,但有些ufunc
的确可以返回多个数组。modf
就是一个例子,它是 Python 内置函数divmod
的矢量化版本,它会返回浮点数数组的小数和整数部分:
In [66]: arr = np.random.randn(7) * 5
In [67]: arr
Out[67]:
array([-1.62682089, -2.86349396, 3.07305674, 5.86658913, 4.95387596
,
4.22672489, 13.27870532])
In [68]: remainder, whole_part = np.modf(arr)
In [69]: remainder
Out[69]:
array([-0.62682089, -0.86349396, 0.07305674, 0.86658913, 0.95387596
,
0.22672489, 0.27870532])
In [70]: whole_part
Out[70]: array([-1., -2., 3., 5., 4., 4., 13.])
Ufuncs
可以接受一个out
可选参数,这样就能在数组原地进行操作:
In [71]: arr
Out[71]:
array([-1.62682089, -2.86349396, 3.07305674, 5.86658913, 4.95387596
,
4.22672489, 13.27870532])
In [72]: np.sqrt(arr)
Out[72]:
array([ nan, nan, 1.75301362, 2.42210428, 2.22573043,
2.05590002, 3.64399579])
In [73]: np.sqrt(arr,arr)
Out[73]:
array([ nan, nan, 1.75301362, 2.42210428, 2.22573043,
2.05590002, 3.64399579])
In [74]: arr
Out[74]:
array([ nan, nan, 1.75301362, 2.42210428, 2.22573043,
2.05590002, 3.64399579])
表 4-3 和表 4-4 分别列出了一些一元和二元 ufunc。
4.3 利用数组进行数据处理
NumPy 数组使你可以将许多种数据处理任务表述为简洁的数组表达式(否则需要编写循环)。用数组表达式代替循环的做法,通常被称为矢量化。一般来说,矢量化数组运算要比等价的纯 Python 方式快上一两个数量级(甚至更多),尤其是各种数值计算。在后面内容中(见附录A)我将介绍广播,这是一种针对矢量化计算的强大手段。
作为简单的例子,假设我们想要在一组值(网格型)上计算函数sqrt(x^2+y^2)
。np.meshgrid
函数接受两个一维数组,并产生两个二维矩阵(对应于两个数组中所有的(x,y)对):
In [75]: points = np.arange(-5, 5, 0.01) # 1000 equally spaced points
In [76]: xs, ys = np.meshgrid(points, points)
In [77]: ys
Out[77]:
array([[-5. , -5. , -5. , ..., -5. , -5. , -5. ],
[-4.99, -4.99, -4.99, ..., -4.99, -4.99, -4.99],
[-4.98, -4.98, -4.98, ..., -4.98, -4.98, -4.98],
...,
[ 4.97, 4.97, 4.97, ..., 4.97, 4.97, 4.97],
[ 4.98, 4.98, 4.98, ..., 4.98, 4.98, 4.98],
[ 4.99, 4.99, 4.99, ..., 4.99, 4.99, 4.99]])
# 现在,对该函数的求值运算就好办了,把这两个数组当做两个浮点数那样编写表达式即可:
In [78]: z = np.sqrt(xs **2 + ys ** 2)
In [79]: z
Out[79]:
array([[7.07106781, 7.06400028, 7.05693985, ..., 7.04988652, 7.0569398
5,
7.06400028],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.0498581
5,
7.05692568],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.0427835
4,
7.04985815],
...,
[7.04988652, 7.04279774, 7.03571603, ..., 7.0286414 , 7.0357160
3,
7.04279774],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.0427835
4,
7.04985815],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.0498581
5,
7.05692568]])
作为第 9 章的先导,我用 matplotlib 创建了这个二维数组的可视化:
In [160]: import matplotlib.pyplot as plt
In [161]: plt.imshow(z, cmap=plt.cm.gray); plt.colorbar()
Out[161]: <matplotlib.colorbar.Colorbar at 0x7f715e3fa630>
In [162]: plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")
Out[162]: <matplotlib.text.Text at 0x7f715d2de748>
见图 4-3。这张图是用matplotlib
的imshow
函数创建的。
将条件逻辑表述为数组运算
numpy.where
函数是三元表达式x if condition else y
的矢量化版本。假设我们有一个布尔数组和两个值数组:
In [165]: xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
In [166]: yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
In [167]: cond = np.array([True, False, True, True, False])
假设我们想要根据 cond
中的值选取xarr
和yarr
的值:当cond
中的值为True
时,选取xarr
的值,否则从yarr
中选取。列表推导式的写法应该如下所示:
In [168]: result = [(x if c else y)
.....: for x, y, c in zip(xarr, yarr, cond)]
In [169]: result
Out[169]: [1.1000000000000001, 2.2000000000000002, 1.3, 1.3999999999999999, 2.5]
这有几个问题。
- 第一,它对大数组的处理速度不是很快(因为所有工作都是由纯Python 完成的)。
- 第二,无法用于多维数组。若使用
np.where
,则可以将该功能写得非常简洁:
In [87]: result = np.where(cond, xarr, yarr)
In [88]: result
Out[88]: array([1.1, 2.2, 1.3, 1.4, 2.5])
np.where
的第二个和第三个参数不必是数组,它们都可以是标量值。在数据分析工作中,where
通常用于根据另一个数组而产生一个新的数组。假设有一个由随机数据组成的矩阵,你希望将所有正值替换为 2,将所有负值替换为-2。若利用np.where
,则会非常简单:
In [89]: arr = np.random.randn(4,4)
In [90]: arr
Out[90]:
array([[-0.88280159, 1.08426248, -0.6446468 , 0.14938064],
[-1.52576928, -0.69221969, -0.04114702, 0.35673447],
[-0.06152351, 1.39585055, -0.08607346, 0.35066278],
[-0.23933849, -2.13141676, -0.88463398, -1.07845589]])
In [91]: arr > 0
Out[91]:
array([[False, True, False, True],
[False, False, False, True],
[False, True, False, True],
[False, False, False, False]])
In [92]: np.where(arr >0,2,-2)
Out[92]:
array([[-2, 2, -2, 2],
[-2, -2, -2, 2],
[-2, 2, -2, 2],
[-2, -2, -2, -2]])
使用np.where
,可以将标量和数组结合起来。例如,我可用常数2替换arr中所有正的值:
In [93]: np.where(arr > 0, 2, arr)
Out[93]:
array([[-0.88280159, 2. , -0.6446468 , 2. ],
[-1.52576928, -0.69221969, -0.04114702, 2. ],
[-0.06152351, 2. , -0.08607346, 2. ],
[-0.23933849, -2.13141676, -0.88463398, -1.07845589]])
数学和统计方法
可以通过数组上的一组数学函数对整个数组或某个轴向的数据进行统计计算。sum、mean以及标准差std等聚合计算(aggregation,通常叫做约简(reduction))既可以当做数组的实例方法调用,也可以当做顶级NumPy函数使用。
这里,我生成了一些正态分布随机数据,然后做了聚类统计:
In [177]: arr = np.random.randn(5, 4)
In [178]: arr
Out[178]:
array([[ 2.1695, -0.1149, 2.0037, 0.0296],
[ 0.7953, 0.1181, -0.7485, 0.585 ],
[ 0.1527, -1.5657, -0.5625, -0.0327],
[-0.929 , -0.4826, -0.0363, 1.0954],
[ 0.9809, -0.5895, 1.5817, -0.5287]])
In [179]: arr.mean()
Out[179]: 0.19607051119998253
In [180]: np.mean(arr)
Out[180]: 0.19607051119998253
In [181]: arr.sum()
Out[181]: 3.9214102239996507
mean
和 sum
这类的函数可以接受一个axis
选项参数,用于计算该轴向上的统计值,最终结果是一个少一维的数组:
In [182]: arr.mean(axis=1)
Out[182]: array([ 1.022 , 0.1875, -0.502 , -0.0881, 0.3611])
In [183]: arr.sum(axis=0)
Out[183]: array([ 3.1693, -2.6345, 2.2381, 1.1486])
这里,arr.mean(1)
是“计算行的平均值”,arr.sum(0)
是“计算每列的和”。
其他如cumsum
和cumprod
之类的方法则不聚合,而是产生一个由中间结果组成的数组:
In [184]: arr = np.array([0, 1, 2, 3, 4, 5, 6, 7])
In [185]: arr.cumsum() # 新生成的数组的每个位置的数是原来的数字和的累加
Out[185]: array([ 0, 1, 3, 6, 10, 15, 21, 28])
在多维数组中,累加函数(如 cumsum)返回的是同样大小的数组,但是会根据每个低维的切片沿着标记轴计算部分聚类:
In [94]: arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [95]: arr
Out[95]:
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [96]: arr.cumsum(axis=0) # 纵向元素的累计和
Out[96]:
array([[ 0, 1, 2],
[ 3, 5, 7],
[ 9, 12, 15]], dtype=int32)
In [97]: arr.cumprod(axis=1) # 横轴元素的累计积
Out[97]:
array([[ 0, 0, 0],
[ 3, 12, 60],
[ 6, 42, 336]], dtype=int32)
表4-5列出了全部的基本数组统计方法。后续章节中有很多例子都会用到这些方法。
用于布尔型数组的方法
在上面这些方法中,布尔值会被强制转换为 1(True)和 0(False)。因此,sum 经常被用来对布尔型数组中的 True 值计数:
In [98]: arr = np.random.randn(100)
In [99]: (arr > 0).sum()
Out[99]: 57
另外还有两个方法any
和all
,它们对布尔型数组非常有用。any
用于测试数组中是否存在一个或多个True
,而all
则检查数组中所有值是否都是True
:
In [100]: bools = np.array([False, False, True, False])
In [101]: bools.any()
Out[101]: True
In [102]: bools.all()
Out[102]: False
这两个方法也能用于非布尔型数组,所有非 0 元素将会被当做 True。
排序
- 跟 Python 内置的列表类型一样,NumPy 数组也可以通过 sort 方法就地排序:
In [195]: arr = np.random.randn(6)
In [196]: arr
Out[196]: array([ 0.6095, -0.4938, 1.24 , -0.1357, 1.43 , -0.8469])
In [197]: arr.sort()
In [198]: arr
Out[198]: array([-0.8469, -0.4938, -0.1357, 0.6095, 1.24 , 1.43 ])
- 多维数组可以在任何一个轴向上进行排序,只需将轴编号传给 sort 即可:
In [199]: arr = np.random.randn(5, 3)
In [200]: arr
Out[200]:
array([[ 0.6033, 1.2636, -0.2555],
[-0.4457, 0.4684, -0.9616],
[-1.8245, 0.6254, 1.0229],
[ 1.1074, 0.0909, -0.3501],
[ 0.218 , -0.8948, -1.7415]])
In [201]: arr.sort(1)
In [202]: arr
Out[202]:
array([[-0.2555, 0.6033, 1.2636],
[-0.9616, -0.4457, 0.4684],
[-1.8245, 0.6254, 1.0229],
[-0.3501, 0.0909, 1.1074],
[-1.7415, -0.8948, 0.218 ]])
- 顶级方法
np.sort
返回的是数组的已排序副本,而就地排序则会修改数组本身。计算数组分位数最简单的办法是对其进行排序,然后选取特定位置的值:
In [105]: large_arr = np.random.randn(1000)
In [106]: large_arr.sort()
In [107]: large_arr[int(0.05 * len(large_arr))]
Out[107]: -1.7052992051341627
更多关于 NumPy 排序方法以及诸如间接排序之类的高级技术,请参阅附录A。在pandas 中还可以找到一些其他跟排序有关的数据操作(比如根据一列或多列对表格型数据进行排序)。
唯一化以及其它的集合逻辑
NumPy 提供了一些针对一维 ndarray 的基本集合运算。最常用的可能要数np.unique
了,它用于找出数组中的唯一值并返回已排序的结果:
In [206]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
In [207]: np.unique(names)
Out[207]:
array(['Bob', 'Joe', 'Will'],
dtype='<U4')
In [208]: ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4])
In [209]: np.unique(ints)
Out[209]: array([1, 2, 3, 4])
- 拿跟
np.unique
等价的纯 Python 代码来对比一下:
In [210]: sorted(set(names))
Out[210]: ['Bob', 'Joe', 'Will']
- 另一个函数
np.in1d
用于测试一个数组中的值在另一个数组中的成员资格,返回一个布尔型数组:
In [110]: values = np.array([6, 0, 0, 3, 2, 5, 6])
In [111]: np.in1d(values,[2,3,6])
Out[111]: array([ True, False, False, True, True, False, True])
4.4 用于数组的文件输入输出
NumPy 能够读写磁盘上的文本数据或二进制数据。这一小节只讨论 NumPy 的内置二进制格式,因为更多的用户会使用 pandas 或其它工具加载文本或表格数据(见第6章)。
np.save
和np.load
是读写磁盘数组数据的两个主要函数。默认情况下,数组是以未压缩的原始二进制格式保存在扩展名为.npy
的文件中的:
In [117]: arr = np.arange(10)
In [118]: np.save('some_array', arr)
# 如果文件路径末尾没有扩展名 .npy,则该扩展名会被自动加上。
# 然后就可以通过 np.load 读取磁盘上的数组:
In [119]: np.load('some_array.npy')
Out[119]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# 通过 np.savez 可以将多个数组保存到一个未压缩文件中,将数组以关键字参数的形式传入即可:
In [120]: np.savez('array_archive.npz', a=arr, b=arr)
# 加载 .npz 文件时,你会得到一个类似字典的对象,该对象会对各个数组进行延迟加载:
In [121]: arch = np.load('array_archive.npz')
In [122]: arch['b']
Out[122]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# 如果数据压缩的很好,就可以使用 numpy.savez_compressed:
In [123]: np.savez_compressed('arrays_compressed.npz', a=arr, b=arr)
4.5 线性代数
线性代数(如矩阵乘法、矩阵分解、行列式以及其他方阵数学等)是任何数组库的重要组成部分。不像某些语言(如 MATLAB),通过*
对两个二维数组相乘得到的是一个元素级的积,而不是一个矩阵点积。因此,NumPy 提供了一个用于矩阵乘法的dot
函数(既是一个数组方法也是numpy
命名空间中的一个函数):
In [223]: x = np.array([[1., 2., 3.], [4., 5., 6.]])
In [224]: y = np.array([[6., 23.], [-1, 7], [8, 9]])
In [225]: x
Out[225]:
array([[ 1., 2., 3.],
[ 4., 5., 6.]])
In [226]: y
Out[226]:
array([[ 6., 23.],
[ -1., 7.],
[ 8., 9.]])
In [227]: x.dot(y)
Out[227]:
array([[ 28., 64.],
[ 67., 181.]])
# x.dot(y)等价于np.dot(x, y):
In [228]: np.dot(x, y)
Out[228]:
array([[ 28., 64.],
[ 67., 181.]])
# 一个二维数组跟一个大小合适的一维数组的矩阵点积运算之后将会得到一个一维数组:
In [229]: np.dot(x, np.ones(3))
Out[229]: array([ 6., 15.])
# @符(类似Python 3.5)也可以用作中缀运算符,进行矩阵乘法:
In [230]: x @ np.ones(3)
Out[230]: array([ 6., 15.])
numpy.linalg
中有一组标准的矩阵分解运算以及诸如求逆和行列式之类的东西。它们跟 MATLAB 和 R 等语言所使用的是相同的行业标准线性代数库,如 BLAS、LAPACK、Intel MKL(Math Kernel Library,可能有,取决于你的 NumPy 版本)等:
In [124]: from numpy.linalg import inv, qr
In [125]: x = np.random.randn(5,5)
In [126]: mat = x.T.dot(x)
In [127]: inv(mat) # 计算方阵的逆
Out[127]:
array([[ 1.11432402, 2.38334321, 0.52453223, -1.35311044, 2.6920607
9],
[ 2.38334321, 6.94750325, 1.37627684, -3.76811382, 8.8297655
],
[ 0.52453223, 1.37627684, 0.56841399, -0.78134949, 1.3267284
2],
[-1.35311044, -3.76811382, -0.78134949, 2.13934199, -4.8158717
8],
[ 2.69206079, 8.8297655 , 1.32672842, -4.81587178, 13.6137554
1]])
In [128]: mat.dot(inv(mat))
Out[128]:
array([[ 1.00000000e+00, 1.77635684e-15, 2.22044605e-16,
0.00000000e+00, 0.00000000e+00],
[-3.10862447e-15, 1.00000000e+00, -1.77635684e-15,
1.77635684e-15, -5.32907052e-15],
[ 0.00000000e+00, -3.55271368e-15, 1.00000000e+00,
0.00000000e+00, -3.55271368e-15],
[-3.55271368e-15, 7.10542736e-15, -1.77635684e-15,
1.00000000e+00, -3.55271368e-15],
[ 2.22044605e-16, -8.88178420e-16, -1.11022302e-16,
0.00000000e+00, 1.00000000e+00]])
In [129]: q,r = qr(mat) # 计算 QR 分解
In [130]: r
Out[130]:
array([[ -6.41269138, -0.85145237, -3.14028915, -10.84629176,
-1.72071891],
[ 0. , -6.22211764, 0.5383217 , -10.13224216,
0.41002201],
[ 0. , 0. , -5.01465665, -4.20806697,
-1.0074459 ],
[ 0. , 0. , 0. , -1.32253775,
-0.50939452],
[ 0. , 0. , 0. , 0. ,
0.05817306]])
表达式X.T.dot(X)
计算 X
和它的转置X.T
的点积。
表4-7 中列出了一些最常用的线性代数函数。
4.6 伪随机数生成
numpy.random
模块对 Python 内置的 random 进行了补充,增加了一些用于高效生成多种概率分布的样本值的函数。例如,你可以用 normal 来得到一个标准正态分布的 4×4 样本数组:
In [131]: samples = np.random.normal(size=(4,4))
In [132]: samples
Out[132]:
array([[-0.26441984, 0.92011017, 1.43525863, 1.57755484],
[ 0.45573755, -0.02369316, 0.74781731, 0.33000411],
[-0.0575799 , 0.7124306 , -0.20261719, 0.09799322],
[-1.08571569, 0.22937284, 0.07188612, -1.97462298]])
而 Python 内置的 random 模块则只能一次生成一个样本值。从下面的测试结果中可以看出,如果需要产生大量样本值,numpy.random 快了不止一个数量级:
In [133]: from random import normalvariate
In [134]: N = 1000000
In [135]: %timeit samples = [normalvariate(0,1) for _ in range(N)]
955 ms ± 15.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [136]: %timeit np.random.normal(size=N)
36 ms ± 108 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
我们说这些都是伪随机数,是因为它们都是通过算法基于随机数生成器种子,在确定性的条件下生成的。你可以用 NumPy 的 np.random.seed
更改随机数生成种子:
In [244]: np.random.seed(1234)
numpy.random
的数据生成函数使用了全局的随机种子。要避免全局状态,你可以使用 numpy.random.RandomState
,创建一个与其它隔离的随机数生成器:
In [245]: rng = np.random.RandomState(1234)
In [246]: rng.randn(10)
Out[246]:
array([ 0.4714, -1.191 , 1.4327, -0.3127, -0.7206, 0.8872, 0.8596,
-0.6365, 0.0157, -2.2427])
表4-8列出了 numpy.random
中的部分函数。在下一节中,我将给出一些利用这些函数一次性生成大量样本值的范例
4.7 示例:随机漫步
我们通过模拟随机漫步来说明如何运用数组运算。先来看一个简单的随机漫步的例子:从 0 开始,步长 1 和 -1 出现的概率相等。
下面是一个通过内置的 random 模块以纯 Python 的方式实现 1000 步的随机漫步:
In [137]: import random
In [138]: position = 0
In [139]: walk = [position]
In [140]: steps = 1000
In [141]: for i in range(steps):
...: step = 1 if random.randint(0,1) else -1
...: position += step
...: walk.append(position)
图4-4是根据前100个随机漫步值生成的折线图:
In [249]: plt.plot(walk[:100])
不难看出,这其实就是随机漫步中各步的累计和,可以用一个数组运算来实现。因此,我用np.random
模块一次性随机产生 1000 个“掷硬币”结果(即两个数中任选一个),将其分别设置为 1 或 -1,然后计算累计和:
In [142]: nsteps = 1000
In [143]: draws = np.random.randint(0,2,size=nsteps)
In [144]: steps = np.where(draws > 0 ,1,-1)
In [146]: walk = steps.cumsum()
# 有了这些数据之后,我们就可以沿着漫步路径做一些统计工作了,比如求取最大值和最小值:
In [148]: walk.min()
Out[148]: 0
In [149]: walk.max()
Out[149]: 55
现在来看一个复杂点的统计任务——首次穿越时间,即随机漫步过程中第一次到达某个特定值的时间。假设我们想要知道本次随机漫步需要多久才能距离初始 0点至少 10 步远(任一方向均可)。np.abs(walk)>= 10
可以得到一个布尔型数组,它表示的是距离是否达到或超过 10,而我们想要知道的是第一个 10 或-10的索引。可以用argmax
来解决这个问题,它返回的是该布尔型数组第一个最大值的索引(True 就是最大值):
In [257]: (np.abs(walk) >= 10).argmax()
Out[257]: 37
注意,这里使用 argmax
并不是很高效,因为它无论如何都会对数组进行完全扫描。在本例中,只要发现了一个 True,那我们就知道它是个最大值了。
一次模拟多个随机漫步
如果你希望模拟多个随机漫步过程(比如5000个),只需对上面的代码做一点点修改即可生成所有的随机漫步过程。只要给numpy.random
的函数传入一个二元元组就可以产生一个二维数组,然后我们就可以一次性计算 5000 个随机漫步过程(一行一个)的累计和了:
In [150]: nwalks = 5000
In [151]: nsteps = 1000
In [153]: draws = np.random.randint(0,2, size=(nwalks,nsteps))
In [154]: steps = np.where(draws >0 ,1,-1)
In [156]: walks = steps.cumsum(1)
In [157]: walks
Out[157]:
array([[ 1, 0, 1, ..., -8, -9, -10],
[ 1, 0, -1, ..., 20, 21, 20],
[ 1, 0, 1, ..., 0, -1, -2],
...,
[ -1, 0, -1, ..., 16, 17, 16],
[ 1, 2, 3, ..., 0, -1, -2],
[ -1, 0, 1, ..., -16, -17, -16]], dtype=int32)
现在,我们来计算所有随机漫步过程的最大值和最小值:
In [158]: walks.max()
Out[158]: 140
In [159]: walks.min()
Out[159]: -137
# 得到这些数据之后,我们来计算 30 或 -30 的最小穿越时间。
# 这里稍微复杂些,因为不是 5000 个过程都到达了 30。我们可以用 any 方法来对此进行检查:
n [158]: walks.max()
Out[158]: 140
In [159]: walks.min()
Out[159]: -137
In [160]: hints30 = (np.abs(walks) >= 30).any(1)
In [161]: hints30
Out[161]: array([ True, True, False, ..., True, False, False])
In [162]: hints30.sum()
Out[162]: 3430
# 然后我们利用这个布尔型数组选出那些穿越了30(绝对值)的随机漫步(行),并调用 argmax 在轴 1 上获取穿越时间:
In [163]: crossing_times = (np.abs(walks[hints30]) >= 30).argmax(1)
In [164]: crossing_times.mean()
Out[164]: 501.35043731778427
# 请尝试用其他分布方式得到漫步数据。只需使用不同的随机数生成函数即可,
# 如 normal用于生成指定均值和标准差的正态分布数据:
In [271]: steps = np.random.normal(loc=0, scale=0.25,
.....: size=(nwalks, nsteps))
4.8 结论
虽然本书剩下的章节大部分是用 pandas 规整数据,我们还是会用到相似的基于数组的计算。在附录A中,我们会深入挖掘 NumPy 的特点,进一步学习数组的技巧。