第四章 NumPy基础：数组与向量化计算

  • NumPy，是Numerical Python的简称，是目前Python数值计算中最重要的基础包，其数组对象作为数据交换的通用语
  • 主要内容
  • ndarray，一种高效多维数组，提供了基于数组的便捷算术操作以及灵活的广播功能。
  • 对所有数据进行快速的矩阵计算，而无须编写循环程序。
  • 对硬盘中数组数据进行读写的工具，并对内存映射文件进行操作。
  • 线性代数、随机数生成以及傅里叶变换功能。
  • 用于连接NumPy到C、C++和FORTRAN语言类库的C语言API。
  • NumPy处理数组数据非常有效
  • NumPy在内部将数据存储在连续的内存块上
  • NumPy可以针对全量数组进行复杂计算而不需要写Python循环

4.1 NumPy ndarray：多维数组对象

  • ndarray是一个通用的多维同类数据容器，包含的每一个元素均为相同类型
  • shape属性：数组每一维度的数量
  • dtype属性：数据的类型

生成ndarray

  • array函数：接收任意的序列型对象，也包括其他数组
import numpy as np
data1 = [6,7.5,8,0,1]
arr1 = np.array(data1)
arr1  # array([6. , 7.5, 8. , 0. , 1. ])
data2 = [[1,2,3,4],[5,6,7,8]]
arr2 = np.array(data2)
arr2
"""
array([[1, 2, 3, 4],
       [5, 6, 7, 8]])
"""
arr2.ndim # 2
arr2.shape # (2, 4)

除非显式地指定，np.array会自动推断生成数组的数据类型，数据类型被存储在特殊元数据dtype中

  • 其他方式创建数组
  • zeros一次性创造全0数组
  • ones一次性创造全1数组
  • empty创建一个没有初始化数值的数组
  • 不要用np.empty来生成一个全0数组，有时候可能会返回未初始化的垃圾数值
  • 创建高维数组：为shape传递一个元组
  • arange是Python内建函数range的数组版
np.arange(15)
# array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])
  • 由于NumPy专注于数值计算，默认数据类型是float64
  • 数组生成函数

ndarray的数据类型

  • 数据类型dtype：元数据（表示数据的数据），一个特殊的对象，表示ndarray为某一种类型数据所申明的内存块信息
  • NumPy数据类型（需要在内存或硬盘上做更深入探索时需要了解）
  • 使用astype方法转换数据类型

使用astype会生成一个新的数组，与原来类型相同也会

  • 使用astype方法显示转换
arr = np.array([1,2,3,4,5])
arr.dtype # dtype('int32')
float_arr = arr.astype(np.float64)
float_arr.dtype # dtype('float64')
  • 使用另一个数组的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.])
  • 使用类型代码传入数据类型
empty_uint32 = np.empty(8,dtype='u4')
  • NumPy修正numpy.string_类型作字符串的大小或删除输入时，可能不发生警告，Pandas在处理非数值数据时有更直观的开箱型操作

NumPy数组算术

  • NumPy数组具有向量化特性，进行批量操作无须任何for循环。两个等尺寸数组之间的操作都应用元素操作的方式。
  • 同尺寸数组之间比较，不同尺寸间的数组会用到广播机制
arr = np.array([[1,2,3],[4,5,6]])
arr2 = np.array([[0,4,1],[7,2,12]])
arr2 > arr
"""
array([[False,  True, False],
       [ True, False,  True]])
"""

基础索引与切片

  • 传数值给数组的切片，数值会被传递到整个切片
  • 区别于Python的内建列表，数组的切片是原数组的视图，数据不是被复制，修改会体现到原数组上。（Numpy被设计处理非常大的数组，如果总是复制会引起很多内存问题）
  • 想要拷贝不是数组视图：使用copy()方法，arr[5:8].copy()
  • 获取元素
arr2d[0][2] # 3
arr2d[0,2] # 3
  • 数组子集选择返回的都是视图
  • 数组的切片索引
  • 多组切片
  • 单独一个冒号表示选择整个轴上的数组
  • 对切片表达式赋值时，整个切片都会重新赋值

布尔索引

  • numpy.random中的randn函数生成随机正态分布数据
  • 在索引数组时可以传入布尔值数组data[names == ‘Bob’]
  • 布尔值数组的长度必须和数组轴索引长度一致。你甚至还可以用切片或整数值（或整数值的序列，后续会介绍）对布尔值数组进行混合和匹配
  • 当布尔值数组的长度不正确时，不会报错，使用该特性要小心。
  • 选择name == ‘Bob’的行，并索引各列
data[names == 'Bob',2:]
data[names == 'Bob',3]
  • 选择除了‘Bob’以外的其他数据，可以使用!=表达式，或者在表达式前用~取反
name != 'Bob'
data[~(name == 'Bob')]
# ～符号可以在你想要对一个通用条件进行取反时使用
cond = names == 'Bob'
data[~cond]
  • 对多个布尔值条件进行联合，需要用&和|
  • 用布尔值索引选择数据时，会生成数据的拷贝。
  • Python的关键字and和or对布尔值数组并没有用，请使用&（and）和|（or）来代替。

神奇索引

神奇索引：NumPy中用于描述使用整数数组进行数据索引

  • 选出符合特定顺序的子集，传递一个包含所需顺序的列表或数组来完成；使用负索引会从尾部进行
  • 传递多个索引数组时会根据每个索引元组对应的元素选出一个一维数组
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],
       [20, 23, 21, 22],
       [28, 31, 29, 30],
       [ 8, 11,  9, 10]])
"""

数组转置和换轴

  • 转置是一种特殊的数据重组形式，返回底层数据视图，无需复制
  • 数组有transpose方法，也有特殊的T属性
  • 计算矩阵内积np.dot
arr = np.random.randn(6,3)
arr
"""
array([[ 2.87744336,  0.42238688, -0.9002813 ],
       [-0.01005598, -0.25550216, -0.44153358],
       [-2.12971346,  0.03663679,  1.32297699],
       [-0.19629232,  1.71053511,  0.13675151],
       [ 1.04189789, -0.46177364,  0.52251651],
       [-0.67053805,  0.96864843, -0.41905599]])
"""
np.dot(arr.T,arr)
"""
array([[ 2.82236803,  1.11445376, -0.67992213],
       [ 1.11445376,  4.67568549, -2.74324947],
       [-0.67992213, -2.74324947,  9.72932224]])
"""
  • 针对更高维的数组，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转置特殊案例：swapaxes方法接收一对轴编号作为参数，对轴进行调整用于重组数据(该方法返回视图，不复制)

a.swapaxes(x,y)，是将n维数组中两个维度进行调换，其中x，y的值为a.shape值（2，3，4）元组中的索引值（下标）

arr
"""
array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7]],

       [[ 8,  9, 10, 11],
        [12, 13, 14, 15]]])
"""
arr.shape # (2,2,4)
arr.swapaxes(1,2).shape # (2, 4, 2)
arr.swapaxes(1,2)
"""
array([[[ 0,  4],
        [ 1,  5],
        [ 2,  6],
        [ 3,  7]],

       [[ 8, 12],
        [ 9, 13],
        [10, 14],
        [11, 15]]])
"""

4.2 通用函数：快速的逐元素数组函数

通用函数，ufunc，在ndarray数据中进行逐元素操作的函数，是对简单函数的向量化封装。

  • 一元通用函数：sqrt和exp函数
  • 二元通用函数：add和maximum，接收两个数组并返回一个数组作为结果
  • 返回多个数组：比如modf，Python内建函数divmod的向量化版本，返回一个浮点值数组的小数部分和整数部分
import numpy as np
arr = np.random.randn(7) * 5
arr
"""
array([-1.14629029,  1.38707903, -1.30572811, -4.18866457, -4.95329571,
       -3.07819829,  2.87052789])
"""
remainder,whole_part = np.modf(arr)
remainder
"""
array([ 0.28959604,  0.54493711, -0.52703059,  0.92737401,  0.27334593,
        0.22031531, -0.75180789])
"""
whole_part
"""
array([ 6., -7., -2., -2., -0., -5.,  2.])
"""

4.3 使用数组进行面向对象数组编程

  • np.meshgrid函数接收两个一维数组，并根据两个数组的所有(x, y)对生成一个二维矩阵

将条件逻辑作为数组操作

  • numpy.where函数是三元表达式x if condition else y的向量化版本，不仅可以避免很多问题，速度更快。
  • where的典型用法是根据一个数组来生成一个新的数组
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])
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]
# 使用np.where方法
result = np.where(cond,xarr,yarr)
result # array([1.1, 2.2, 1.3, 1.4, 2.5])
  • 传递给np.where的数组既可以是同等大小的数组，也可以是标量
# 将arr中所有正值替换为常数2
arr = np.random.randn(4,4)
np.where(arr>2,2,arr)
"""
array([[-1.54828495, -0.11515345,  0.64979016, 1.1025353],
       [ 1.26570479, -0.38811792, -0.05638353, -0.23453037],
       [ 0.26468309, -0.33101883, -0.30085075, -1.63699033],
       [-0.5281036 ,  0.9839786 ,  0.64820938,  0.92989919]])
"""

数学和统计学方法

  • 聚合函数：sum、mean和std
  • 可选参数axis：计算给定轴向上的统计值，形成下降一维度的数组
  • arr.mean(1)表示“计算每一列的平均值”
  • arr.sum(0)表示“计算行轴向的累和”
  • cumsum和cumprod不会聚合，产生一个中间结果
arr = np.array([[0,1,2],[3,4,5],[6,7,8]])
arr
"""
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)
"""
  • 基础数组统计方法

布尔值数组的方法

  • any和all：any检查数组中是否至少有一个True，all检查是否每个值都是True

也适用于非布尔值数组，所有的非0元素都是True

排序

  • 和Python内建列表相似，使用sort方法按位置排序
  • 在多维数组中传递axis值，沿着轴向对一个维度数据进行排序
  • 顶层np.sort方法返回排好序的数组拷贝

唯一值与其他集合逻辑

  • NumPy包含一些针对一维ndarray的基础集合操作
  • 常用方法np.unique：返回数组中唯一值排序后形成的数组
names = np.array(['Bob','Joe','Will','Bob','Will','Joe','Joe'])
np.unique(names)
# array(['Bob', 'Joe', 'Will'], dtype='<U4')
# 纯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.4 使用数组进行文件输入和输出

  • 读取表格数据更多使用pandas
  • np.save和np.load是两大工具函数，默认以未压缩格式进行存储，后缀名是.npy
  • 操作代码
arr = np.arange(10)
np.save('some_array',arr) # .npy后缀名会被自动加上
# 载入数据
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)
# 载入一个．npy文件，获得字典型对象，可以方便地载入单个数组
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)

4.5 线性代数

  • 和Matlab等相比，NumPy的*是矩阵的逐元素乘积，而不是矩阵的点乘积
  • NumPy的数组方法和numpy命名空间中都有一个函数dot，用于矩阵的操作
x = np.array([[1,2,3],[4,5,6]])
y = np.array([[6,23],[-1,7],[8,9]])
# 下面两个函数结果相同
x.dot(y)
np.dot(x,y)
# 特殊符号@也作为中缀操作符，用于点乘矩阵操作
np.ones(3) # array([1., 1., 1.])
x @ np.ones(3) # array([ 6., 15.])

numpy.linalg

  • umpy.linalg拥有一个矩阵分解的标准函数集，以及其他常用函数，例如求逆和行列式求解。
from numpy.linalg import inv,qr
X = np.random.randn(5,5)
# X和它的转置矩阵X.T的点乘积
mat = X.T.dot(X)
inv(mat)
"""
array([[ 0.68887032, -0.10901956,  1.11151606, -0.34320821, -0.03279823],
       [-0.10901956,  0.2706275 ,  0.14261334,  0.00853147,  0.11737419],
       [ 1.11151606,  0.14261334,  4.19833265, -1.61815629,  0.41947259],
       [-0.34320821,  0.00853147, -1.61815629,  1.06411834, -0.21598402],
       [-0.03279823,  0.11737419,  0.41947259, -0.21598402,  0.20790649]])
"""
mat.dot(inv(mat))
q,r = qr(mat)
r
"""
array([[-5.07518215, -1.18301817,  1.8596161 ,  3.73253357,  5.72862095],
       [ 0.        , -1.4987164 ,  2.27029981,  0.47419322,  2.88815351],
       [ 0.        ,  0.        , -3.40688982, -1.45302848, -4.59955239],
       [ 0.        ,  0.        ,  0.        , -0.56543743, -0.04838235],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.97738804]])
"""
  • 常用的numpy.linalg函数

4.6 伪随机数生成

  • Python内建的random模块一次只能生成一个值，numpy.random模块填补了这项不足，可以高效地生成多种概率分布下的完整样本值数组。
  • 生成的是伪随机数：因为它们是由具有确定性行为的算法根据随机数生成器中的随机数种子生成的
  • np.random.seed更改NumPy的随机数种子：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.7 示例：随机漫步

  • 使用np.random模块随机漫步
nsteps = 1000
draws = np.random.randint(0,2,size=nsteps)
steps = np.where(draws > 0,1,-1)
walk = steps.cumsum()
# 统计漫步数据
walk.min() # -28
# walk.max() # 23
  • 找到第一次在某个方向连续走了10步的位置
# np.abs(walk) >= 10返回一个布尔值数组，表示是否连续在同一个方向走了10步
# argmax返回布尔值数组中最大值的第一个位置（True就是最大值）
# argmax效率不高，因为它总是完整的扫描整个数组
(np.abs(walk) >= 10).argmax() # 387

一次性模拟多次随机漫步

  • 二维的抽取数组，一次性地跨行算出全部5,000个随机步的累积和
nwalks = 5000
nsteps = 1000
draws = np.random.randint(0,2,size=(nwalks,nsteps))
steps = np.where(draws>0,1,-1)
walks = steps.cumsum(1)
walks
  • 随机步中计算出30或-30的最小穿越时间
  • any方法检查
hits30 = (np.abs(walks)>=30).any(1)
hits30 # array([False,  True,  True, ...,  True, False,  True])
hits30.sum() # 达到30或-30的数字 3367
  • 使用布尔值数组选出绝对步数超过30步所在的行，使用argmx从轴向1上获取穿越时间
crossing_times = (np.abs(walks[hits30])>=30).argmax(1)
crossing_times.mean()
# 508.3691713691714