第四篇 NumPy基础:数组和⽮量计算

NumPy(Numerical Python的简称)是Python数值计算最重要的基础包。⼤多数提供科学计算的包都是⽤NumPy的数组作为构建基础。
NumPy的部分功能如下:
  ndarray,⼀个具有⽮量算术运算和复杂⼴播能⼒的快速且节省空间的多维数组。
  ⽤于对整组数据进⾏快速运算的标准数学函数(⽆需编写循环)。
  ⽤于读写磁盘数据的⼯具以及⽤于操作内存映射⽂件的⼯具。
  线性代数、随机数⽣成以及傅⾥叶变换功能。
  ⽤于集成由C、C++、Fortran等语⾔编写的代码的A C API。
由于NumPy提供了⼀个简单易⽤的C API,因此很容易将数据传递给由低级语⾔编写的外部库,外部库也能以NumPy数组的形式将数据返回给Python。这个功能使Python成为⼀种包装C/C++/Fortran历史代码库的选择,并使被包装库拥有⼀个动态的、易⽤的接⼝。

⼤部分数据分析应⽤,最关注的功能主要集中在:
  ⽤于数据整理和清理、⼦集构造和过滤、转换等快速的⽮量化数组运算。
  常⽤的数组算法,如排序、唯⼀化、集合运算等。
  ⾼效的描述统计和数据聚合/摘要运算。
  ⽤于异构数据集的合并/连接运算的数据对⻬和关系型数据运算。
  将条件逻辑表述为数组表达式(⽽不是带有if-elif-else分⽀的循环)。
  数据的分组运算(聚合、转换、函数应⽤等)。
Python的⾯向数组计算起源1995年,JimHugunin创建了Numeric库。在2005年,Travis Oliphant从Numeric和Numarray项⽬整合出了NumPy项⽬,进⽽所有社区都集合到了这个框架下。

NumPy之于数值计算特别重要的原因之⼀,是因为它可以⾼效处理⼤数组的数据。这是因为:NumPy是在⼀个连续的内存块中存储数据,独⽴于其他Python内置对象。
  NumPy的C语⾔编写的算法库可以操作内存,⽽不必进⾏类型检查或其它前期⼯作。
  ⽐起Python的内置序列,NumPy数组使⽤的内存更少。
NumPy可以在整个数组上执⾏复杂的计算,⽽不需要Python的for循环。要搞明⽩具体的性能差距,考察⼀个包含⼀百万整数的数组,和⼀个等价的Python列表:
import numpy as np
my_arr = np.arange(1000000)
my_list = list(range(1000000))
%time for _ in range(10): my_arr2 = my_arr * 2         # 输出:Wall time: 30 ms
%time for _ in range(10): my_list2 = [x * 2 for x in my_list]   # 输出:Wall time: 1.24 s
基于NumPy的算法要⽐纯Python快10到100倍(甚⾄更快),并且使⽤的内存更少。

 

一、NumPy的ndarray:⼀种多维数组对象
  NumPy最重要的⼀个特点就是其N维数组对象(即ndarray),该对象是⼀个快速⽽灵活的⼤数据集容器。你可以利⽤这种数组对整块数据执⾏⼀些数学运算,其语法跟标量元素之间的运算⼀样。要明⽩Python是如何利⽤与标量值类似的语法进⾏批次计算,我先引⼊NumPy,然后⽣成⼀个包含随机数据的⼩数组:
import numpy as np
data = np.random.randn(2, 3)  # Generate some random data
data               # 输出如下:
array([[-0.77451155, 0.44650751, -0.04153134],
   [ 2.94286194, -1.39901653, -0.07859503]])
data * 10           # 输出如下:,该操作不会改变原数组data,所有的元素都乘以10
array([[ -7.7451155 , 4.46507507, -0.41531343],
   [ 29.4286194 , -13.99016527, -0.78595028]])
data + data          # 输出如下:,该操作不会改变原数组data,每个元素都与⾃身相加。
array([[-1.5490231 , 0.89301501, -0.08306269],
   [ 5.88572388, -2.79803305, -0.15719006]])
ndarray是⼀个通⽤的同构数据多维容器,也就是说,其中的所有元素必须是相同类型的。
每个数组都有⼀个shape(⼀个表示各维度⼤⼩的元组)和⼀个dtype(⼀个⽤于说明数组数据类型的对象)
data.shape     # 输出:(2,3),表示2行3列
data.dtype     # 输出:dtype('float64')
下面看到的“数组”、“NumPy数组”、"ndarray"时,都指的ndarray对象。

1、创建ndarray

创建数组最简单的办法就是使⽤array函数。它接受⼀切序列型的对象(包括其他数组),然后产⽣⼀个新的含有传⼊数据的NumPy数组。以⼀个列表的转换为例:
data1 = [6, 7.5, 8, 0, 1]
arr1 = np.array(data1)     # python列表转换为Numpy数组
arr1             # 输出:array([ 6. , 7.5, 8. , 0. , 1. ])

嵌套序列(⽐如由⼀组等⻓列表组成的列表)将会被转换为⼀个多维数组:
data2 = [[1, 2, 3, 4], [5, 6, 7, 8]]
arr2 = np.array(data2)
In [24]: arr2         # 输出如下:
array([[1, 2, 3, 4],
   [5, 6, 7, 8]])
因为data2是列表的列表,NumPy数组arr2的两个维度的shape是从data2引⼊的。可以⽤属性ndimshape验证:
arr2.ndim          # 输出:2
arr2.shape          # 输出:(2, 4)
除⾮特别说明,np.array会尝试为新建的这个数组推断出⼀个较为合适的数据类型。数据类型保存在⼀个特殊的dtype对象中。⽐如说,在上⾯的两个例⼦中,我们有:
arr1.dtype          # 输出:dtype('float64')
arr2.dtype            # 输出:dtype('int32')
除np.array之外,还有⼀些函数也可以新建数组。⽐如,zerosones分别可以创建指定⻓度或形状的全0或全1数组。empty可以创建⼀个没有任何具体值的数组。要⽤这些⽅法创建多维数组,只需传⼊⼀个表示形状的元组即可:
np.zeros(10)       # 输出:array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
np.zeros((3, 6))        # 输出如下:
array([[0., 0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0., 0.]])
np.empty((2, 3, 2))     # 输出如下:
array([[[0., 0.],
   [0., 0.],
   [0., 0.]],
     [[0., 0.],
   [0., 0.],
   [0., 0.]]])
np.empty会返回的都是⼀些未初始化的垃圾值。认为返回全0是不安全的。
arange是Python内置函数range的数组版:
np.arange(15)     # 输出:array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14])
Numpy数组创建常用函数如下表4-1所示:

 

2、ndarray的数据类型
dtype(数据类型)是⼀个特殊的对象,它含有ndarray将⼀块内存解释为特定数据类型所需的信息:
arr1 = np.array([1, 2, 3], dtype=np.float64)
arr2 = np.array([1, 2, 3], dtype=np.int32)
arr1.dtype                # 输出:dtype('float64')
arr2.dtype                # 输出:dtype('int32')
dtype是NumPy灵活交互其它系统的源泉之⼀。多数情况下,它们直接映射到相应的机器表示,这使得“读写磁盘上的⼆进制数据流”以及“集成低级语⾔代码(如C、Fortran)”等⼯作变得更加简单。

数值型dtype的命名⽅式相同:⼀个类型名(如float或int),后⾯跟⼀个⽤于表示各元素位⻓的数字。
标准的双精度浮点值(即Python中的float对象)需要占⽤8字节(即64位)。因此,该类型在NumPy中就记作float64。
NumPy的数据类型如下表4-2所示:

 

 

可通过ndarray的astype⽅法明确地将⼀个数组从⼀个dtype转换成另⼀个dtype:
arr = np.array([1, 2, 3, 4, 5])
arr.dtype          # dtype('int64')
float_arr = arr.astype(np.float64)
float_arr.dtype       # dtype('float64')

如果将浮点数转换成整数,则⼩数部分将会被截取删除:
arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
arr.astype(np.int32)      # array([ 3, -1, -2, 0, 12, 10], dtype=int32)

如果某字符串数组表示的全是数字,也可以⽤astype将其转换为数值形式:
numeric_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_)
numeric_strings.astype(float)   # array([ 1.25, -9.6 , 42. ])
注意:使⽤numpy.string_类型时,⼀定要⼩⼼,因为NumPy的字符串数据是⼤⼩固定的,发⽣截取时,不会发出警告。转换过程因为某种原因⽽失败了(⽐如某个不能被转换为float64的字符串),就会引发⼀个ValueError。

数组的dtype还有另⼀个属性:将本数组的dtype设置为另一个数组的dtype
int_array = np.arange(10)
calibers = np.array([.22, .270, .357, .380, .44, .50], dtype=np.float64)
int_array.astype(calibers.dtype)           # array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.])
int_array.astype(calibers.dtype).dtype == 'float64'   # 输出:True

还可以⽤简洁的类型代码来表示dtype:

empty_uint32 = np.empty(8, dtype='u4')
empty_uint32     # 输出如下:
array([ 0, 1075314688, 0, 1075707904, 0,
   1075838976, 0, 1072693248], dtype=uint32)

调⽤astype总会创建⼀个新的数组(⼀个数据的备份),即使新的dtype与旧的dtype相同。

3、NumPy数组的运算
数组不⽤编写循环即可对数据执⾏批量运算。称为⽮量化(vectorization)。⼤⼩相等的数组之间的任何算术运算都会将运算应⽤到元素级:
arr = np.array([[1., 2., 3.], [4., 5., 6.]])   # arr = np.arange(1,7).reshape(2,3),用该方法创建同等shape的整数数组
arr * arr     # 输出如下:
array([[ 1., 4., 9.],
   [ 16., 25., 36.]])
arr - arr     # 输出如下:
array([[ 0., 0., 0.],
   [ 0., 0., 0.]])
数组与标量的算术运算会将标量值传播到各个元素:
1 / arr     # 输出如下:
array([[ 1. , 0.5 , 0.3333],
   [ 0.25 , 0.2 , 0.1667]])
arr ** 0.5   # 输出如下:
array([[ 1. , 1.4142, 1.7321],
   [ 2. , 2.2361, 2.4495]])
⼤⼩相同的数组之间的⽐较会⽣成布尔值数组:
arr2 = np.array([[0., 4., 1.], [7., 2., 12.]])
arr2 > arr     # 输出如下:元素级比较
array([[False, True, False],
   [ True, False, True]], dtype=bool)
不同⼤⼩的数组之间的运算叫做⼴播(broadcasting)

4、基本的索引和切⽚

⼀维数组很简单。从表⾯上看,它们跟Python列表的功能差不多。
arr = np.arange(10)
arr     # 输出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
arr[5:8] = 12   # 切片赋值,自动传播到整个选区
arr       # 输出:array([ 0, 1, 2, 3, 4, 12, 12, 12, 8, 9])
当你将⼀个标量值赋值给⼀个切⽚时(如arr[5:8]=12),该值会⾃动传播到整个选区。跟列表最重要的区别在于,数组切⽚是原始数组的视图。这意味着数据不会被复制,视图上的任何修改都会直接反映到源数组上。作为例⼦,先创建⼀个arr的切⽚:

arr_slice = arr[5:8]
当修改arr_slice中的值,变动也会体现在原始数组arr中:
arr_slice[1] = 12345      # 原始数组中的值也会改变
arr             # 输出:array([ 0, 1, 2, 3, 4, 12, 12345, 12, 8, 9])
切⽚[ : ]会给数组中的所有值赋值:
arr_slice[:] = 64
arr       # array([ 0, 1, 2, 3, 4, 64, 64, 64, 8, 9])
注意:如果要得到的是ndarray切⽚的⼀份副本⽽⾮视图,就需要明确地进⾏复制操作,例如 arr[5:8].copy()

对于⾼维度数组,能做的事情更多。在⼀个⼆维数组中,各索引位置上的元素不再是标量⽽是⼀维数组:
arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr2d[2]   # array([7, 8, 9])
因此,可以对各个元素进⾏递归访问,但这样需要做的事情有点多。你可以传⼊⼀个以逗号隔开的索引列表来选取单个元素。也就是说,下⾯两种⽅式是等价的:
arr2d[0][2]    # 输出:3
arr2d[0, 2]    # 输出:3
arr2d的索引方式中,轴0代表行,轴2代表列

在多维数组中,如果省略了后⾯的索引,则返回对象会是⼀个维度低⼀点的ndarray(它含有⾼⼀级维度上的所有数据)。因此,在2×2×3数组arr3d中:
arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
arr3d[0]是⼀个2×3数组:
arr3d[0]    # 输出如下:
array([[1, 2, 3],
   [4, 5, 6]])

标量值和数组都可以被赋值给 arr3d[0]:
old_values = arr3d[0].copy()
arr3d[0] = 42
arr3d     # 输出如下:
array([[[42, 42, 42],
    [42, 42, 42]],
   [[ 7, 8, 9],
    [10, 11, 12]]])
arr3d[0] = old_values    # 赋回原值
arr3d          # 输出如下:
array([[[ 1, 2, 3],
    [ 4, 5, 6]],
   [[ 7, 8, 9],
   [10, 11, 12]]])

相似的,arr3d[1,0]可以访问索引以(1,0)开头的那些值(以⼀维数组的形式返回):
arr3d[1, 0]      # 输出:array([7, 8, 9])
虽然是⽤两步进⾏索引的,表达式是相同的:
x = arr3d[1]     # 先赋值给x
x          # 输出如下:
array([[ 7, 8, 9],
   [10, 11, 12]])
x[0]         # 输出:array([7, 8, 9])
在上⾯所有这些选取数组⼦集的例⼦中,返回的数组都是视图。

 

5、切⽚索引

ndarray的切⽚语法跟Python列表这样的⼀维对象差不多:对于之前的⼆维数组arr2d,其切⽚⽅式稍显不同
arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr2d[:2]     # 输出如下:
array([[1, 2, 3],
   [4, 5, 6]])
切⽚是沿着⼀个轴向选取元素的。表达式arr2d[:2]可以被认为是“选取arr2d的前两⾏”。
可以⼀次传⼊多个切⽚,就像传⼊多个索引那样:
arr2d[:2, 1:]   # 输出如下:
array([[2, 3],
   [5, 6]])
像这样进⾏切⽚时,只能得到相同维数的数组视图。通过将整数索引和切⽚混合,可以得到低维度的切⽚。
选取第⼆⾏的前两列:
arr2d[1, :2]     # 输出:array([4, 5])
相似的,还可以选择第三列的前两⾏:
arr2d[:2, 2]     # 输出:array([3, 6])

只有冒号”表示选取整个轴,因此你可以像下⾯这样只对⾼维轴进⾏切⽚:
arr2d[:, :1]      # 输出如下:
array([[1],
   [4],
   [7]])

对切⽚表达式的赋值操作也会被扩散到整个选区:
arr2d[:2, 1:] = 0
arr2d        # 输出如下:
array([[1, 0, 0],
   [4, 0, 0],
   [7, 8, 9]])

 

6、布尔型索引,看一个例子:
假设我们有⼀个⽤于存储数据的数组以及⼀个存储姓名的数组(含有重复项)。先使⽤numpy.random中的randn函数⽣成⼀些正态分布的随机数据
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
data = np.random.randn(7, 4)
names     # 输出:array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='<U4')
data      # 输出如下:
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"的⽐较运算将会产⽣⼀个布尔型数组:
names == 'Bob'     # 输出:array([ True, False, False, True, False, False, False])
这个布尔型数组可⽤于数组索引:
data[names == 'Bob']   # 输出如下:(如果要取相应的列,可这样做:data[:,names == 'Bob'],前提是长度要一致)
array([[ 0.0929, 0.2817, 0.769 , 1.2464],
   [ 1.669 , -0.4386, -0.5397, 0.477 ]])

布尔型数组的⻓度必须跟被索引的轴⻓度⼀致。此外,还可以将布尔型数组跟切⽚、整数(或整数序列)混合使⽤。
注意:如果布尔型数组的⻓度不对,布尔型选择就会出错,因此⼀定要⼩⼼。

下⾯的例⼦,选取了names == 'Bob'的⾏,并索引了列:布尔型数组与切片混合使用
data[names == 'Bob', 2:]  # 输出如下:
array([[ 0.769 , 1.2464],
   [-0.5397, 0.477 ]])
data[names == 'Bob', 3]   # 输出:array([ 1.2464, 0.477 ])

要选择除"Bob"以外的其他值,既可以使⽤不等于符号(!=),也可以通过(~)对条件进⾏否定
names != 'Bob'       # 输出:array([False, True, True, False, True, True, True])
data[~(names == 'Bob')]   # 输出如下:
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]])
(~)操作符⽤来反转条件很好⽤:
cond = names == 'Bob'
data[~cond]     # 输出与前面的一样

选取这三个名字中的两个需要组合应⽤多个布尔条件,使⽤&(与)、|(或)之类的布尔算术运算符即可:
mask = (names == 'Bob') | (names == 'Will')
mask       # 输出:array([ True, False, True, True, True, False, False])
data[mask]     # 输出如下:
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,我们只需:
data[data < 0] = 0
data       # 输出如下:
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. ]])

通过⼀维布尔数组设置整⾏或列的值也很简单:
data[names != 'Joe'] = 7
data         # 输出如下:
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. ]])

 

7、花式索引(复杂索引,Fancy indexing)

花式索引(Fancy indexing)是⼀个NumPy术语,它指的是利⽤整数数组进⾏索引。假设有⼀个8×4数组:
arr = np.empty((8, 4))
for i in range(8):
  arr[i] = i
arr     # 输出如下:
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即可:
arr[[4, 3, 0, 6]]
array([[4., 4., 4., 4.],
  [3., 3., 3., 3.],
  [0., 0., 0., 0.],
  [6., 6., 6., 6.]])
这段代码确实达到我们的要求了!使⽤负数索引将会从末尾开始选取⾏
arr[[-3, -5, -7]]
array([[5., 5., 5., 5.],
  [3., 3., 3., 3.],
  [1., 1., 1., 1.]])
⼀次传⼊多个索引数组会有⼀点特别。它返回的是⼀个⼀维数组,其中的元素对应各个索引元组
arr = np.arange(32).reshape((8, 4))
arr     # 输出如下:
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]])
arr[[1, 5, 7, 2], [0, 3, 1, 2]]     #输出:array([ 4, 23, 29, 10])
最终选出的是元素(1,0)、(5,3)、(7,1)和(2,2)。⽆论数组是多少维的,花式索引总是⼀维的。
花式索引的⾏为可能会跟预期不⼀样,选取矩阵的⾏列⼦集应该是矩形区域的形式才对。下⾯是得到该结果的⼀个办法:
arr[[1, 5, 7, 2]][:, [0, 3, 1, 2]]   # 输出如下:
array([[ 4, 7, 5, 6],      # 该行选取的元素是:(1, 0),(1, 3),(1, 1),(1, 2)
  [20, 23, 21, 22],     # 该行选取的元素是:(5, 0),(5, 3),(5, 1),(5, 2),以此类推
  [28, 31, 29, 30],
  [ 8, 11, 9, 10]])
花式索引跟切⽚不⼀样,它总是将数据复制到新数组中。

 

8、数组转置和轴对换
转置是重塑的⼀种特殊形式,它返回的是源数据的视图(不会进⾏任何复制操作)。
数组不仅有transpose⽅法,还有⼀个特殊的T属性
arr = np.arange(15).reshape((3,5))
arr         # 输出如下:
array([[ 0, 1, 2, 3, 4],
  [ 5, 6, 7, 8, 9],
  [10, 11, 12, 13, 14]])
arr.T       # 输出如下:使用T属性转置
array([[ 0, 5, 10],
  [ 1, 6, 11],
  [ 2, 7, 12],
  [ 3, 8, 13],
  [ 4, 9, 14]])
在进⾏矩阵计算时,经常需要⽤到该操作,⽐如利⽤np.dot计算矩阵内积
arr = np.random.randn(6, 3)
arr       # 输出如下:
array([[-2.13468628, 1.06582408, -0.72489158],
  [-1.17032574, -0.82376495, 0.07941442],
  [-0.29019922, -1.98512296, 0.12906019],
  [ 0.14798076, 0.9989999 , -0.83814468],
  [-1.61693841, -0.23485318, -0.38723287],
  [ 0.09921124, -1.21635853, -2.44492694]])
np.dot(arr.T, arr)   # 输出如下:
array([[ 8.65699441, -0.32814613, 1.67656039],
  [-0.32814613, 8.28796773, 1.1333181 ],
  [ 1.67656039, 1.1333181 , 7.37853451]])

对于⾼维数组,transpose需要得到⼀个由轴编号组成的元组才能对这些轴进⾏转置
arr = np.arange(16).reshape((2, 2, 4))
arr      # 输出如下:
array([[[ 0, 1, 2, 3],
  [ 4, 5, 6, 7]],
  [[ 8, 9, 10, 11],
  [12, 13, 14, 15]]])
arr.transpose((1, 0, 2))   # 输出如下:
array([[[ 0, 1, 2, 3],
  [ 8, 9, 10, 11]],
  [[ 4, 5, 6, 7],
  [12, 13, 14, 15]]])
这⾥,第⼀个轴被换成了第⼆个,第⼆个轴被换成了第⼀个,最后⼀个轴不变。

简单的转置可以使⽤.T,它其实就是进⾏轴对换⽽已。ndarray还有⼀个swapaxes⽅法,它需要接受⼀对轴编号:
arr.swapaxes(1, 2)   # 输出如下:
array([[[ 0, 4],
   [ 1, 5],
   [ 2, 6],
   [ 3, 7]],

   [[ 8, 12],
   [ 9, 13],
   [10, 14],
   [11, 15]]])
swapaxes也是返回源数据的视图(不会进⾏任何复制操作)。

 

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

 

通⽤函数(即ufunc)是⼀种对ndarray中的数据执⾏元素级运算的函数。你可以将其看做简单函数(接受⼀个或多个标量值,并产⽣⼀个或多个标量值)的⽮量化包装器。许多ufunc都是简单的元素级变体,如sqrt和exp
arr = np.arange(10)
arr       # 输出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
np.sqrt(arr)   # 输出:(相当于arr ** 0.5)
array([0. , 1. , 1.41421356, 1.73205081, 2. ,
    2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
这些都是⼀元(unary)ufunc。另外⼀些(如add或maximum)接受2个数组(因此也叫⼆元(binary)ufunc),并返回⼀个结果数组:
x = np.random.randn(8)
y = np.random.randn(8)
x     # 输出如下:
array([ 1.20051688, -1.65710862, -0.53963283, -1.21187257, 0.63507697,
  -0.37209111, -0.93943269, -1.00426802])
y     # 输出如下:
array([ 1.05467508, 1.4219286 , -0.92969225, -1.04814943, -0.82260389,
  0.91427391, 0.2504812 , 0.97563071])
np.maximum(x, y)   # 输出如下:实际比较的是(x[0], y[0]), (x[1], y[1])...
array([ 1.20051688, 1.4219286 , -0.53963283, -1.04814943, 0.63507697,
  0.91427391, 0.2504812 , 0.97563071])
numpy.maximum计算了x和y中元素级别最⼤的元素。

modf就是Python内置函数divmod的⽮量化版本,它会返回浮点数数组的⼩数和整数部分
arr = np.random.randn(7) * 5
arr     # 原始数据如下:
array([-3.52389072, -3.4832977 , 2.56810837, 1.2339455 , 6.10703696,
  7.61798681, 3.05483746])
remainder, whole_part = np.modf(arr)
remainder      # 输出小数部分
array([-0.52389072, -0.4832977 , 0.56810837, 0.2339455 , 0.10703696,
  0.61798681, 0.05483746])
whole_part     # 输出整数部分
array([-3., -3., 2., 1., 6., 7., 3.])

Ufuncs接受out选项参数,可以让它们在数组的原地进⾏操作:
arr       # 输出如下:
array([ 3.43072847, 10.82286699, -3.63291583, 4.07716536, -1.60187118, -1.24558177, -5.24618548])
np.sqrt(arr)   # 输出如下:
array([1.85222258, 3.28981261, nan, 2.01919919, nan, nan, nan])
arr       # 输出如下:
array([1.85222258, 3.28981261, nan, 2.01919919, nan, nan, nan])

一元ufunc(一元通用函数)如下表4-3所示:

二元ufunc(二元通用函数)如下表4-4所示:

 

三、利⽤数组进⾏数据处理


⽤数组表达式代替循环的做法,通常被称为⽮量化。⽮量化数组运算要⽐等价的纯Python⽅式快上⼀两个数量级(甚⾄更多),尤其是各种数值计算。在后⾯将介绍⼴播,这是⼀种针对⽮量化计算的强⼤⼿段。看一个例子:假设我们想要在⼀组值(⽹格型)上计算函数sqrt(x^2+y^2)。

np.meshgrid函数接受两个⼀维数组,并产⽣两个⼆维矩阵(对应于两个数组中所有的(x,y)对):
points = np.arange(-5, 5, 0.01)    # 1000 equally spaced points
xs, ys = np.meshgrid(points, points)
ys               # 输出如下:
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]])
现在,把这两个数组当做两个浮点数那样编写表达式即可:
z = np.sqrt(xs ** 2 + ys ** 2)
z           # 输出如下:
array([[7.07106781, 7.06400028, 7.05693985, ..., 7.04988652, 7.05693985,
    7.06400028],
  [7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
    7.05692568],
  [7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
    7.04985815],
  ...,
  [7.04988652, 7.04279774, 7.03571603, ..., 7.0286414 , 7.03571603,
    7.04279774],
  [7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
    7.04985815],
  [7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
    7.05692568]])
%matplotlib
import matplotlib.pyplot as plt
plt.imshow(z, cmap=plt.cm.gray); plt.colorbar()   # 输出:<matplotlib.colorbar.Colorbar at 0x1f488261e80>
plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")    # 输出:Text(0.5,1,'Image plot of $\\sqrt{x^2 + y^2}$ for a grid of values')

 

1、将条件逻辑表述为数组运算
numpy.where函数是三元表达式x if condition else y的⽮量化版本。假设我们有⼀个布尔数组和两个值数组:
xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
cond = np.array([True, False, True, True, False])
根据cond中的值选取xarr和yarr的值:当cond中的值为True时,选取xarr的值,否则从yarr中选取。列表推导式的写法应该如下所示:
result = [(x if c else y) for x, y, c in zip(xarr, yarr, cond)]
result     # 输出:[1.1, 2.2, 1.3, 1.4, 2.5]
第⼀,它对⼤数组的处理速度不是很快(因为所有⼯作都是由纯Python完成的)。
第⼆,⽆法⽤于多维数组。若使⽤np.where,则可以将该功能写得⾮常简洁:
result1 = np.where(cond, xarr, yarr)
result1     # 输出:array([1.1, 2.2, 1.3, 1.4, 2.5])
np.where的第⼆个和第三个参数不必是数组,它们都可以是标量值。在数据分析⼯作中,where通常⽤于根据另⼀个数组⽽产⽣⼀个新的数组。

假设有⼀个由随机数据组成的矩阵,你希望将所有正值替换为2,将所有负值替换为-2。若利⽤np.where,则会⾮常简单:
arr = np.random.randn(4,4)
arr       # 输出如下:
array([[ 0.80426609, 1.0974934 , 1.52641017, 0.26308459],
  [ 0.39494499, -0.09829609, 2.93032116, 0.77040912],
  [-1.35818382, 0.78914328, 1.38073343, 0.98920804],
  [ 0.93070469, 0.32393877, -0.76748938, -0.94859212]])
arr > 0     # 输出如下:
array([[ True, True, True, True],
  [ True, False, True, True],
  [False, True, True, True],
  [ True, True, False, False]])
np.where(arr > 0, 2, -2)   # 输出如下:
array([[ 2, 2, 2, 2],
  [ 2, -2, 2, 2],
  [-2, 2, 2, 2],
  [ 2, 2, -2, -2]])
使⽤np.where,可以将标量和数组结合起来。例如,我可⽤常数2替换arr中所有正的值:
np.where(arr > 0, 2, arr)     # 大于0的值设置为2
array([[ 2. , 2. , 2. , 2. ],
  [ 2. , -0.09829609, 2. , 2. ],
  [-1.35818382, 2. , 2. , 2. ],
  [ 2. , 2. , -0.76748938, -0.94859212]])
传递给where的数组⼤⼩可以不相等,甚⾄可以是标量值。

 

2、数学和统计⽅法

可以通过数组上的⼀组数学函数对整个数组或某个轴向的数据进⾏统计计算。sum、mean以及标准差std等聚合计算(aggregation,通常叫做约简(reduction))既可以当做数组的实例⽅法调⽤,也可以当做顶级NumPy函数使⽤。下面例子,先生成一些正态分布随机数据,然后做了聚类统计:
arr = np.random.randn(5,4)
arr        # 输出如下:
array([[-0.1133632 , 0.58949215, 0.275524 , -0.3855296 ],
  [ 2.55070219, 0.63003176, -0.34253814, -0.01309244],
  [ 1.04023817, -1.9999175 , -0.40581791, -0.36883092],
  [-0.71064938, -1.40922625, -1.18840513, -0.92986452],
  [-0.05223665, 0.5629013 , 0.69139223, -1.45967786]])
arr.mean()      # 输出:-0.15194338563310547
np.mean(arr)      # 输出:-0.15194338563310547
arr.sum()      # 输出:-3.0388677126621095

mean和sum这类的函数可以接受⼀个axis选项参数,⽤于计算该轴向上的统计值,最终结果是⼀个少⼀维的数组:
arr.mean(axis=1)   # 输出:array([ 0.09153084, 0.70627584, -0.43358204, -1.05953632, -0.06440525])
arr.sum(axis=0)    # 输出:array([ 2.71469112, -1.62671854, -0.96984495, -3.15699534])
这⾥,arr.mean(1)是“计算⾏的平均值”,arr.sum(0)是“计算每列的和”。

其他如cumsum和cumprod之类的⽅法则不聚合,⽽是产⽣⼀个由中间结果组成的数组:
arr = np.array([0, 1, 2, 3, 4, 5, 6, 7])
arr.cumsum()     # 输出:array([ 0, 1, 3, 6, 10, 15, 21, 28], dtype=int32),累加

在多维数组中,累加函数(如cumsum)返回的是同样⼤⼩的数组,但是会根据每个低维的切⽚沿着标记轴计算部分聚类:
arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
arr.cumsum(axis=0)      # 纵向累计求和,输出如下:
array([[ 0, 1, 2],
  [ 3, 5, 7],
  [ 9, 12, 15]], dtype=int32)
arr.cumprod(axis=1)      # 横向累计求积,输出如下:
array([[ 0, 0, 0],
  [ 3, 12, 60],
  [ 6, 42, 336]], dtype=int32)

全部的基本数组统计方法如下表4-5所示:

 

3、⽤于布尔型数组的⽅法
在上⾯这些⽅法中,布尔值会被强制转换为1(True)和0(False)。因此,sum经常被⽤来对布尔型数组中的True值计数:
arr = np.random.randn(100)
(arr > 0).sum()     # 输出:一个随机数,大于0的个数
另外还有两个⽅法any和all,它们对布尔型数组⾮常有⽤。any⽤于测试数组中是否存在⼀个或多个True,⽽all则检查数组中所有值是否都是True:
bools = np.array([False, False, True, False])
bools.any()       # 输出:True
bools.all()        # 输出:False
这两个⽅法也能⽤于⾮布尔型数组,所有⾮0元素将会被当做True。

 

4、排序

NumPy数组也可以通过sort⽅法就地排序:要改变原始数组
arr = np.random.randn(6)
arr          # 输出如下:
array([ 1.37750254, 0.36255115, -0.87588226, -0.80512566, -1.33751372, 0.18521148])
arr.sort()        # 就地排序
arr            # 输出如下:
array([-1.33751372, -0.87588226, -0.80512566, 0.18521148, 0.36255115, 1.37750254])

多维数组可以在任何⼀个轴向上进⾏排序,只需将轴编号传给sort即可:
arr = np.random.randn(5, 3)
arr        # 输出如下:
array([[-0.70497864, 0.75979962, -0.04453997],
  [ 1.44800768, -1.67338973, 1.1078619 ],
  [ 0.40954131, 1.03518218, -0.49360008],
  [ 0.30501957, 1.29319448, 1.02856608],
  [-1.16897745, -1.51113507, 0.9717796 ]])
arr.sort(1)      # 对每一行的数进行排序
arr
array([[-0.70497864, -0.04453997, 0.75979962],
  [-1.67338973, 1.1078619 , 1.44800768],
  [-0.49360008, 0.40954131, 1.03518218],
  [ 0.30501957, 1.02856608, 1.29319448],
  [-1.51113507, -1.16897745, 0.9717796 ]])
顶级⽅法np.sort返回的是数组的已排序副本,⽽就地排序则会修改数组本身。计算数组分位数最简单的办法是对其进⾏排序,然后选取特定位置的值:
large_arr = np.random.randn(1000)     # 生成1000个随机数
large_arr.sort()              # 从小到大排序
large_arr[int(0.05 * len(large_arr))]      # 5% 分位数(quantile)
# 输出:-1.5311513550102103

 

5、唯⼀化以及其它的集合逻辑
NumPy提供了⼀些针对⼀维ndarray的基本集合运算。最常⽤的可能要数np.unique了,
它⽤于找出数组中的唯⼀值并返回已排序的结果:
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
# 去掉重复元素,并按字母排序,
np.unique(names) # 输出:array(['Bob', 'Joe', 'Will'], dtype='<U4')
ints = np.array([3, 3, 3, 2, 2, 1, 1, 4, 4]) # 对数字排序
np.unique(ints) # 输出:array([1, 2, 3, 4])
拿np.unique跟等价的纯Python代码来对⽐⼀下:
sorted(set(names)) # 输出:['Bob', 'Joe', 'Will']

另⼀个函数np.in1d⽤于测试⼀个数组中的值在另⼀个数组中的成员资格,返回⼀个布尔型数组:
values = np.array([6, 0, 0, 3, 2, 5, 6])
np.in1d(values, [2, 3, 6]) # 输出:array([ True, False, False, True, True, False, True])

NumPy数组的集合运算(表4-6):


四、⽤于数组的⽂件输⼊输出


NumPy能够读写磁盘上的⽂本数据或⼆进制数据。先讨论NumPy的内置⼆进制格式,因为更多的⽤户会使⽤pandas或其它⼯具加载⽂本或表格数据。

np.save和np.load是读写磁盘数组数据的两个主要函数。默认情况下,数组是以未压缩的原始⼆进制格式保存在扩展名为.npy的⽂件中的:例如:
arr = np.arange(10)
np.save('some_array', arr)
如果⽂件路径末尾没有扩展名.npy,则该扩展名会被⾃动加上。然后就可以通过np.load读取磁盘上的数组:
np.load('some_array.npy')   # 输出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

通过np.savez可以将多个数组保存到⼀个未压缩⽂件中,将数组以关键字参数的形式传⼊即可:
np.savez('array_archive.npz', a=arr, b=arr)
加载.npz⽂件时,你会得到⼀个类似字典的对象,该对象会对各个数组进⾏延迟加载:
arch = np.load('array_archive.npz')
arch['b']      # 输出:array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

如果数据压缩的很好,就可以使⽤numpy.savez_compressed
np.savez_compressed('arrays_compressed.npz', a=arr, b=arr)

 

五、线性代数


线性代数(如矩阵乘法、矩阵分解、⾏列式以及其他⽅阵数学等)是任何数组库的重要组成部分。MATLAB通过 对两个⼆维数组相乘得到的是⼀个元素级的积,⽽不是⼀个矩阵点积。NumPy提供了⼀个⽤于矩阵乘法的dot函数(既是⼀个数组⽅法也是numpy命名空间中的⼀个函数):
x = np.array([[1., 2., 3.], [4., 5., 6.]])
y = np.array([[6., 23.], [-1, 7], [8, 9]])
x.dot(y)        # 输出如下:得到的结果是行数等于x的行,列数等于y的列
array([[ 28., 64.],
  [ 67., 181.]])
x.dot(y)等价于np.dot(x, y):
np.dot(x, y)     # 输出同上

⼀个⼆维数组跟⼀个⼤⼩合适的⼀维数组的矩阵点积运算之后将会得到⼀个⼀维数组:
np.dot(x, np.ones(3))      # 输出:array([ 6., 15.])
x数组的每个元素乘以第二个数组的每个元素后相加的结果

@符(类似Python 3.5)也可以⽤作中缀运算符,进⾏矩阵乘法
x @ np.ones(3)        # 输出:array([ 6., 15.])

numpy.linalg中有⼀组标准的矩阵分解运算以及诸如求逆和⾏列式之类的东⻄。它们跟MATLAB和R等语⾔所使⽤的是相同的⾏业标准线性代数库,如BLAS、LAPACK、Intel MKL(MathKernel Library,可能有,取决于你的NumPy版本)等:
from numpy.linalg import inv, qr
X = np.random.randn(5, 5)
mat = X.T.dot(X)
inv(mat)          # 输出如下:
array([[ 0.28133188, -0.0977005 , -0.47041506, 0.01876554, 0.32981926],
  [-0.0977005 , 0.31194322, 0.15049957, -0.27441188, -0.58712427],
  [-0.47041506, 0.15049957, 1.1281112 , -0.07145714, -0.59076183],
  [ 0.01876554, -0.27441188, -0.07145714, 0.56844553, 0.74076293],
  [ 0.32981926, -0.58712427, -0.59076183, 0.74076293, 1.68251546]])
mat.dot(inv(mat))     # 输出如下:
array([[ 1.00000000e+00, -7.03456949e-17, -2.26629202e-16, -3.58528577e-17, 2.46988638e-16],
  [-4.69384500e-16, 1.00000000e+00, -1.58143616e-17, -3.57503814e-16, -5.44789041e-16],
  [-2.31542198e-16, 1.52236389e-16, 1.00000000e+00,-1.57566693e-16, -4.74249967e-16],
  [ 2.37545108e-16, 1.70674192e-16, 1.02588110e-16, 1.00000000e+00, -5.48276527e-16],
  [-2.69371029e-16, 1.04852032e-16, -5.72195720e-17, -7.39959808e-17, 1.00000000e+00]])
q, r = qr(mat)
r         # 输出如下:
array([[-16.80324568, -0.94736212, -5.62018513, -5.74117222, 3.61353809],
   [ 0., -10.27117765, -0.92342637, 0.583159 , -4.36651116],
   [ 0.  , 0. , -1.83415284, 4.48040783, -2.86906256],
   [ 0. , 0. , 0. , -3.01189889, 1.35175951],
   [ 0. , 0. , 0. , 0. , 0.48899415]])
表达式X.T.dot(X)计算X和它的转置X.T的点积。

常用的numpy.linalg函数(线性代数函数)表4-7:

 

六、伪随机数⽣成


numpy.random模块对Python内置的random进⾏了补充,增加了⼀些⽤于⾼效⽣成多种概率分布的样本值的函数。例如,可⽤normal来得到⼀个标准正态分布的4×4样本数组:
samples = np.random.normal(size=(4, 4))
samples      # 输出如下:
array([[-1.22523734, 1.0020195 , 0.62450895, -0.30984473],
  [ 0.17302941, -1.1195615 , -1.4692486 , 0.47059667],
  [ 1.12656905, 0.66808849, -0.69090706, 0.13304373],
  [ 0.46783217, -1.33587931, 1.19216589, 0.13784556]])
⽽Python内置的random模块则只能⼀次⽣成⼀个样本值。从下⾯的测试结果中可以看出,如果需要产⽣⼤量样本值,numpy.random快了不⽌⼀个数量级:
from random import normalvariate
N = 1000000
%timeit samples = [normalvariate(0, 1) for _ in range(N)]   # 输出如下:
886 ms ± 5.22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.random.normal(size=N)              # 输出如下:
40.4 ms ± 229 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

说它是伪随机数,是因为它们都是通过算法基于随机数⽣成器种⼦,在确定性的条件下⽣成的
你可以⽤NumPy的np.random.seed更改随机数⽣成种⼦:
np.random.seed(1234)

numpy.random的数据⽣成函数使⽤了全局的随机种⼦。要避免全局状态,可以使⽤numpy.random.RandomState,创建⼀个与其它隔离的随机数⽣成器:
rng = np.random.RandomState(1234)
rng.randn(10)             # 输出如下:
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873,
  0.88716294, 0.85958841, -0.6365235 , 0.01569637, -2.24268495])

下面列出了numpy.random中的部分函数(表4-8):

 

七、示例:随机漫步


通过模拟随机漫步来说明如何运⽤数组运算。⼀个简单的随机漫步的例⼦:从0开始,步⻓1和-1出现的概率相等。
下⾯是⼀个通过内置的random模块以纯Python的⽅式实现1000步的随机漫步:
import random
position = 0
walk = [ position ]
steps = 1000
for i in range(steps):
  step = 1 if random.randint(0, 1) else -1      # random.randint(0, 1)随机选取0或1
  position += step
  walk.append(position)
根据前100个随机漫步值⽣成的折线图:
%matplotlib
import matplotlib.pyplot as plt
plt.plot(walk[:100])    # 输出如下图照片所示

从输出照片可以看出,这其实就是随机漫步中各步的累计和,可以⽤⼀个数组运算来实现。⽤np.random模块⼀次性随机产⽣1000个“掷硬币”结果(即两个数中任选⼀个),将其分别设置为1或 -1,然后计算累计和:
nsteps = 1000
draws = np.random.randint(0, 2, size=nsteps)
steps = np.where(draws > 0, 1, -1)
walk = steps.cumsum()    # cumsum()方法是累加

有了这些数据就可以沿着漫步路径做⼀些统计⼯作,⽐如求取最⼤值和最⼩值:
walk.min()      # 输出:-9,每次结果不一样
walk.max()     # 输出:60

现在来看⼀个复杂点的统计任务——⾸次穿越时间,即随机漫步过程中第⼀次到达某个特定值的时间。假设想要知道本次随机漫步需要多久才能距离初始0点⾄少10步远(任⼀⽅向均可)。np.abs(walk)>=10可以得到⼀个布尔型数组,它表示的是距离是否达到或超过10,⽽我们想要知道的是第⼀个10或 -10的索引。可以⽤argmax来解决这个问题,它返回的是该布尔型数组第⼀个最⼤值的索引(True就是最⼤值):
(np.abs(walk) >= 10).argmax()     # 输出:20
使⽤argmax并不是很⾼效,它会对数组进⾏完全扫描。在本例中,只要发现⼀个True,就知道它是个最⼤值。

 

1、⼀次模拟多个随机漫步
模拟多个随机漫步过程(⽐如5000个),只需对上⾯的代码做⼀点点修改即可⽣成所有的随机漫步过程。只要给numpy.random的函数传⼊⼀个⼆元元组就可以产⽣⼀个⼆维数组,然后我们就可以⼀次性计算5000个随机漫步过程(⼀⾏⼀个)的累计和了:
nwalks = 5000
nsteps = 1000
draws = np.random.randint(0, 2, size=(nwalks, nsteps))    # 在[0,1]之间选取,size参数指定要生成的矩阵行数和列数
steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(1)       # 输出如下:
array([[ 1, 0, 1, ..., 42, 43, 44],
  [ -1, -2, -1, ..., 40, 41, 42],
  [ -1, -2, -1, ..., -42, -43, -44],
  ...,
  [ -1, 0, 1, ..., 66, 65, 64],
  [ -1, -2, -1, ..., 0, -1, 0],
  [ -1, -2, -1, ..., 34, 35, 34]], dtype=int32)
现在,来计算所有随机漫步过程的最⼤值和最⼩值:
walks.max()       # 输出:118
walks.min()       # 输出:-130
得到这些数据之后,来计算30或 -30的最⼩穿越时间。这⾥稍微复杂些,因为不是5000个过程都到达了30。可以⽤any⽅法来对此进⾏检查:
hits30 = (np.abs(walks) >= 30).any(1)  # 返回每一行的“满足条件的第一个True,均不满足条件就返回False”,组成的一维数组
hits30         # array([ True, True, True, ..., True, False, True])
hits30.sum()      # 输出:3381,达到绝对值30次以上的个数
然后利⽤这个布尔型数组选出那些穿越了30(绝对值)的随机漫步(⾏),并调⽤argmax在轴1上获取穿越时间:
crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
crossing_times.mean()   # 输出:502.03756285122745

尝试⽤其他分布⽅式得到漫步数据。只需使⽤不同的随机数⽣成函数即可,如normal⽤于⽣成指定均值和标准差的正态分布数据:
steps = np.random.normal(loc=0, scale=0.25, size=(nwalks, nsteps))

posted @ 2018-12-05 15:48  远方那一抹云  阅读(675)  评论(0编辑  收藏  举报