目录
NumPy(Numerical Python),大多数科学计算的包都是用NumPy的数组作为构建基础。
部分功能:
- ndarray,一个具有矢量算术运算和复杂广播能力的快速且节省空间的多维数组。
- 用于对整组数据进行快速运算的标准数学函数(无需编写循环)。
- 用于读写磁盘数据的工具以及用于操作内存映射文件的工具。
- 线性代数、随机数生成以及傅里叶变换功能。
- 用于集成由C、C++、Fortran等语言编写的代码的工具。
对于大部分数据分析应用而言,我最关注的功能主要集中在:
- 用于数据整理和清理、子集构造和过滤、转换等快速的矢量化数组运算。
- 常用的数组算法,如排序、唯一化、集合运算等。
- 高效的描述统计和数据聚合/摘要运算。
- 用于异构数据集的合并/连接运算的数据对齐和关系型数据运算。
- 将条件逻辑表述为数组表达式(而不是带有if-elif-else分支的循环)。
- 数据的分组运算(聚合、转换、函数应用等)。
NumPy能高效处理大数组的数据的原因:
- NumPy是在一个连续的内存块中存储数据,独立于其他Python内置对象。NumPy的C语言编写的算法库可以操作内存,而不必进行类型检查或其它前期工作。比起Python的内置序列,NumPy数组使用的内存更少。
-
可以在整个数组上执行复杂的计算,而不需要Python的for循环。
4.1 NumPy的ndarray:一种多维数组对象
N维数组对象(即ndarray)是一个快速而灵活的大数据集容器。用这种数组对整块数据执行数学运算语法跟标量元素之间的运算一样。
In [12]: import numpy as np
In [13]: data = np.random.randn(2, 3) In [14]: data
Out[14]: array([[-0.2047, 0.4789, -0.5194], [-0.5557, 1.9658, 1.3934]])
In [15]: data * 10
Out[15]: array([[ -2.0471, 4.7894, -5.1944], [ -5.5573, 19.6578, 13.9341]])
In [16]: data + data
Out[16]: array([[-0.4094, 0.9579, -1.0389], [-1.1115, 3.9316, 2.7868]])
ndarray是一个通用的同构数据多维容器,所有元素必须是相同类型。每个数组都有一个shape(一个表示各维度大小的元组)和一个dtype(一个用于说明数组数据类型的对象):
In [17]: data.shape
Out[17]: (2, 3)
In [18]: data.dtype
Out[18]: dtype('float64')
创建ndarray
使用array函数创建数组,接受一切序列型的对象,然后产生一个新的含有传入数据的NumPy数组
In [19]: data1 = [6, 7.5, 8, 0, 1]
In [20]: arr1 = np.array(data1)
In [21]: arr1
Out[21]: array([ 6. , 7.5, 8. , 0. , 1. ])
嵌套序列(比如由一组等长列表组成的列表)将会被转换为一个多维数组:
In [22]: data2 = [[1, 2, 3, 4], [5, 6, 7, 8]]
In [23]: arr2 = np.array(data2)
In [24]: arr2
Out[24]:
array([[1, 2, 3, 4], [5, 6, 7, 8]])
可以用属性ndim和shape验证:NumPy数组arr2的两个维度
In [25]: arr2.ndim
Out[25]: 2
In [26]: arr2.shape
Out[26]: (2, 4)
np.array会尝试为新建的这个数组推断出一个较为合适的数据类型。数据类型保存在一个特殊的dtype对象中。
In [27]: arr1.dtype
Out[27]: dtype('float64')
In [28]: arr2.dtype
Out[28]: dtype('int64')
其他函数也可以新建数组。zeros和ones分别可以创建指定长度或形状的全0或全1数组。empty可以创建一个没有任何具体值的数组。
In [29]: np.zeros(10)
Out[29]: array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [30]: np.zeros((3, 6))
Out[30]: array([[ 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0.]])
In [31]: np.empty((2, 3, 2))
Out[31]: array([[[ 0., 0.], [ 0., 0.], [ 0., 0.]], [[ 0., 0.], [ 0., 0.], [ 0., 0.]]])
arange是Python内置函数range的数组版:
In [32]: np.arange(15)
Out[32]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
表4-1列出了一些数组创建函数。如果没有特别指定,NumPy数据类型基本都是float64(浮点数)。
ndarray的数据类型
dtype(数据类型)含有ndarray将一块内存解释为特定数据类型所需的信息:
In [33]: arr1 = np.array([1, 2, 3], dtype=np.float64)
In [34]: arr2 = np.array([1, 2, 3], dtype=np.int32)
In [35]: arr1.dtype
Out[35]: dtype('float64')
In [36]: arr2.dtype
Out[36]: dtype('int32')
数值型dtype的命名方式相同:一个类型名(如float或int),后面跟一个用于表示各元素位长的数字。标准的双精度浮点值(即Python中的float对象)需要占用8字节(即64位)。因此,该类型在NumPy中就记作float64。表4-2列出了NumPy所支持的全部数据类型。
可以通过ndarray的astype方法明确地将一个数组从一个dtype转换成另一个dtype:
In [37]: arr = np.array([1, 2, 3, 4, 5])
In [38]: arr.dtype
Out[38]: dtype('int64')
In [39]: float_arr = arr.astype(np.float64)
In [40]: float_arr.dtype
Out[40]: dtype('float64')
笔记:调用astype总会创建一个新的数组(一个数据的备份),即使新的dtype与旧的dtype相同。
NumPy数组的运算
数组很重要,不用编写循环即可对数据执行批量运算。NumPy用户称其为矢量化(vectorization)。大小相等的数组之间的任何算术运算都会将运算应用到元素级:
数组与标量的算术运算会将标量值传播到各个元素:
大小相同的数组之间的比较会生成布尔值数组:
不同大小的数组之间的运算叫做广播(broadcasting)。
基本的索引和切片
将一个标量值赋值给一个切片时(如arr[5:8]=12),该值会自动传播(“广播”)到整个选区。跟列表的区别:数组切片是原始数组的视图。这意味着数据不会被复制,视图上的任何修改都会直接反映到源数组上。
图4-1说明了二维数组的索引方式。轴0作为行,轴1作为列。
在多维数组中,如果省略了后面的索引,则返回对象会是一个维度低一点的ndarray
切片索引
ndarray的切片语法跟Python列表这样的一维对象差不多:
对于之前的二维数组arr2d,其切片方式稍显不同:
切片是沿着一个轴向选取元素的。表达式arr2d[:2]可以被认为是“选取arr2d的前两行”。
切片从0开始包括起始元素,不包含结束元素。start或stop都可以被省略,省略之后,分别默认序列的开头和结尾:因此,结果中包含的元素个数是stop - start。所以arr2d[:2]里是指0,1两行。之所以是行,默认第一个维度是行即第0轴。加了冒号后就指0轴和1轴两个维度了。
可以一次传入多个切片,就像传入多个索引那样:
arr2d[:2, 1:]即表示前两行的第2列以后的数据。(包含第二列,因为是起始元素,第因为是从0开始,所以是第二列)
arr2d[1, :2]即第二行的前两列。
“只有冒号”表示选取整个轴
arr2d[:, :1]即所有行第一个列。从第0列开始,不包含第一列。
对切片表达式的赋值操作也会被扩散到整个选区:
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 [98]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
In [99]: data = np.random.randn(7, 4) In
[100]: names
Out[100]: array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='<U4')
In [101]: data
Out[101]: array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.0072, -1.2962, 0.275 , 0.2289], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 1.669 , -0.4386, -0.5397, 0.477 ], [ 3.2489, -1.0212, -0.5771, 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [-0.7135, -0.8312, -2.3702, -1.8608]])
假设每个名字都对应data数组中的一行,而我们想要选出对应于名字”Bob”的所有行。跟算术运算一样,数组的比较运算(如==)也是矢量化的。因此,对names和字符串”Bob”的比较运算将会产生一个布尔型数组:
In [102]: names == 'Bob'
Out[102]: array([ True, False, False, True, False, False, False], dtype=bool)
这个布尔型数组可用于数组索引:
In [103]: data[names == 'Bob']
Out[103]:
array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.669 , -0.4386, -0.5397, 0.477 ]])
布尔型数组长度必须跟被索引的轴长度一致。
可以将布尔数组跟切片、整数混合使用:
In [104]: data[names == 'Bob', 2:]
Out[104]: array([[ 0.769 , 1.2464], [-0.5397, 0.477 ]]) In [105]: data[names == 'Bob', 3]
Out[105]: array([ 1.2464, 0.477 ])
要选择除了“Bob”以外的其他值,既可以使用不等于符号(!=),也可以通过~对条件进行否定
In [106]: names != 'Bob'
Out[106]: array([False, True, True, False, True, True, True], dtype=bool)
In [107]: data[~(names == 'Bob')] Out[107]:
array([[ 1.0072, -1.2962, 0.275 , 0.2289], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 3.2489, -1.0212, -0.5771, 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [-0.7135, -0.8312, -2.3702, -1.8608]])
选取这三个名字中的两个需要组合应用多个布尔条件,使用&(和)、|(或)之类的布尔算术运算符即可:
In [110]: mask = (names == 'Bob') | (names == 'Will')
In [111]: mask Out[111]: array([ True, False, True, True, True, False, False], dtype=bool)
In [112]: data[mask]
Out[112]:
array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.3529, 0.8864, -2.0016, -0.3718], [ 1.669 , -0.4386, -0.5397, 0.477 ], [ 3.2489, -1.0212, -0.5771, 0.1241]])
通过布尔型索引选取数组中的数据,将总是创建数据的副本,即使返回一模一样的数组也是如此。
注意:Python关键字and和or在布尔型数组中无效。要是用&与|。
通过布尔型数组设置值。例如将data中的所有负值都设置为0,只需:
In [113]: data[data < 0] = 0
In [114]: data
Out[114]:
array([[ 0.0929, 0.2817, 0.769 , 1.2464], [ 1.0072, 0. , 0.275 , 0.2289], [ 1.3529, 0.8864, 0. , 0. ], [ 1.669 , 0. , 0. , 0.477 ], [ 3.2489, 0. , 0. , 0.1241], [ 0.3026, 0.5238, 0.0009, 1.3438], [ 0. , 0. , 0. , 0. ]])
数组转置和轴对换
转置是重塑的一种特殊形式,它返回的是源数据的视图(不会进行任何复制操作)。数组不仅有transpose方法,还有一个特殊的T属性:
In [126]: arr = np.arange(15).reshape((3, 5))
In [127]: arr
Out[127]:
array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14]])
In [128]: arr.T
Out[128]:
array([[ 0, 5, 10], [ 1, 6, 11], [ 2, 7, 12], [ 3, 8, 13], [ 4, 9, 14]])
在进行矩阵计算时,经常需要用到该操作,比如利用np.dot计算矩阵内积:
In [129]: arr = np.random.randn(6, 3)
In [130]: arr
Out[130]:
array([[-0.8608, 0.5601, -1.2659], [ 0.1198, -1.0635, 0.3329], [-2.3594, -0.1995, -1.542 ], [-0.9707, -1.307 , 0.2863], [ 0.378 , -0.7539, 0.3313], [ 1.3497, 0.0699, 0.2467]])
In [131]: np.dot(arr.T, arr) #3X6 dot 6X3 = 3X3
Out[131]:
array([[ 9.2291, 0.9394, 4.948 ], [ 0.9394, 3.7662, -1.3622], [ 4.948 , -1.3622, 4.3437]])
对于高维数组,transpose需要得到一个由轴编号组成的元组才能对这些轴进行转置:
In [132]: arr = np.arange(16).reshape((2, 2, 4))
In [133]: arr
Out[133]:
array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]], [[ 8, 9, 10, 11], [12, 13, 14, 15]]])
In [134]: arr.transpose((1, 0, 2)) Out[134]:
array([[[ 0, 1, 2, 3], [ 8, 9, 10, 11]], [[ 4, 5, 6, 7], [12, 13, 14, 15]]])
这里,第一个轴被换成了第二个,第二个轴被换成了第一个,最后一个轴不变。
例如可以只看前2维,则可将每一行当成一个数压缩成2x2的[[ 0], [ 4]], [[ 8], [12]]即第一行为0,8,第二行为4,12换轴后则为[[ 0], [ 8]], [[ 4], [12]]再拓展开即结果。
In [135]: arr
Out[135]:
array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]],
[[ 8, 9, 10, 11], [12, 13, 14, 15]]])
In [136]: arr.swapaxes(1, 2)
Out[136]:
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 [21]: arr=np.arange(10)
In [22]: np.sqrt(arr)
Out[22]:
array([ 0. , 1. , 1.41421356, 1.73205081, 2. ,
2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
In [23]: np.exp(arr)
Out[23]:
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 [25]: x=np.random.randn(8)
In [26]: y=np.random.randn(8)
In [27]: x
Out[27]:
array([-1.47208377, 1.77425071, 1.60968406, 2.51201825, -0.42125923,
-0.1886985 , 0.22943125, -1.33311273])
In [28]: y
Out[28]:
array([-0.79881692, -0.79049443, 0.4106048 , -0.65989715, 1.19432883,
2.66407523, -1.08633441, -0.86861555])
In [30]: np.maximum(x,y) #元素级的最大值
Out[30]:
array([-0.79881692, 1.77425071, 1.60968406, 2.51201825, 1.19432883,
2.66407523, 0.22943125, -0.86861555])
这里,numpy.maximum计算了x和y中元素级别最大的元素。
有些ufunc可以返回多个数组。例如modf,它是Python内置函数divmod的矢量化版本,它会返回浮点数数组的小数和整数部分:
In [31]: arr=np.random.randn(7)*5
In [32]: np.modf(arr)
Out[32]:
(array([-0.6184109 , -0.69467218, 0.61799246, 0.51558266, -0.26625316,
0.25942319, 0.15090508]),
array([ -5., -0., 1., 3., -6., 11., 6.]))
Ufuncs接受out选项参数,可以让它们在数组的原地进行操作:
In [151]: arr
Out[151]:
array([-3.2623, -6.0915, -6.663 , 5.3731, 3.6182, 3.45 , 5.0077])
In [152]: np.sqrt(arr)
Out[152]:
array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378])
In [153]: np.sqrt(arr, arr)
Out[153]:
array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378]) In [154]: arr Out[154]: array([ nan, nan, nan, 2.318 , 1.9022, 1.8574, 2.2378])
4.3 利用数组进行数据处理
numpy数组使你可以将许多数据处理表达任务表述为简洁的数组表达式。(否则需要编写循环)。
用数组表达式代替循环的做法,通常被称之为 矢量化。
一般来说,矢量化的数组要比等价的python方式上快的多,尤其是各种数值的计算。
假如,我们想在一组值(网格型)上计算sqrt(x^2+y^2) 。
np.meshgrid函数接受两个一维数组,并产生两个二维矩阵(对应于两个数组中所有的(x,y)对)
In [155]: points = np.arange(-5, 5, 0.01) # 1000 equally spaced points
In [156]: xs, ys = np.meshgrid(points, points)
In [157]: ys
Out[157]:
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 [158]: z = np.sqrt(xs ** 2 + ys ** 2)
In [159]: z
Out[159]:
array([[ 7.0711, 7.064 , 7.0569, ..., 7.0499, 7.0569, 7.064 ],
[ 7.064 , 7.0569, 7.0499, ..., 7.0428, 7.0499, 7.0569],
[ 7.0569, 7.0499, 7.0428, ..., 7.0357, 7.0428, 7.0499],
...,
[ 7.0499, 7.0428, 7.0357, ..., 7.0286, 7.0357, 7.0428],
[ 7.0569, 7.0499, 7.0428, ..., 7.0357, 7.0428, 7.0499],
[ 7.064 , 7.0569, 7.0499, ..., 7.0428, 7.0499, 7.0569]])
将条件逻辑表述为数组运算
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]
问题:1.它对大数组的处理速度不是很快2.无法用于多维数组。
若使用np.where,则可以将该功能写得非常简洁
In [170]: result = np.where(cond, xarr, yarr)
In [171]: result
Out[171]: array([ 1.1, 2.2, 1.3, 1.4, 2.5])
np.where的第二个和第三个参数不必是数组,可以是标量值。在数据分析工作中,where通常用于根据另一个数组而产生一个新的数组。假设有一个由随机数据组成的矩阵,你希望将所有正值替换为2,将所有负值替换为-2。若利用np.where,则会非常简单:
In [172]: arr = np.random.randn(4, 4)
In [173]: arr
Out[173]:
array([[-0.5031, -0.6223, -0.9212, -0.7262], [ 0.2229, 0.0513, -1.1577, 0.8167], [ 0.4336, 1.0107, 1.8249, -0.9975], [ 0.8506, -0.1316, 0.9124, 0.1882]])
In [174]: arr > 0
Out[174]:
array([[False, False, False, False], [ True, True, False, True], [ True, True, True, False], [ True, False, True, True]], dtype=bool)
In [175]: np.where(arr > 0, 2, -2)
Out[175]:
array([[-2, -2, -2, -2], [ 2, 2, -2, 2], [ 2, 2, 2, -2], [ 2, -2, 2, 2]])
np.where可以将标量和数组结合起来。例如,用常数2替换arr中所有正的值:
In [176]: np.where(arr > 0, 2, arr) # set only positive values to 2
Out[176]:
array([[-0.5031, -0.6223, -0.9212, -0.7262], [ 2. , 2. , -1.1577, 2. ], [ 2. , 2. , 2. , -0.9975], [ 2. , -0.1316, 2. , 2. ]])
传递给where的数组大小可以不相等,甚至可以是标量值。
数学和统计方法
可以通过数组上的一组数学函数对整个数组或某个轴向的数据进行统计计算。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)#把1轴即x轴上每一行的数据压缩成1个数
Out[182]: array([ 1.022 , 0.1875, -0.502 , -0.0881, 0.3611])
In [183]: arr.sum(axis=0)#把1轴即y轴上每一列的数据压缩成1个数
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 [186]: arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [187]: arr
Out[187]:
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [188]: arr.cumsum(axis=0)
Out[188]:
array([[ 0, 1, 2], [ 3, 5, 7], [ 9, 12, 15]])
In [189]: arr.cumprod(axis=1)
Out[189]:
array([[ 0, 0, 0], [ 3, 12, 60], [ 6, 42, 336]])
表4-5列出了全部的基本数组统计方法。后续章节中有很多例子都会用到这些方法。
用于布尔型数组的方法
在上面这些方法中,布尔值会被强制转换为1(True)和0(False)。因此,sum经常被用来对布尔型数组中的True值计数:
In [190]: arr = np.random.randn(100)
In [191]: (arr > 0).sum() # Number of positive values
Out[191]: 42
any用于测试数组中是否存在一个或多个True,而all则检查数组中所有值是否都是True:
In [192]: bools = np.array([False, False, True, False])
In [193]: bools.any()
Out[193]: True
In [194]: bools.all()
Out[194]: 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 [203]: large_arr = np.random.randn(1000)
In [204]: large_arr.sort()
In [205]: large_arr[int(0.05 * len(large_arr))] # 5% quantile
Out[205]: -1.5311513550102103
唯一化以及其它的集合逻辑
针对一维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 [211]: values = np.array([6, 0, 0, 3, 2, 5, 6])
In [212]: np.in1d(values, [2, 3, 6])
Out[212]: array([ True, False, False, True, True, False, True], dtype=bool)
NumPy中的集合函数请参见表4-6。
4.4 用于数组的文件输入输出
NumPy能够读写磁盘上的文本数据或二进制数据。
np.save和np.load是读写磁盘数组数据的两个主要函数。默认情况下,数组是以未压缩的原始二进制格式保存在扩展名为.npy的文件中的:
In [213]: arr = np.arange(10)
In [214]: np.save('some_array', arr)
如果文件路径末尾没有扩展名.npy,则该扩展名会被自动加上。然后就可以通过np.load读取磁盘上的数组:
In [215]: np.load('some_array.npy')
Out[215]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
np.savez可以将多个数组保存到一个未压缩文件中,将数组以关键字参数的形式传入即可:
In [216]: np.savez('array_archive.npz', a=arr, b=arr)
加载.npz文件时,你会得到一个类似字典的对象,该对象会对各个数组进行延迟加载:
In [217]: arch = np.load('array_archive.npz')
In [218]: arch['b']
Out[218]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
如果数据压缩的很好,就可以使用numpy.savez_compressed:
In [219]: 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 [231]: from numpy.linalg import inv, qr
In [232]: X = np.random.randn(5, 5)
In [233]: mat = X.T.dot(X)
In [234]: inv(mat)
Out[234]:
array([[ 933.1189, 871.8258, -1417.6902, -1460.4005, 1782.1391], [ 871.8258, 815.3929, -1325.9965, -1365.9242, 1666.9347], [-1417.6902, -1325.9965, 2158.4424, 2222.0191, -2711.6822], [-1460.4005, -1365.9242, 2222.0191, 2289.0575, -2793.422 ], [ 1782.1391, 1666.9347, -2711.6822, -2793.422 , 3409.5128]])
In [235]: mat.dot(inv(mat))
Out[235]:
array([[ 1., 0., -0., -0., -0.], [-0., 1., 0., 0., 0.], [ 0., 0., 1., 0., 0.], [-0., 0., 0., 1., -0.], [-0., 0., 0., 0., 1.]])
In [236]: q, r = qr(mat)
In [237]: r
Out[237]:
array([[-1.6914, 4.38 , 0.1757, 0.4075, -0.7838], [ 0. , -2.6436, 0.1939, -3.072 , -1.0702], [ 0. , 0. , -0.8138, 1.5414, 0.6155], [ 0. , 0. , 0. , -2.6445, -2.1669], [ 0. , 0. , 0. , 0. , 0.0002]])
表达式X.T.dot(X)计算X和它的转置X.T的点积。
表4-7中列出了一些最常用的线性代数函数。
4.6 伪随机数生成
numpy.random模块对Python内置的random进行了补充,增加了一些用于高效生成多种概率分布的样本值的函数。例如,你可以用normal来得到一个标准正态分布的4×4样本数组:
In [238]: samples = np.random.normal(size=(4, 4))
In [239]: samples
Out[239]:
array([[ 0.5732, 0.1933, 0.4429, 1.2796],
[ 0.575 , 0.4339, -0.7658, -1.237 ],
[-0.5367, 1.8545, -0.92 , -0.1082],
[ 0.1525, 0.9435, -1.0953, -0.144 ]])
而Python内置的random模块则只能一次生成一个样本值。从下面的测试结果中可以看出,如果需要产生大量样本值,numpy.random快了不止一个数量级:
In [240]: from random import normalvariate
In [241]: N = 1000000
In [242]: %timeit samples = [normalvariate(0, 1) for _ in range(N)]
1.77 s +- 126 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)
In [243]: %timeit np.random.normal(size=N)
61.7 ms +- 1.32 ms 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 [247]: import random
.....: position = 0
.....: walk = [position]
.....: steps = 1000
.....: 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 [251]: nsteps = 1000
In [252]: draws = np.random.randint(0, 2, size=nsteps)
In [253]: steps = np.where(draws > 0, 1, -1)
In [254]: walk = steps.cumsum()
有了这些数据之后,我们就可以沿着漫步路径做一些统计工作了,比如求取最大值和最小值:
In [255]: walk.min()
Out[255]: -3
In [256]: walk.max()
Out[256]: 31
现在来看一个复杂点的统计任务——首次穿越时间,即随机漫步过程中第一次到达某个特定值的时间。假设我们想要知道本次随机漫步需要多久才能距离初始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 [258]: nwalks = 5000
In [259]: nsteps = 1000
In [260]: draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1
In [261]: steps = np.where(draws > 0, 1, -1)
In [262]: walks = steps.cumsum(1)
In [263]: walks
Out[263]:
array([[ 1, 0, 1, ..., 8, 7, 8],
[ 1, 0, -1, ..., 34, 33, 32],
[ 1, 0, -1, ..., 4, 5, 4],
...,
[ 1, 2, 1, ..., 24, 25, 26],
[ 1, 2, 3, ..., 14, 13, 14],
[ -1, -2, -3, ..., -24, -23, -22]])
现在,我们来计算所有随机漫步过程的最大值和最小值:
In [264]: walks.max()
Out[264]: 138
In [265]: walks.min()
Out[265]: -133
得到这些数据之后,我们来计算30或-30的最小穿越时间。这里稍微复杂些,因为不是5000个过程都到达了30。我们可以用any方法来对此进行检查:
In [266]: hits30 = (np.abs(walks) >= 30).any(1)
In [267]: hits30
Out[267]: array([False, True, False, ..., False, True, False], dtype=bool)
In [268]: hits30.sum() # Number that hit 30 or -30
Out[268]: 3410
然后我们利用这个布尔型数组选出那些穿越了30(绝对值)的随机漫步(行),并调用argmax在轴1上获取穿越时间:
In [269]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
In [270]: crossing_times.mean()
Out[270]: 498.88973607038122
请尝试用其他分布方式得到漫步数据。只需使用不同的随机数生成函数即可,如normal用于生成指定均值和标准差的正态分布数据:
In [271]: steps = np.random.normal(loc=0, scale=0.25,
.....: size=(nwalks, nsteps))
4.8 结论
虽然本书剩下的章节大部分是用pandas规整数据,我们还是会用到相似的基于数组的计算。在附录A中,我们会深入挖掘NumPy的特点,进一步学习数组的技巧。