利用Python进行数据分析学习笔记三:NumPy基础:数组和矢量计算

一、Numpy基础

NumPy之于数值计算特别重要的原因之⼀,是因为它可以⾼效处理⼤数组的数据。这是因为:
NumPy是在⼀个连续的内存块中存储数据,独⽴于其他Python内置对象。NumPy的C语⾔编写的算法库可以操作内存,⽽不必进⾏类型检查或其它前期⼯作。⽐起Python的内置序列,NumPy数组使⽤的内存更少。
NumPy可以在整个数组上执⾏复杂的计算,⽽不需要Python的for循环。

要搞明⽩具体的性能差距,考察⼀个包含⼀百万整数的数组,和⼀个等价的Python列表:

In [7]: import numpy as np
In [8]: my_arr = np.arange(1000000)
In [9]: my_list = list(range(1000000))

各个序列分别乘以2:

In [10]: %time for _ in range(10): my_arr2 = my_arr * 2
CPU times: user 20 ms, sys: 50 ms, total: 70 ms
Wall time: 72.4 ms
In [11]: %time for _ in range(10): my_list2 = [x * 2 for x in my_list]
CPU times: user 760 ms, sys: 290 ms, total: 1.05 s
Wall time: 1.05 s

基于NumPy的算法要⽐纯Python快10到100倍(甚⾄更快),并且使⽤的内存更少。

1、Numpy的ndarray:一种多维数组对象

NumPy最重要的⼀个特点就是其N维数组对象(即ndarray),该对象是⼀个快速⽽灵活的⼤数据集容器。你可以利⽤这种数组对整块数据执⾏⼀些数学运算,其语法跟标量元素之间的运算⼀样。

要明⽩Python是如何利⽤与标量值类似的语法进⾏批次计算,我先引⼊NumPy,然后⽣成⼀个包含随机数据的⼩数组:

In [12]: import numpy as np
# Generate some random data
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],
116
 [ -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]])

第⼀个例⼦中,所有的元素都乘以10。第⼆个例⼦中,每个元素都与⾃身相加。

笔记:在本章及全书中,我会使⽤标准的NumPy惯⽤法 import numpy as np 。你当然也可以在代码中使⽤ from numpy import * ,但不建议这么做。 numpy 的命名空间很⼤,包含许多函数,其中⼀些的名字与Python的内置函
数重名(⽐如min和max)。

darray是⼀个通⽤的同构数据多维容器,也就是说,其中的所有元素必须是相同类型的。每个数组都有⼀个shape(⼀个表示各维度⼤⼩的元组)和⼀个dtype(⼀个⽤于说明数组数据类型的对象):

In [17]: data.shape
Out[17]: (2, 3)
In [18]: data.dtype
Out[18]: dtype('float64')

当你在本书中看到“数组”、“NumPy数组”、"ndarray"时,基本上都指的是同⼀样东⻄,即ndarray对象。

2、创建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]])

因为data2是列表的列表,NumPy数组arr2的两个维度的shape是从data2引⼊的。可以⽤属性ndim和shape验证:

ndim指维度:两个维度

In [25]: arr2.ndim
Out[25]: 2
In [26]: arr2.shape
Out[26]: (2, 4)

(2,4)代表有两个维度,每个维度4个元素
在这里插入图片描述
除⾮特别说明(稍后将会详细介绍),np.array会尝试为新建的这个数组推断出⼀个较为合适的数据类型。数据类型保存在⼀个特殊的dtype对象中。⽐如说,在上⾯的两个例⼦中,我们有:

In [27]: arr1.dtype
Out[27]: dtype('float64')
In [28]: arr2.dtype
Out[28]: dtype('int64')

除np.array之外,还有⼀些函数也可以新建数组。⽐如,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)) #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.]]])

注意:认为np.empty会返回全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])

下表列出了⼀些数组创建函数。由于NumPy关注的是数值计算,因此,如果没有特别指定,数据类型基本都是float64(浮点数)。
在这里插入图片描述

3、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是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('int64')
In [39]: float_arr = arr.astype(np.float64)
In [40]: float_arr.dtype
Out[40]: dtype('float64')

在这里插入图片描述

本质上是生成了新的一个ndarray

在本例中,整数被转换成了浮点数。如果将浮点数转换成整数,则⼩数部分将会被截取删除:

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], dtype=int32)

如果某字符串数组表示的全是数字,也可以⽤astype将其转换为数值形式:

In [44]: numeric_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_)
In [45]: numeric_strings.astype(float)
Out[45]: array([ 1.25, -9.6 , 42. ])

注意:使⽤numpy.string_类型时,⼀定要⼩⼼,因为NumPy的字符串数据是⼤⼩固定的,发⽣截取时,不会发出警告。pandas提供了更多⾮数值数据的便利的处理⽅法。
就是numpy.string_,最后有一小横杠!

如果转换过程因为某种原因⽽失败了(⽐如某个不能被转换为float64的字符串),就会引发⼀个ValueError。
这⾥,我⽐较懒,写的是float⽽不是np.float64;NumPy很聪明,它会将Python类型映射到等价的dtype上。
数组的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相同。

4、NumPy数组的运算

数组很重要,因为它使你不⽤编写循环即可对数据执⾏批量运算。NumPy⽤户称其为⽮量化(vectorization)。 (很重要!)⼤⼩相等的数组之间的任何算术运算都会将运算应⽤到元素级:

In [51]: arr = np.array([[1., 2., 3.], [4., 5., 6.]])
In [52]: arr
Out[52]:
array([[ 1., 2., 3.],
[ 4., 5., 6.]])
In [53]: arr * arr
Out[53]:
array([[ 1., 4., 9.],
[ 16., 25., 36.]])
In [54]: arr - arr
Out[54]:
array([[ 0., 0., 0.],
[ 0., 0., 0.]])

数组与标量的算术运算会将标量值传播到各个元素:

In [55]: 1 / arr
Out[55]:
array([[ 1. , 0.5 , 0.3333],
[ 0.25 , 0.2 , 0.1667]])
In [56]: arr ** 0.5 #几次方用 ** 不用^
Out[56]:
array([[ 1. , 1.4142, 1.7321],
[ 2. , 2.2361, 2.4495]])

⼤⼩相同的数组之间的⽐较会⽣成布尔值数组:

In [57]: arr2 = np.array([[0., 4., 1.], [7., 2., 12.]])
In [58]: arr2
Out[58]:
array([[ 0., 4., 1.],
[ 7., 2., 12.]])
In [59]: arr2 > arr
Out[59]:
array([[False, True, False],
[ True, False, True]], dtype=bool)

(比较是一个个,每个元素进行比较!)

不同⼤⼩的数组之间的运算叫做⼴播(broadcasting),将在附录A中对其进⾏详细讨论。本书的内容不需要对⼴播机制有多深的理解。

5、基本的索引和切片

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),该值会⾃动传播(也就说后⾯将会讲到的“⼴播”)到整个选区。 跟列表最重要的区别在于,数组切⽚是原始数组的视图。这意味着数据不会被复制,视图上的任何修改都会直接反映到源数组上。

(对一个切片进行赋值时会改变原始数组中的数据,改变的哪个数字由切片的开始、结束序号决定(实际是开始 - 结束序号-1))

作为例⼦,先创建⼀个arr的切⽚:

In [66]: arr_slice = arr[5:8]
In [67]: arr_slice
Out[67]: array([12, 12, 12])

现在,当我修稿arr_slice中的值,变动也会体现在原始数组arr中:

In [68]: arr_slice[1] = 12345
In [69]: arr
Out[69]: array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8,9])

切⽚[ : ]会给数组中的所有值赋值:

In [70]: arr_slice[:] = 64
In [71]: arr
Out[71]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])

如果你刚开始接触NumPy,可能会对此感到惊讶(尤其是当你曾经⽤过其他热衷于复制数组数据的编程语⾔)。由于NumPy的设计⽬的是处理⼤数据,所以你可以想象⼀下,假如NumPy坚持要将数据复制来复制去的话会产⽣何等的性能和内存问题。

注意:如果你想要得到的是ndarray切⽚的⼀份副本⽽⾮视图,就需要明确地进⾏复制操作,例如 arr[5:8].copy() 。

对于⾼维度数组,能做的事情更多。在⼀个⼆维数组中,各索引位置上的元素不再是标量⽽是⼀维数组:

In [72]: arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
In [73]: arr2d[2]
Out[73]: array([7, 8, 9])

因此,可以对各个元素进⾏递归访问,但这样需要做的事情有点多。你可以传⼊⼀个以逗号隔开的索引列表来选取单个元素。也就是说,下⾯两种⽅式是等价的:

In [74]: arr2d[0][2]
Out[74]: 3
In [75]: arr2d[0, 2]
Out[75]: 3

在多维数组中,如果省略了后⾯的索引,则返回对象会是⼀个维度低⼀点的ndarray(它含有⾼⼀级维度上的所有数据)。因此,在2×2×3数组arr3d中:

In [76]: arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
In [77]: arr3d
Out[77]:
array([[[ 1, 2, 3],
[ 4, 5, 6]],
[[ 7, 8, 9],
[10, 11, 12]]])

arr3d[0]是⼀个2×3数组:

In [78]: arr3d[0]
Out[78]:
array([[1, 2, 3],
[4, 5, 6]])

标量值和数组都可以被赋值给arr3d[0]:
(赋标量值时,则会让数组中所有的元素值均为该值)

In [79]: old_values = arr3d[0].copy()
In [80]: arr3d[0] = 42
In [81]: arr3d
Out[81]:
array([[[42, 42, 42],
[42, 42, 42]],
[[ 7, 8, 9],
[10, 11, 12]]])
In [82]: arr3d[0] = old_values
In [83]: arr3d
Out[83]:
array([[[ 1, 2, 3],
[ 4, 5, 6]],
[[ 7, 8, 9],
[10, 11, 12]]])

相似的,arr3d[1,0]也可以访问相应的数据(以⼀维数组的形式返回):

In [84]: arr3d[1, 0]
Out[84]: array([7, 8, 9])

虽然是⽤两步进⾏索引的,表达式是相同的:

In [85]: x = arr3d[1]
In [86]: x
Out[86]:
array([[ 7, 8, 9],
[10, 11, 12]])
In [87]: x[0]
Out[87]: array([7, 8, 9])

注意,在上⾯所有这些选取数组⼦集的例⼦中,返回的数组都是视图。(对视图的修改会导致原始数组的改变)

6、切片索引

ndarray的切⽚语法跟Python列表这样的⼀维对象差不多:

In [88]: arr
Out[88]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
In [89]: arr[1:6]
Out[89]: array([ 1, 2, 3, 4, 64])

对于之前的⼆维数组arr2d,其切⽚⽅式稍显不同:

In [90]: arr2d
Out[90]:
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
In [91]: arr2d[:2]
Out[91]:
array([[1, 2, 3],[4, 5, 6]])

可以看出,它是沿着第0轴(即第⼀个轴)切⽚的。也就是说,切⽚是沿着⼀个轴向选取元素的。表达式arr2d[:2]可以被认为是“选取arr2d的前两⾏”。你可以⼀次传⼊多个切⽚,就像传⼊多个索引那样:(即对前面的切片再次进行切片)

In [92]: arr2d[:2, 1:]
Out[92]:
array([[2, 3],
[5, 6]])

像这样进⾏切⽚时,只能得到相同维数的数组视图。通过将整数索引和切⽚混合,可以得到低维度的切⽚。
例如,我可以选取第⼆⾏的前两列:

In [93]: arr2d[1, :2]
Out[93]: array([4, 5])

相似的,还可以选择第三列的前两⾏:

In [94]: arr2d[:2, 2]
Out[94]: array([3, 6])

注意,“只有冒号”表示选取整个轴,因此你可以像下⾯这样只对⾼维轴进⾏切⽚:

In [95]: arr2d[:, :1]
Out[95]:
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]])

7、布尔型索引

来看这样⼀个例⼦,假设我们有⼀个⽤于存储数据的数组以及⼀个存储姓名的数组(含有重复项)。在这⾥,我将使⽤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)

这个布尔型数组可⽤于数组索引:(即只输出True的那一行数据)

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 [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 [104]: data[names == 'Bob', 2:]
Out[104]:
array([[ 0.769 , 1.2464],
[-0.5397, 0.477 ]])
In [105]: data[names == 'Bob', 3] #只选取第3列(从0开始)
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 [108]: cond = names == 'Bob'
In [109]: data[~cond]
Out[109]:
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 #将所有负值设置为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. ]])

通过⼀维布尔数组设置整⾏或列的值也很简单:

In [115]: data[names != 'Joe'] = 7
In [116]: data
Out[116]:
array([[ 7. , 7. , 7. , 7. ],
[ 1.0072, 0. , 0.275 , 0.2289],
[ 7. , 7. , 7. , 7. ],
[ 7. , 7. , 7. , 7. ],
[ 7. , 7. , 7. , 7. ],
[ 0.3026, 0.5238, 0.0009, 1.3438],
[ 0. , 0. , 0. , 0. ]])

8、花式索引

花式索引(Fancy indexing)是⼀个NumPy术语,它指的是利⽤整数数组进⾏索引。假设我们有⼀个8×4数组:

In [117]: arr = np.empty((8, 4))
In [118]: for i in range(8):
.....: arr[i] = i
In [119]: arr
Out[119]:
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 [120]: arr[[4, 3, 0, 6]]
Out[120]:
array([[ 4., 4., 4., 4.],
[ 3., 3., 3., 3.],
[ 0., 0., 0., 0.],
[ 6., 6., 6., 6.]])

这段代码确实达到我们的要求了!使⽤负数索引将会从末尾开始选取⾏:

In [121]: arr[[-3, -5, -7]]
Out[121]:
array([[ 5., 5., 5., 5.],
[ 3., 3., 3., 3.],
[ 1., 1., 1., 1.]])

⼀次传⼊多个索引数组会有⼀点特别。它返回的是⼀个⼀维数组,其中的元素对应各个索引元组:

In [122]: arr = np.arange(32).reshape((8, 4)) #reshape是指将数据重组成n行m列的数组
In [123]: arr
Out[123]:
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 [124]: arr[[1, 5, 7, 2], [0, 3, 1, 2]]
Out[124]: array([ 4, 23, 29, 10])

最终选出的是元素(1,0)、(5,3)、(7,1)和(2,2)。⽆论数组是多少维的,花式索引总是⼀维的。这个花式索引的⾏为可能会跟某些⽤户的预期不⼀样(包括我在内),选取矩阵的⾏列⼦集应该是矩形区域的形式才对。下⾯是得到该结果的⼀个办法:

In [125]: arr[[1, 5, 7, 2]][:, [0, 3, 1, 2]]
Out[125]:
array([[ 4, 7, 5, 6],
[20, 23, 21, 22],
 [28, 31, 29, 30],
[ 8, 11, 9, 10]])

记住,花式索引跟切⽚不⼀样,它总是将数据复制到新数组中。

9、数组转置和轴对称

转置是重塑的⼀种特殊形式,它返回的是源数据的视图(不会进⾏任何复制操作)。数组不仅有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)
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]]])

这⾥,第⼀个轴被换成了第⼆个,第⼆个轴被换成了第⼀个,最后⼀个轴不变。

简单的转置可以使⽤.T,它其实就是进⾏轴对换⽽已。ndarray还有⼀个swapaxes⽅法,它需要接受⼀对轴编号:

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也是返回源数据的视图(不会进⾏任何复制操作)。

二、通用函数:快速的元素级数组函数

通⽤函数(即ufunc)是⼀种对ndarray中的数据执⾏元素级运算的函数。你可以将其看做简单函数(接受⼀个或多个标量值,并产⽣⼀个或多个标量值)的⽮量化包装器。许多ufunc都是简单的元素级变体,如sqrt和exp:

In [137]: arr = np.arange(10)
In [138]: arr
Out[138]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [139]: np.sqrt(arr)
Out[139]:
array([ 0. , 1. , 1.4142, 1.7321, 2. , 2.2361, 2.4495,
2.6458, 2.8284, 3. ])
In [140]: np.exp(arr)
Out[140]:
array([ 1. , 2.7183, 7.3891, 20.0855, 54.5982,
148.4132, 403.4288, 1096.6332, 2980.958 , 8103.0839])

这些都是⼀元(unary)ufunc。另外⼀些(如add或maximum)接受2个数组(因此也叫⼆元(binary)ufunc),并返回⼀个结果数组:

In [141]: x = np.random.randn(8)
In [142]: y = np.random.randn(8)
In [143]: x
Out[143]:
array([-0.0119, 1.0048, 1.3272, -0.9193, -1.5491, 0.0222, 0.7584,
-0.6605])
In [144]: y
Out[144]:
array([ 0.8626, -0.01 , 0.05 , 0.6702, 0.853 , -0.9559, -0.0235,
-2.3042])
In [145]: np.maximum(x, y)
Out[145]:
array([ 0.8626, 1.0048, 1.3272, 0.6702, 0.853 , 0.0222, 0.7584,
-0.6605])

这⾥,numpy.maximum计算了x和y中元素级别最⼤的元素。虽然并不常⻅,但有些ufunc的确可以返回多个数组。modf就是⼀个例⼦,它是Python内置函数divmod的⽮量化版本,它会返回浮点数数组的⼩数和整数部分:

In [146]: arr = np.random.randn(7) * 5
In [147]: arr
Out[147]: array([-3.2623, -6.0915, -6.663 , 5.3731, 3.6182, 3.45 , 5.0077])
In [148]: remainder, whole_part = np.modf(arr)
In [149]: remainder
Out[149]: array([-0.2623, -0.0915, -0.663 , 0.3731,
0.6182, 0.45 , 0.0077])
In [150]: whole_part
Out[150]: array([-3., -6., -6., 5., 3., 3., 5.])

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])

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

三、利用数组进行数据处理

NumPy数组使你可以将许多种数据处理任务表述为简洁的数组表达式(否则需要编写循环)。⽤数组表达式代替循环的做法,通常被称为⽮量化。 ⼀般来说,⽮量化数组运算要⽐等价的纯Python⽅式快上⼀两个数量级(甚⾄更多),尤其是各种数值计算。 在后⾯内容中(⻅附录A)我将介绍⼴播,这是⼀种针对⽮量化计算的强⼤⼿段。

作为简单的例⼦,假设我们想要在⼀组值(⽹格型)上计算函数 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]])

作为第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函数创建的。
在这里插入图片描述

1、将条件逻辑表述为数组运算

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 [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的数组⼤⼩可以不相等,甚⾄可以是标量值。

2、数学和统计方法

可以通过数组上的⼀组数学函数对整个数组或某个轴向的数据进⾏统计计算。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)是“计算每列的和”。

axis=0代表列,axis=1代表行(第一个、二个轴)

其他如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列出了全部的基本数组统计⽅法。后续章节中有很多例⼦都会⽤到这些⽅法。
在这里插入图片描述
在这里插入图片描述

3、用于布尔型数组的方法

在上⾯这些⽅法中,布尔值会被强制转换为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和all,它们对布尔型数组⾮常有⽤。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。

4、排序

跟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

在pandas中还可以找到⼀些其他跟排序有关的数据操作(⽐如根据⼀列或多列对表格型数据进⾏排序)。

5、唯一值以及其它的集合逻辑

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 [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。
在这里插入图片描述

6、用于数组的文件输入输出

NumPy能够读写磁盘上的⽂本数据或⼆进制数据。这⼀⼩节只讨论NumPy的内置⼆进制格式,因为更多的⽤户会使⽤pandas或其它⼯具加载⽂本或表格数据(⻅第6章)。

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)

7、线性代数

线性代数(如矩阵乘法、矩阵分解、⾏列式以及其他⽅阵数学等)是任何数组库的重要组成部分。不像某些语⾔(如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中列出了⼀些最常⽤的线性代数函数。

在这里插入图片描述

8、伪随机数生成

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中的部分函数。在下⼀节中,我将给
出⼀些利⽤这些函数⼀次性⽣成⼤量样本值的范例。
在这里插入图片描述
在这里插入图片描述

9、示例:随机漫步

我们通过模拟随机漫步来说明如何运⽤数组运算。先来看⼀个简单的随机漫步的例⼦:从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])

在这里插入图片描述

图4-4 简单的随机漫步

不难看出,这其实就是随机漫步中各步的累计和,可以⽤⼀个数组运算来实现。因此,我⽤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, nstep
  • 2
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

@u@

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值