1. NumPy Basics: Arrays and Vectorized Computation

NumPy 是 Numerical Python 的简称,是高性能计算和数据分析的基础包。本书中几乎所有高级工具都是建立在它的基础之上,下面是它所能做的一些事情:

  • ndarray,快速和节省空间的多维数组,提供数组化的算术运算和高级的 广播 功能。
  • 使用标准数学函数对整个数组的数据进行快速运算,而不需要编写循环。
  • 读取/写入磁盘上的阵列数据和操作存储器映像文件的工具。
  • 线性代数,随机数生成,和傅里叶变换的能力。
  • 集成C,C++,Fortran代码的工具。

从生态系统的角度看,最后一点是最为重要的。因为NumPy 提供了易用的C API,它可以很容易的将数据传递到使用低级语言编写的外部库,也可以使外部库返回NumPy数组数据到Python。 这一特性使得Python成为包装传统的C/C++/Fortran代码库,并给它们一个动态的、易于使用的接口的首选语言。

虽然NumPy本身并没有提供非常高级的数据分析功能,但是了解NumPy的数组和面向数组的计算将会帮助你高效的使用类似于pandas这样的工具。 如果你是Python新手,并且只希望使用pandas来处理你手边的数据,随时可以略过这一章。 更多的NumPy的特性例如广播,请见 第12章

对于大多数的数据分析应用来说,我关注的主要功能是:

  • 快速的矢量化数组操作:数据切割和清除,子集和过滤,转化和任何其它类型的计算
  • 通用的数组算法,例如:sorting,unique和set操作
  • 有效的描述性统计和聚集/汇总数据
  • 数据对齐、关系数据的合并操作、异构数据的拼接操作
  • 使用数组表达式来表示条件逻辑,而不是用带有 if-elif-else 分支的循环来表示
  • 组间数据的操作(聚合,转换,功能应用)。关于这一点详见 第5章

虽然NumPy提供了这些操作的计算功能,但你或许希望使用pandas作为大多数数据分析的基础(特别是结构化或表格数据),因为它提供了一个丰富的,高级的接口使得常见的数据任务非常简洁和简单。 pandas也提供了更多的一些特定领域的功能,如时间数组操作,这是NumPy所没有的。

在本章和全书,我始终使用 import numpy as np 。当然,你喜欢在代码中使用 from numpy import * 来避免写 np ,但我要提醒你反对这种习惯。

1.1. NumPy ndarray:多维数组对象

NumPy的一个关键特性是它的N维数组对象(ndarray),它在Python中是一个大型数据集的快速的,灵活的容器。 数组使你能够在整个数据块上进行数学运算,且与对应的纯量元素间操作有相似的语法:

In [8]: data
Out[8]:
array([[ 0.9526, -0.246 , -0.8856],
       [ 0.5639, 0.2379, 0.9104]])
In [9]: data * 10                           In [10]: data + data
Out[9]:                                     Out[10]:
array([[ 9.5256, -2.4601, -8.8565],         array([[ 1.9051, -0.492 , -1.7713],
       [ 5.6385, 2.3794, 9.104 ]])                 [ 1.1277, 0.4759, 1.8208]])

ndarray是一个同种类数据的多维容器,也就是说,它的所有元素都是同类型的。每一个数组都有一个 shape (表示它每一维大小的元组)和 dtype (一个描述数组数据类型的对象):

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

本章将介绍ndarray的基础知识,并足以应对本书剩下的部分。 虽然对于许多的数据分析应用来说不必要对NumPy有深入的理解,但是精通面向数组编程和思想是成为一名科学的Python大师的关键一步。

每当你在正文中看见“array”, “NumPy array”, or “ndarray”,除了很少的列外之外,它们都指的是同一个东西:ndarray对象。

1.1.1. 创建ndarray

最简单的创建数组的方式是使用 array 函数。它接受任何数组对象(包括其它数组),产生一个包含所传递的数据的新NumPy数组。例如,列表就是一个很好的用于转换的候选:

In [13]: data1 = [6, 7.5, 8, 0, 1]
In [14]: arr1 = np.array(data1)
In [15]: arr1
Out[15]: array([ 6. , 7.5, 8. , 0. , 1. ])

嵌套序列,如等长列表的列表,将会转化为一个多维数组:

In [16]: 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]])
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('int64')

np.array 之外,还有许多函数来创建新的数组。例如, zerosones 使用给定的长度或形状分别的创建0s 和 1s数组。 empty 会创建一个没有使用特定值来初始化的数组。给这些方法传递一个元组作为形状来创建高维数组:

In [23]: np.zeros(10)
Out[23]: array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [24]: np.zeros((3, 6))
Out[24]:
array([[ 0., 0., 0., 0., 0., 0.],
       [ 0., 0., 0., 0., 0., 0.],
       [ 0., 0., 0., 0., 0., 0.]])
In [25]: np.empty((2, 3, 2))
Out[25]:
array([[[ 4.94065646e-324, 4.94065646e-324],
        [ 3.87491056e-297, 2.46845796e-130],
        [ 4.94065646e-324, 4.94065646e-324]],
       [[ 1.90723115e+083, 5.73293533e-053],
        [ -2.33568637e+124, -6.70608105e-012],
        [ 4.42786966e+160, 1.27100354e+025]]])
假定 np.array 会返回一个全零的数组是不安全的。在许多情况下,如前所示,它将返回未初始化的垃圾值。

arange 是Python内建 range 函数的数组版本:

In [26]: np.arange(15)
Out[26]: array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])

表格4-1 是一个用于构建数组的标准函数的清单。

数组构建函数
函数 描述
array 转换输入数据(列表,数组或其它序列类型)到一个ndarray,可以推断一个dtype或明确的设置一个dtype。默认拷贝输入数据。
asarray 转换输入为一个ndarray,当输入已经是一个ndarray时就不拷贝。
arange 同内建的range函数,但不返回列表而是一个ndarray
ones, ones_like 根据提供的shape和dtype产生一个全1的数组。ones_like使用另一歌数组为入参,产生一个shape和dtype都相同的数组。
zeros, zeros_like 同ones和ones_like,但是生成全0的数组
empty, enpty_like 通过分配新内存来构造新的数组,但不同与ones 和 zeros,不初始任何值。
eye, identity 生成一个NxN的单位方阵(对角线上为1,其它地方为0)

1.1.2. ndarray的数据类型

数据类型或dtype是一个特别的对象,保存了ndarray如何解释一块内存为特定类型数据的信息:

In [27]: arr1 = np.array([1, 2, 3], dtype=np.float64)
In [28]: arr2 = np.array([1, 2, 3], dtype=np.int32)
In [29]: arr1.dtype
Out[29]: dtype('float64')
In [30]: arr2.dtype
Out[30]: dtype('int32')

Dtypes是使NumPy如此强大和灵活的一部分。在大多数情况下,它们直接映射到底层的机器表示,这是的很容易地读取和写入二进制流到磁盘上,也能链接低级语言,如C 或Fortran编写的代码。数值表示的dtypes以相同的方式命名:一个类型名,如 foltint ,后面跟着一个表示数字有多少位的数字。一个标准的双精度浮点值(它是Python的 float 对象的底层表示)占据8字节或64位。因此,这一类型在NumPy中被认为是 float64 。见 表格4-2 是一个NumPy支持的全部数据类型的清单。

不要为了记忆NumPy的dtypes而烦恼,尤其你是一个新用户。通常只需要关心 你所处理数据的普通类型(浮点、复数、整形、布尔型、字符窜或一般的Python对象)。 当你需要更多的控制数据如何存储到内存和磁盘,特别是大的数据集,知道你所控制的存储类型是很好的。
NumPy数据类型
类型 类型码 描述
类型 类型码 描述
int8, uint8 i1, u1 有符号和无符号8位(1字节)整数类型
int16, uint16 i2, u2 有符号和无符号16位整数类型
int32, uint32 i4, u4 有符号和无符号32位整数类型
int64, uint64 i8, u8 有符号和无符号64位整数类型
float16 f2 半精度浮点类型
float32 f4 or f 标准精度浮点。与C的 float 兼容
float64, float128 f8 or d 标准双精度浮点。与C的 double 和Python 的 folat 对象兼容
float128 f16 or g 扩展精度浮点
complex64, complex128, complex256 c8, c16, c32 分别使用两个32,64,128位浮点表示的复数
bool ? 布尔值,存储 TrueFalse
object O Python对象类型
string_ S 定长字符窜类型(每字符一字节)。例如,为了生成长度为10的字符窜,使用 ‘S10’
unicode_ f16 or g 扩展精度浮点(字节书依赖平台)。同 string_ 有相同的语义规范(例如:U10

你可以使用ndarray的 astype 方法显示的把一个数组的dtype转换或 投射 到另外的类型:

In [31]: arr = np.array([1, 2, 3, 4, 5])
In [32]: arr.dtype
Out[32]: dtype('int64')
In [33]: float_arr = arr.astype(np.float64)
In [34]: float_arr.dtype
Out[34]: dtype('float64')

在这个例子中,整形被转换到浮点型。如果把浮点数转换到整形dtype,小数部分将会被截断:

In [35]: arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
In [36]: arr
Out[36]: array([ 3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
In [37]: arr.astype(np.int32)
Out[37]: array([ 3, -1, -2, 0, 12, 10], dtype=int32)

你可能有一个字符窜数组表示的数字,可以使用 astype 把它们转换到数字的形式:

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

如果因为某些原因(如一个字符窜不能转换到 float64 )转换失败了,将会引起一个 TypeError 。正如你所看见的,我有一些懒,使用 float 而不是 np.float64 ;NumPy会足够聪明的把Python的类型对应到等价的dtypes。

你也可以使用dtype的另一个属性:

In [40]: int_array = np.arange(10)
In [41]: calibers = np.array([.22, .270, .357, .380, .44, .50], dtype=np.float64)
In [42]: int_array.astype(calibers.dtype)
Out[42]: array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])

你也可以使用速记的类型码字符窜来指定一个dtype:

In [43]: empty_uint32 = np.empty(8, dtype='u4')
In [44]: empty_uint32
Out[44]:
array([ 0, 0, 65904672, 0, 64856792, 0,
        39438163, 0], dtype=uint32)
调用 astype 总是会创建一个新的数组(原数据的拷贝),即使是新的dtype和原来的dtype相同。
值得牢记的是浮点数,如那些是 float64float32 的数组,是唯一能够接近分数的。在复杂的计算中,可能会产生 浮点错误 ,计较时到了一定的小数位数时才有效。

1.1.3. 数组和纯量间的操作

数组非常重要,因为它们使你不使用循环就可以在数据上进行一系列操作。 这通常被叫做矢量化。相同大小的数组间的算术运算,其操作作用在对应的元素上:

In [45]: arr = np.array([[1., 2., 3.], [4., 5., 6.]])

In [46]: arr
Out[46]:
array([[ 1.,  2.,  3.],
       [ 4.,  5.,  6.]])

In [47]: arr * arr                 In [48]:arr - arr
Out[47]:                           Out[48]:
array([[  1.,   4.,   9.],        array([[ 0., 0., 0.],
       [ 16.,  25.,  36.]])              [ 0., 0., 0.]])

纯量的算术操作正如你期望的一样,把操作值作用于每一个元素:

In [49]: 1 / arr                        In [50]: arr ** 0.5
Out[49]:                                Out[50]:
array([[ 1. , 0.5 , 0.3333],            array([[ 1. , 1.4142, 1.7321],
       [ 0.25 , 0.2 , 0.1667]])                [ 2. , 2.2361, 2.4495]])

在不同大小的数组见的操作被叫做 broadcasting ,将在 第12章 详细讨论。深入的了解broadcasting在本书的多数地方是不必要的。

1.1.4. 基本的索引和切片

NumPy的索引是一个内容丰富的主题,因为有许多方法可以使你在你的数据中选取一个子集或单个元素。一维的数组很简单,表面上它们的行为类似于Python的列表:

In [51]: arr = np.arange(10)
In [52]: arr
Out[52]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [53]: arr[5]
Out[53]: 5
In [54]: arr[5:8]
Out[54]: array([5, 6, 7])
In [55]: arr[5:8] = 12
In [56]: arr
Out[56]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

如你所见,当你给一个切片赋一纯量值,如 arr[5:8] = 12 所示,该值被传送(或 传播 )到整个选择区域。与列表的第一个重要的区别是数组的切片在原来的数组上(不生成新的数组)。这意味着数据不会被拷贝,且对切片的任何修改都会影响源数组:

In [57]: arr_slice = arr[5:8]
In [58]: arr_slice[1] = 12345
In [59]: arr
Out[59]: array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8, 9])
In [60]: arr_slice[:] = 64
In [61]: arr
Out[61]: array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])

如果你是使用NumPy的新手,这一点回事你感到惊讶,尤其当你使用过其它数组编程语言,它们非常热衷于拷贝数据。请记住,NumPy是设计用来处理大数据的情况,你可以想象如果NumPy坚持使用拷贝数据将会出现的性能和内存问题。

如果你想有数组切片的一个拷贝,你需要明显的拷贝数组;例如 arr[5:8].copy()

对于高维数组,你会有更多选项。在两维的数组,每一个索引的元素将不再是一个纯量,而是一个一维数组:

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

因此,单个元素可以递归的访问,但是这会做多一点的工作。不过,你可以使用一个逗号分隔的索引列表来选择单个元素。因此,下面的操作是等价的:

In [64]: arr2d[0][2]
Out[64]: 3
In [65]: arr2d[0, 2]
Out[65]: 3

NumPy数组的索引,是在二维数组上的索引图例。

_images/index-2d.png

NumPy数组的索引

在多维数组中,如果你省略了后面的索引,返回的对象将会是一个较低维的ndarray,它包括较高维度的所有数据。因此,在 2*2*3 的数组 arr3d

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

arr3d[0] 是一个 2*3 的数组:

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

纯量值和数组都可以给 arr3d[0] 赋值:

In [69]: old_values = arr3d[0].copy()
In [70]: arr3d[0] = 42

In [71]: arr3d
Out[71]:
array([[[42, 42, 42],
        [42, 42, 42]],
       [[ 7, 8, 9],
        [10, 11, 12]]])

In [72]: arr3d[0] = old_values

In [73]: arr3d
Out[73]:
array([[[ 1, 2, 3],
        [ 4, 5, 6]],
       [[ 7, 8, 9],
        [10, 11, 12]]])

类似的, arr3d[1, 0] 给你那些索引以 (1, 0) 开始的值,形成了一个1维数组:

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

请注意,在所有的情况下,被选中的子节返回的数组总是数组视窗。

1.1.4.1. 带切片的索引

如同一维对象,例如Python的列表,ndarrys可以使用熟悉的语法来切片:

In [75]: arr[1:6]
Out[75]: array([ 1, 2, 3, 4, 64])

较高维的对象给你更多的选择,你可以切割一个或多个坐标坐标轴,并且可以混合整数。对上面的2维数组, arr2d ,对它的切片有些不同:

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

正如你所见,它沿着 0 坐标坐标轴(第一个坐标坐标轴)切片。因此,一个切片沿着一个坐标坐标轴向选择一个范围的元素。你可以传递多个切片,就像你传递多个索引一样:

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

象这样切片时,你得到的总是相同维数的数组视窗。通过混合整形索引和切片,你可以得到较低维的切片:

In [79]: arr2d[1, :2]             In [80]: arr2d[2, :1]
Out[79]: array([4, 5])            Out[80]: array([7])

两维数组切片 图解。注意,一个单一的冒号意味着取整个坐标/坐标轴,因此,你可以只切割更高维的坐标轴,做法如下:

In [81]: arr2d[:, :1]
Out[81]:
array([[1],
       [4],
       [7]])

当然,给一个切片表达式赋值会对整个选择赋值:

In [82]: arr2d[:2, 1:] = 0
_images/slicing-2d.png

两维数组切片

1.1.5. 布尔索引

让我们来考虑一个例子,我们有一些数据在一个数组中和一个有重复名字的数组。我将会在这使用 numpy.random 中的 randn 函数来产生一些随机的正态分布的数据:

In [83]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
In [84]: data = randn(7, 4)
In [85]: names
Out[85]:
array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'],
      dtype='|S4')
In [86]: data
Out[86]:
array([[-0.048 , 0.5433, -0.2349, 1.2792],
       [-0.268 , 0.5465, 0.0939, -2.0445],
       [-0.047 , -2.026 , 0.7719, 0.3103],
       [ 2.1452, 0.8799, -0.0523, 0.0672],
       [-1.0023, -0.1698, 1.1503, 1.7289],
       [ 0.1913, 0.4544, 0.4519, 0.5535],
       [ 0.5994, 0.8174, -0.9297, -1.2564]])

假设每一个名字都和 data 数组中的一行对应。如果我们想要选择与 ‘Bob’ 名字对应的所有行。象算术运算一样,数组的比较操作(例如 == )也可以矢量化。因此, namesBob 字符窜的比较会产生一个布尔数组:

In [87]: names == 'Bob'
Out[87]: array([ True, False, False, True, False, False, False], dtype=bool)

当索引数组时可以传递这一布尔数组:

In [88]: data[names == 'Bob']
Out[88]:
array([[-0.048 , 0.5433, -0.2349, 1.2792],
       [ 2.1452, 0.8799, -0.0523, 0.0672]])

布尔数组必须和它索引的坐标轴的长度相同。你甚至可以把布尔数组和切片或整数(或者整数序列,关于这一点后面会更多介绍)混合和匹配起来:

In [89]: data[names == 'Bob', 2:]
Out[89]:
array([[-0.2349, 1.2792],
       [-0.0523, 0.0672]])
In [90]: data[names == 'Bob', 3]
Out[90]: array([ 1.2792, 0.0672])

为了选择除了 ‘Bob’ 之外的所有东西,你可以使用 != 或用 - 对条件表达式取反:

In [91]: names != 'Bob'
Out[91]: array([False, True, True, False, True, True, True], dtype=bool)

In [92]: data[-(names == 'Bob')]
Out[92]:
array([[-0.268 , 0.5465, 0.0939, -2.0445],
       [-0.047 , -2.026 , 0.7719, 0.3103],
       [-1.0023, -0.1698, 1.1503, 1.7289],
       [ 0.1913, 0.4544, 0.4519, 0.5535],
       [ 0.5994, 0.8174, -0.9297, -1.2564]])

使用布尔算术操作符如 & (and) 和 | (or)来结合多个布尔条件,下面是从三个名字中选取两个的操作:

In [93]: mask = (names == 'Bob') | (names == 'Will')
In [94]: mask
Out[94]: array([True, False, True, True, True, False, False], dtype=bool)
In [95]: data[mask]
Out[95]:
array([[-0.048 , 0.5433, -0.2349, 1.2792],
       [-0.047 , -2.026 , 0.7719, 0.3103],
       [ 2.1452, 0.8799, -0.0523, 0.0672],
       [-1.0023, -0.1698, 1.1503, 1.7289]])

通过布尔索引从一个数组中选取数据 总是 会创建数据的一份拷贝,即使是返回的数组没有改变。

Python的 andor 关键字不能与布尔数组一起工作。

通过布尔数组设置值工作于一种种常识性的方式。为了设置 data 中所有的负值为0,我们只需要:

In [96]: data[data < 0] = 0
In [97]: data
Out[97]:
array([[ 0. , 0.5433, 0. , 1.2792],
       [ 0. , 0.5465, 0.0939, 0. ],
       [ 0. , 0. , 0.7719, 0.3103],
       [ 2.1452, 0.8799, 0. , 0.0672],
       [ 0. , 0. , 1.1503, 1.7289],
       [ 0.1913, 0.4544, 0.4519, 0.5535],
       [ 0.5994, 0.8174, 0. , 0. ]])

使用一维布尔数组设置整行或列也非常简单:

In [98]: data[names != 'Joe'] = 7
In [99]: data
Out[99]:
array([[ 7. , 7. , 7. , 7. ],
       [ 0. , 0.5465, 0.0939, 0. ],
       [ 7. , 7. , 7. , 7. ],
       [ 7. , 7. , 7. , 7. ],
       [ 7. , 7. , 7. , 7. ],
       [ 0.1913, 0.4544, 0.4519, 0.5535],
       [ 0.5994, 0.8174, 0. , 0. ]])

1.1.6. Fancy索引

Fancy 索引 是一个术语,被NumPy用来描述使用整形数组索引。假如我们有一个 8*4 的数组:

In [100]: arr = np.empty((8, 4))
In [101]: for i in range(8):
   .....:     arr[i] = i
In [102]: arr
Out[102]:
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 [103]: arr[[4, 3, 0, 6]]
Out[103]:
array([[ 4., 4., 4., 4.],
       [ 3., 3., 3., 3.],
       [ 0., 0., 0., 0.],
       [ 6., 6., 6., 6.]])

很庆幸这个代码做了你所期望的!使用负的索引从结尾选择行:

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

传递多个索引数组有些微的不同;它选取一个一维数组,元素对应与索引的每一个元组:

# 关于reshape在第12章会跟多介绍
In [105]: arr = np.arange(32).reshape((8, 4))
In [106]: arr
Out[106]:
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 [107]: arr[[1, 5, 7, 2], [0, 3, 1, 2]]
Out[107]: array([ 4, 23, 29, 10])

花一点儿时间来看看刚刚发生了什么:元素 (1, 0), (5, 3), (7, 1), 和(2, 2)被选择了。 fancy索引的行为与一些用户(也包括我自己)可能期望的有所不同, 它因该是一个矩形区域,由选取的矩形的行和列组成。这里有一个方法来得到它:

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

另一种方法是使用 np.ix_ 函数,将两个以为整形数组转换为位标,来选取一个正方形区域:

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

注意,fancy索引,不像切片,它总是拷贝数据到一个新的数组。

1.1.7. 转置数组和交换坐标轴

转置是一种特殊形式的变形,类似的它会返回基础数据的一个视窗,而不会拷贝任何东西。数组有 transpose 方法和专门的 T 属性:

In [110]: arr = np.arange(15).reshape((3, 5))
In [111]: arr                          In [112]: arr.T
Out[111]:                              Out[112]:
array([[ 0, 1, 2, 3, 4],               array([[ 0, 5, 10],
       [ 5, 6, 7, 8, 9],                      [ 1, 6, 11],
       [10, 11, 12, 13, 14]])                 [ 2, 7, 12],
                                              [ 3, 8, 13],
                                              [ 4, 9, 14]])

当进行矩阵运算时,你常常会这样做,像下面的例子一样,使用 np.dot 计算内部矩阵来产生 XTX` :

In [113]: arr = np.random.randn(6, 3)
In [114]: np.dot(arr.T, arr)
Out[114]:
array([[ 2.584 , 1.8753, 0.8888],
       [ 1.8753, 6.6636, 0.3884],
       [ 0.8888, 0.3884, 3.9781]])

对于更高维的数组, transpose 接受用于转置的坐标轴的号码的一个元组(for extra mind bending):

In [115]: arr = np.arange(16).reshape((2, 2, 4))
In [116]: arr
Out[116]:
array([[[ 0, 1, 2, 3],
        [ 4, 5, 6, 7]],
       [[ 8, 9, 10, 11],
        [12, 13, 14, 15]]])

In [117]: arr.transpose((1, 0, 2))
Out[117]:
array([[[ 0, 1, 2, 3],
        [ 8, 9, 10, 11]],
       [[ 4, 5, 6, 7],
        [12, 13, 14, 15]]])

使用 .T 的转置,仅仅是交换坐标轴的一个特殊的情况:

In [118]: arr                             In [119]: arr.swapaxes(1, 2)
Out[118]:                                 Out[119]:
array([[[ 0, 1, 2, 3],                    array([[[ 0, 4],
        [ 4, 5, 6, 7]],                           [ 1, 5],
                                                  [ 2, 6],
       [[ 8, 9, 10, 11],                          [ 3, 7]],
        [12, 13, 14, 15]]])
                                                 [[ 8, 12],
                                                 [ 9, 13],
                                                 [10, 14],
                                                 [11, 15]]])

类似的 swapaxes 返回在数据上的一个视窗,而不进行拷贝。

1.2. 通用函数:快速的基于元素的数组函数

一个通用的函数,或者 ufunc ,是一个在ndarrays的数据上进行基于元素的操作的函数。你可以认为它们是对简单函数的一个快速矢量化封装,它们接受一个或多个标量值并产生一个或多个标量值。

许多 ufuncs 都是基于元素的简单变换,像 sqrtexp

In [120]: arr = np.arange(10)
In [121]: np.sqrt(arr)
Out[121]:
array([ 0. , 1. , 1.4142, 1.7321, 2. , 2.2361, 2.4495,
        2.6458, 2.8284, 3. ])
In [122]: np.exp(arr)
Out[122]:
array([ 1. , 2.7183, 7.3891, 20.0855, 54.5982,
       148.4132, 403.4288, 1096.6332, 2980.958 , 8103.0839])

这些归诸于 unary ufuncs。其它的,例如 addmaximum ,接受两个数组(因此,叫做 binary ufuncs)且返回一个数组:

In [123]: x = randn(8)
In [124]: y = randn(8)
In [125]: x
Out[125]:
array([ 0.0749, 0.0974, 0.2002, -0.2551, 0.4655, 0.9222, 0.446 ,
       -0.9337])
In [126]: y
Out[126]:
array([ 0.267 , -1.1131, -0.3361, 0.6117, -1.2323, 0.4788, 0.4315,
       -0.7147])
In [127]: np.maximum(x, y) # element-wise maximum
Out[127]:
array([ 0.267 , 0.0974, 0.2002, 0.6117, 0.4655, 0.9222, 0.446 ,
       -0.7147])

虽然不常见,一个ufunc可以返回多个数组。 nodf 就是一个例子,它是Python内建 divmod 的矢量化的版本:它返回一个副点数数组的分数和整数部分:

In [128]: arr = randn(7) * 5
In [129]: np.modf(arr)
Out[129]:
(array([-0.6808, 0.0636, -0.386 , 0.1393, -0.8806, 0.9363, -0.883 ]),
 array([-2., 4., -3., 5., -3., 3., -6.]))

表格4-3表格4-4 是可用的ufuncs的清单。

Unary ufuncs
函数 描述
abs, fabs 计算基于元素的整形,浮点或复数的绝对值。fabs对于没有复数数据的快速版本
sqrt 计算每个元素的平方根。等价于 arr ** 0.5
square 计算每个元素的平方。等价于 arr ** 2
exp 计算每个元素的指数。
log, log10, log2, log1p 自然对数(基于e),基于10的对数,基于2的对数和 log(1 + x)
sign 计算每个元素的符号:1(positive),0(zero), -1(negative)
ceil 计算每个元素的天花板,即大于或等于每个元素的最小值
floor 计算每个元素的地板,即小于或等于每个元素的最大值
rint 圆整每个元素到最近的整数,保留dtype
modf 分别返回分数和整数部分的数组
isnan 返回布尔数组标识哪些元素是 NaN (不是一个数)
isfinite, isinf 分别返回布尔数组标识哪些元素是有限的(non-inf, non-NaN)或无限的
cos, cosh, sin sinh, tan, tanh regular 和 hyperbolic 三角函数
arccos, arccosh, arcsin, arcsinh, arctan, arctanh 反三角函数
logical_not 计算基于元素的非x的真值。等价于 -arr
Binary universal funcitons
函数 描述
add 在数组中添加相应的元素
substract 在第一个数组中减去第二个数组
multiply 对数组元素相乘
divide, floor_divide 除和地板除(去掉余数)
power 使用第二个数组作为指数提升第一个数组中的元素
maximum, fmax 基于元素的最大值。 fmax 忽略 NaN
minimum, fmin 基于元素的最小值。 fmin 忽略 NaN
mod 基于元素的模(取余)
copysign 拷贝第二个参数的符号到第一个参数
greater, greater_equal, less, less_equal, not_equal 基于元素的比较,产生布尔数组。等价于中缀操作符 >, >=, <, <=, ==, !=
logical_and, logical_or, logical_xor 计算各个元素逻辑操作的真值。等价于中缀操作符 &, |, ^

1.3. 使用数组进行数据处理

使用NumPy可以是你能够使用简明的数组表达式而不是编写循环表达许多种类的数据处理任务。这种使用数组表达式代替显示循环通常被成为“矢量化”。在一般情况下,矢量化数组操作比与之等价的纯Python操作数度快一到两(或更多)个等级,这对任何种类的数值计算有最大的影响。稍后,在chp12index中,我会讲解 broadcasting ,一个矢量化计算的强大方法。

作为一个简单示例,假如我们希望研究函数 sqrt(x\ :sup:`^`\ 2 + \ :sup:`^`\ 2) 穿过一个网格数据。 np.meshgrid 函数接受两个一维数组并产生两个二维矩阵,其值对于两个数组的所有 (x, y) 对:

In [130]: points = np.arange(-5, 5, 0.01) # 1000个等间隔点
In [131]: xs, ys = np.meshgrid(points, points)
In [132]: ys
Out[132]:
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 [134]: import matplotlib.pyplot as plt
In [135]: z = np.sqrt(xs ** 2 + ys ** 2)
In [136]: z
Out[136]:
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]])
In [137]: plt.imshow(z, cmap=plt.cm.gray); plt.colorbar()
Out[137]: <matplotlib.colorbar.Colorbar instance at 0x4e46d40>
In [138]: plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")
Out[138]: <matplotlib.text.Text at 0x4565790>

绘制在网格上的函数,我使用 matplotlib 函数 imshow 创建一个了一个图像,数据来源于上面的函数生成的二维数组。

_images/plot-grid.png

绘制在网格上的函数

1.3.1. 用数组操作来表达条件逻辑

函数 numpy.where 是三元表达式 x if condition else y 的矢量化版本。假如我们有一个布尔数组和两个值数组:

In [140]: xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
In [141]: yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
In [142]: cond = np.array([True, False, True, True, False])

假如我们想要当对应的 cond 值为 True 时,从 xarr 中获取一个值,否则从 yarr 中获取值。使用列表推到来做这件事,可能会像这样:

In [143]: result = [(x if c else y)
   .....: for x, y, c in zip(xarr, yarr, cond)]
In [144]: result
Out[144]: [1.1000000000000001, 2.2000000000000002, 1.3, 1.3999999999999999, 2.5]

这样做会有许多问题。首先,对于大的数组,它不会很快(因为所有的工作都是有纯Python来做的)。其次,对于多维数组,它不能工作。使用 np.where 你可以像这样非常简洁的编写:

In [145]: result = np.where(cond, xarr, yarr)
In [146]: result
Out[146]: array([ 1.1, 2.2, 1.3, 1.4, 2.5])

np.where 的第一个和第二个参数不需要是数组;它们中的一个或两个可以是纯量。 在数据分析中 where 的典型使用是生成一个新的数组,其值基于另一个数组。假如你有一个矩阵,其数据是随机生成的,你想要把其中的正值替换为2,负值替换为-2,使用 np.where 非常容易:

In [147]: arr = randn(4, 4)
In [148]: arr
Out[148]:
array([[ 0.6372, 2.2043, 1.7904, 0.0752],
       [-1.5926, -1.1536, 0.4413, 0.3483],
       [-0.1798, 0.3299, 0.7827, -0.7585],
       [ 0.5857, 0.1619, 1.3583, -1.3865]])
In [149]: np.where(arr > 0, 2, -2)
Out[149]:
array([[ 2, 2, 2, 2],
       [-2, -2, 2, 2],
       [-2, 2, 2, -2],
       [ 2, 2, 2, -2]])

In [150]: np.where(arr > 0, 2, arr) # 仅设置正值为 2
Out[150]:
array([[ 2. , 2. , 2. , 2. ],
       [-1.5926, -1.1536, 2. , 2. ],
       [-0.1798, 2. , 2. , -0.7585],
       [ 2. , 2. , 2. , -1.3865]])

传递到 where 的数组不仅仅只是大小相等的数组或纯量。

使用一些小聪明,你可以使用 where 来表达更复杂的逻辑;考虑这个例子,我有两个布尔数组, cond1cond2 ,并想根据4种布尔值来赋值:

result = []
for i in range(n):
    if cond1[i] and cond2[i]:
        result.append(0)
    elif cond1[i]:
        result.append(1)
    elif cond2[i]:
        result.append(2)
    else:
        result.append(3)

也许可能不会很明显,这个 for 循环可以转换成一个嵌套的 where 表达式:

np.where(cond1 & cond2, 0,
         np.where(cond1, 1,
                  np.where(cond2, 2, 3)))

在这个特殊的例子中,我们还可以利用布尔表达式在计算中被当作0或1这一事实,因此可以使用算术运算来表达:

result = 1 * cond1 + 2 * cond2 + 3 * -(cond1 | cond2)

1.3.2. 数学和统计方法

一组数学函数,计算整个数组或一个轴向上数据的统计,和数组函数一样是容易访问的。聚合(通常被称为 reductions ),如 sunmean ,标准偏差 std 可以使用数组实例的方法,也可以使用顶层NumPy的函数:

In [151]: arr = np.random.randn(5, 4) # 正态分布数据
In [152]: arr.mean()
Out[152]: 0.062814911084854597
In [153]: np.mean(arr)
Out[153]: 0.062814911084854597
In [154]: arr.sum()
Out[154]: 1.2562982216970919

meansun 函数可以有一个可选的 axis 参数,它对给定坐标轴进行统计,结果数组将会减少一个维度:

In [155]: arr.mean(axis=1)
Out[155]: array([-1.2833, 0.2844, 0.6574, 0.6743, -0.0187])
In [156]: arr.sum(0)
Out[156]: array([-3.1003, -1.6189, 1.4044, 4.5712])

cumsumcumprod 这些函数并不聚集,而是产生一个 intermediate results 的数组:

In [157]: arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [158]: arr.cumsum(0)         In [159]: arr.cumprod(1)
Out[158]:                       Out[159]:
array([[ 0, 1, 2],              array([[ 0, 0, 0],
[ 3, 5, 7],                            [ 3, 12, 60],
[ 9, 12, 15]])                         [ 6, 42, 336]])

表格4-5 是一个完整的清单。我们将在稍后的章节中看见关于这些函数的大量例子。

数组构建函数
方法 描述
sum 对数组的所有或一个轴向上的元素求和。零长度的数组的和为灵。
mean 算术平均值。灵长度的数组的均值为NaN。
std, var 标准差和方差,有可选的调整自由度(默认值为n)。
min, max 最大值和最小值
argmin, argmax 索引最小和最大元素。
cumsum 从0元素开始的累计和。
cumprod 从1元素开始的累计乘。

1.3.3. 布尔数组的方法

在上面的方法中布尔值被强制为1( True )和0a( False )。因此, sum 经常被用来作为对一个布尔数组中的 True 计数的手段:

In [160]: arr = randn(100)
In [161]: (arr > 0).sum() # 正值的个数
Out[161]: 44

有两个额外的方法, anyall ,对布尔数组尤其有用。 any 用来测试一个数组中是否有一个或更多的 True ,而 all 用来测试所有的值是否为 True

In [162]: bools = np.array([False, False, True, False])
In [163]: bools.any()
Out[163]: True
In [164]: bools.all()
Out[164]: False

这些方法这些方法也可以工作在非不而数组上,非零元素作为 True

1.3.4. 排序

像Python的内建列表一样,NumPy数组也可以使用 sort 方法就地排序:

In [165]: arr = randn(8)
In [166]: arr
Out[166]:
array([ 0.6903, 0.4678, 0.0968, -0.1349, 0.9879, 0.0185, -1.3147,
       -0.5425])
In [167]: arr.sort()
In [168]: arr
Out[168]:
array([-1.3147, -0.5425, -0.1349, 0.0185, 0.0968, 0.4678, 0.6903,
        0.9879])

多维数组可以通过传递一个坐标轴数到 sort ,对一维截面上的数据进行就地排序:

In [169]: arr = randn(5, 3)
In [170]: arr
Out[170]:
array([[-0.7139, -1.6331, -0.4959],
       [ 0.8236, -1.3132, -0.1935],
       [-1.6748, 3.0336, -0.863 ],
       [-0.3161, 0.5362, -2.468 ],
       [ 0.9058, 1.1184, -1.0516]])
In [171]: arr.sort(1)
In [172]: arr
Out[172]:
array([[-1.6331, -0.7139, -0.4959],
       [-1.3132, -0.1935, 0.8236],
       [-1.6748, -0.863 , 3.0336],
       [-2.468 , -0.3161, 0.5362],
       [-1.0516, 0.9058, 1.1184]])

顶层的 np.sort 函数返回一个经过排序后的数组拷贝,而不是就地修改。一个快速和肮脏的计算一个数组的位数是对它排序并选择一个特定阶层值:

In [173]: large_arr = randn(1000)
In [174]: large_arr.sort()
In [175]: large_arr[int(0.05 * len(large_arr))] # 5% quantile
Out[175]: -1.5791023260896004

关于使用NumPy的排序方法和更高级的技术,如间接排序,请见第12章。其它几种有关排序的数据操作(例如,通过一列或多列对数据表排序)也会在 pandas 中找到。

1.3.5. Unique 和其它集合逻辑

Numpy有一些基本的针对一维ndarrays的集合操作。最常使用的一个可能是 np.unique ,它返回一个数组的经过排序的 unique 值:

In [176]: names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
In [177]: np.unique(names)
Out[177]:
array(['Bob', 'Joe', 'Will'],
      dtype='|S4')
In [178]: ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4])
In [179]: np.unique(ints)
Out[179]: array([1, 2, 3, 4])

np.unique 与纯Python版本比较:

In [180]: sorted(set(names))
Out[180]: ['Bob', 'Joe', 'Will']

另一个函数 np.in1d ,测试一个数组的值和另一个的关系,返回一个布尔数组:

In [181]: values = np.array([6, 0, 0, 3, 2, 5, 6])
In [182]: np.in1d(values, [2, 3, 6])
Out[182]: array([ True, False, False, True, True, False, True], dtype=bool)

表格4-6 是关于集合函数的清单。

数组集合操作
unique(x) 计算x单一的元素,并对结果排序
intersect1d(x, y) 计算x和y相同的元素,并对结果排序
union1d 结合x和y的元素,并对结果排序
in1d(x, y) 得到一个布尔数组指示x中的每个元素是否在y中
setdiff1d(x, y) 差集,在x中但不再y中的集合
setxor1d(x, y) 对称差集,不同时在两个数组中的元素

1.4. 关于数组的文件输入和输出

NumPy能够保存数据到磁盘和从磁盘加载数据,不论数据是文本或二进制的。在后面的章节你可以学到使用pandas提供的工具来加载表格化的数据到内存。

1.4.1. 对磁盘上的二进制格式数组排序

np.savenp.load 是两个主力功能,有效的保存和加载磁盘数据。数组默认保存为未经过压缩的原始二进制数据,文件扩展名为 .npy

In [183]: arr = np.arange(10)
In [184]: np.save('some_array', arr)

如果文件路进并不是以 .npy 结尾,扩展名将会被自动加上。在磁盘上的数组可以使用 np.load 加载:

In [185]: np.load('some_array.npy')
Out[185]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

你可以使用 np.savez 并以关键字参数传递数组来保存多个数组到一个zip的归档文件中:

In [186]: np.savez('array_archive.npz', a=arr, b=arr)

当你加载一个 .npz 文件时,会得到一个字典对象,它懒洋洋的加载单个数组:

In [187]: arch = np.load('array_archive.npz')
In [188]: arch['b']
Out[188]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

1.4.2. 保存和加载文本文件

从文件加载文本是一个相当标准的任务。对一个新人来说,Python的文件加读取和写入函数的景象可能有一点儿混乱,因此我将主要集中在pandas的 read_csvread_table 函数上。有时使用 np.loadtxt 或更专门的 np.genfromtxt 对于加载数据到 vanilla NumPy 数组是很有用的。

这些函数有许多选项,允许你指定不同的分割副,特定列的转换函数,跳过某些行,和其它的事情。以这样一个逗号分割文件(CSV)作为一个简单的例子:

In [191]: !cat array_ex.txt
0.580052,0.186730,1.040717,1.134411
0.194163,-0.636917,-0.938659,0.124094
-0.126410,0.268607,-0.695724,0.047428
-1.484413,0.004176,-0.744203,0.005487
2.302869,0.200131,1.670238,-1.881090
-0.193230,1.047233,0.482803,0.960334

它可以像这样被加载到一个二维数组:

In [192]: arr = np.loadtxt('array_ex.txt', delimiter=',')
In [193]: arr
Out[193]:
array([[ 0.5801, 0.1867, 1.0407, 1.1344],
[ 0.1942, -0.6369, -0.9387, 0.1241],
[-0.1264, 0.2686, -0.6957, 0.0474],
[-1.4844, 0.0042, -0.7442, 0.0055],
[ 2.3029, 0.2001, 1.6702, -1.8811],
[-0.1932, 1.0472, 0.4828, 0.9603]])

np.savatxt 执行相反的操作:写入数组到一个界定文本文件中。 genfromtxtloadtxt 相似,但是她是面向结构数组和缺失数据处理的;更多关于结构数组请见第12章

更多有关读取和写入,特别是表格式的或类电子表格的数据,见后面涉及到pandas和DataFrame对象的章节。

1.5. 线性代数

线性代数,如矩阵乘法,分解,行列式和其它的方阵数学,对任何一个数组库来说都是重要的部分。不像一些语言,如 MATLAB ,使用 * 来乘两个二维数组是基于元素的乘法,而不是矩阵点积。因此,有一个 dot 函数,是数组的一个方法和 numpy 命名空间中的一个函数,用来进行矩阵乘法运算:

In [194]: x = np.array([[1., 2., 3.], [4., 5., 6.]])
In [195]: y = np.array([[6., 23.], [-1, 7], [8, 9]])
In [196]: x                In [197]: y
Out[196]:                  Out[197]:
array([[ 1., 2., 3.],      array([[ 6., 23.],
       [ 4., 5., 6.]])            [ -1., 7.],
                                  [ 8., 9.]])
In [198]: x.dot(y) # equivalently np.dot(x, y)
Out[198]:
array([[ 28., 64.],
       [ 67., 181.]])

在一个二维数组和合适大小的一维数组间的矩阵乘积的结果是一个一维数组:

In [199]: np.dot(x, np.ones(3))
Out[199]: array([ 6., 15.])

numpy.linalg 有一个关于矩阵分解和像转置和行列式等的一个标准集合。它们和其它语言(如: MATLABR )一样都是基于行业标准的 Fortran 库,如 BLSALAPACK ,或可能的 Intel MKL (依赖于你的NumPy的编译)实现的:

In [201]: from numpy.linalg import inv, qr
In [202]: X = randn(5, 5)
In [203]: mat = X.T.dot(X)
In [204]: inv(mat)
Out[204]:
array([[ 3.0361, -0.1808, -0.6878, -2.8285, -1.1911],
       [-0.1808, 0.5035, 0.1215, 0.6702, 0.0956],
       [-0.6878, 0.1215, 0.2904, 0.8081, 0.3049],
       [-2.8285, 0.6702, 0.8081, 3.4152, 1.1557],
       [-1.1911, 0.0956, 0.3049, 1.1557, 0.6051]])
In [205]: mat.dot(inv(mat))
Out[205]:
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 [206]: q, r = qr(mat)
In [207]: r
Out[207]:
array([[ -6.9271, 7.389 , 6.1227, -7.1163, -4.9215],
       [ 0. , -3.9735, -0.8671, 2.9747, -5.7402],
       [ 0. , 0. , -10.2681, 1.8909, 1.6079],
       [ 0. , 0. , 0. , -1.2996, 3.3577],
       [ 0. , 0. , 0. , 0. , 0.5571]])

表格4-7 是一些常用的线性代数常用的函数清单。

科学Python社区希望有一天可以实现矩阵乘法的中缀操作符,提供一个语法上更好的使用 np.dot 的替代。但是现在只能这样做。
常用 numpy.linglg 函数
函数 描述
diag 返回一个方阵的对角线(或非对角线)元素为一个一维数组,或者转换一个一维数组到一个方阵(非对角线元素为零)
dot 矩阵乘积
trace 计算对角线上元素的和
det 计算矩阵行列式
eig 计算方阵的特征值和特征向量
inv 计算方阵转置
pinv 计算方阵 Moore-Penrose pseudo-inverse 的转置
qr 计算 QR 分解
svd 计算奇异值分解(SVD)
solve 求解线性系统方程 Ax = b 的x,其中A是一个方阵
lstsq 计算 y = Xb 的最小二乘解

1.6. 示例:随机游走

这是一个利用数组操作来模拟随机游走的示例程序。让我们先来看一个简单的随机游走的例子,从0开始,步长为1和-1,且以相等的概率出现。一个纯Python方式来实现一个单一的有1000步的随机游走的方式是使用内建的 random 模块:

import random
position = 0
walk = [position]
steps = 1000
for i in xrange(steps):
    step = 1 if random.randint(0, 1) else -1
    position += step
    walk.append(position)

一个简单的随机游走是使用这些随机游走的前100个值的例图。

_images/random-walk.png

一个简单的随机游走

你可能会发现 walk 简单的把随机步长累积起来并且可以可以使用一个数组表达式来计算。因此,我用 np.random 模块去1000次硬币翻转,设置它们为1和-1,并计算累计和:

In [215]: nsteps = 1000
In [216]: draws = np.random.randint(0, 2, size=nsteps)
In [217]: steps = np.where(draws > 0, 1, -1)
In [218]: walk = steps.cumsum()

从这,我们可以开始沿着游走轨迹来提取如最小或做大值的统计信息:

In [219]: walk.min() In [220]: walk.max()
Out[219]: -3 Out[220]: 31

一个更复杂的统计数据是第一交叉时间,随机游走达到一个特定值的步值。这里,我们可能想要知道过了多长时间的随机游走,从任一个方向到达距离原点0至少10步之遥。 ** np.ads(walk) >= 10 ** 会给我们一个布尔数组指示在哪儿游走到达了或超过了10,但是我需要的是第一个10或-10的索引。可以使用 argmax 来计算,它返回布尔数组(最大值为 True)中第一个最大值的索引:

In [221]: (np.abs(walk) >= 10).argmax()
Out[221]: 37

注意在这使用 ragmax 并不是总是高效的,因为它总是对数组做全扫描。在这一特殊情况下,一旦一个 True 出现了,我们就知道它是一个最大值。

1.6.1. 一次模拟许多随机游走

如果你的目标是模拟许多随机游走,如5000个,你可以对上面的代码稍作修改来生成所有的随机游动。 numpy.random 函数,如果通过一个2元组,将产生一个二维数组绘制,我们可以跨越行一次计算5000个随机游动的累计和:

In [222]: nwalks = 5000
In [223]: nsteps = 1000
In [224]: draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1
In [225]: steps = np.where(draws > 0, 1, -1)
In [226]: walks = steps.cumsum(1)
In [227]: walks
Out[227]:
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 [228]: walks.max() In [229]: walks.min()
Out[228]: 138 Out[229]: -133

在这些游走中,让我们来计算到达30或-30的最短时间。这有一点儿狡猾,因为不是所有的5000个游走都能到达30。我们可以使用 any 方法来检测:

In [230]: hits30 = (np.abs(walks) >= 30).any(1)
In [231]: hits30
Out[231]: array([False, True, False, ..., False, True, False], dtype=bool)
In [232]: hits30.sum() # 30或-30的个数
Out[232]: 3410

我们可以使用这个布尔数组来选择这些游走中跨过绝对30的行,并调用 argmax 来取得坐标轴1的交叉时间:

In [233]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
In [234]: crossing_times.mean()
Out[234]: 498.88973607038122

可以大胆的试验其它的分布的步长,而不是相等大小的硬币翻转。你只需要使用一个不同的随机数生成函数,如 normal 来产生相同均值和标准偏差的正态分布:

In [235]: steps = np.random.normal(loc=0, scale=0.25,
   .....: size=(nwalks, nsteps))