3  NumPy基础:数组与向量化计算

3.1 引言与学习目标

学习目标

完成本章后,你应能:

  • 创建具有指定 shapedtypendarray,并根据 strides 解释二维索引对应的内存位置;
  • 判断切片、花式索引与布尔索引返回视图还是副本,并用一次原地修改验证判断;
  • 对给定数组写出广播后的目标形状,识别不兼容维度并修正计算;
  • 使用聚合、布尔过滤、矩阵乘法和线性代数函数计算题设收益与组合方差,并核对形状和标量结果;
  • 从本地上证综指行情计算简单收益、对数收益及其累计和,准确区分累计对数收益与超额收益。

目标—活动—核心习题/答案映射

正式目标 正文活动 核心评价证据
解释 shapedtypestrides 小节 3.3.2.1小节 3.3.4 的内存地址推导 习题 3.1(f) 及答案中的元数据和字节偏移核验
判断视图与副本并原地验证 小节 3.5.1小节 3.7小节 3.6 的对比 习题 3.2(f) 及答案中的 shares_memory、原地修改和断言
预测广播形状并修复不兼容维度 小节 3.4 的形状规则与反例 习题 3.5(a)、(d)—(e) 及答案中的形状断言和修正
使用数组统计与线性代数 小节 3.10.2小节 3.10.3 和转置矩阵活动 习题 3.3—3.5 及答案中的过滤、组合方差与标量核验
区分本地指数的三类收益对象 小节 3.12.2 的真实上证综指路径 习题 3.6 及答案中对教师提供 ndarray 使用 diffcumprodcumsum 与恒等式断言

前置知识与考核边界

读者应会使用第 2 章的列表、切片、循环和函数,并了解均值、方差、矩阵乘法与单期收益率。核心习题评价数组元数据、索引、广播、线性代数和本地指数收益口径。真实指数的 HDF5 读取、Pandas 筛选、日期索引与对象转换由教师脚手架提供,学生只运行并核对输出的日期数组和价格数组;这些 Pandas 操作不要求解释、修改或默写,也不评分。几何布朗运动、风险中性定价、Black–Scholes 与蒙特卡洛期权定价依赖连续时间随机过程和衍生品定价知识,只在章末给出选修路线图,不纳入本章核心评价。

NumPy,全称 Numerical Python,是 Python 数值计算生态中的基础包之一。对于经济学和金融学的学生,本章的直接任务是学会把同质数值组织成数组,并用索引、逐元素运算和聚合回答可核对的问题。pandasstatsmodelsscikit-learn 等库也常与 NumPy 数组交换数据,因此这些基本操作会在后续章节反复出现。

你将在 NumPy 中使用以下几类功能:

  • ndarray: 一种高效的多维数组,提供快速的面向数组的算术运算和灵活的广播(broadcasting)功能。
  • 数学函数: 用于对整个数据数组进行快速运算,而无需编写显式的 Python 循环。
  • 数据 I/O: 用于在磁盘上读/写数组数据以及处理内存映射文件的工具。
  • 线性代数和随机数生成: 提供全面的线性代数运算、随机数生成和傅里叶变换功能。
  • C API: 用于把 NumPy 与 C、C++ 或 FORTRAN 库连接起来的接口;本章只需知道许多数值循环由编译实现执行,不要求使用该接口。

NumPy 本身不负责完整的数据建模流程;本章先建立可观察的数组操作,再解释这些操作背后的形状、类型和内存规则。

NumPy 在经济学和金融学中的核心地位

作为经济学专业的学生,你可能会问:为什么要从数组学起?因为后续的表格、统计模型和机器学习接口最终都需要处理有形状、有类型的数值集合。理解内存布局和向量化有助于减少许多数值任务的 Python 解释器开销;实际收益仍取决于数组规模、dtype、临时数组、BLAS 实现与硬件,必须对具体任务测量后再下结论。

3.2 先做一个最小数组任务

先不讨论指针或步长。下面把四个题设单价放进一维数组,观察“创建—索引—切片—运算”四个动作。数值只是确定性的教学输入,不是市场观测。

import platform  # 识别操作系统,为本章后续真实指数数据读取建立统一入口
import numpy as np  # 导入NumPy并采用全书统一的np别名
DATA_ROOT = 'C:/qiufei/data' if platform.system() == 'Windows' else '/home/ubuntu/r2_data_mount/data'  # 根据操作系统选择规范数据根路径
unit_prices_cny = np.array([12.0, 15.5, 9.8, 20.0])  # 创建四项题设单价的一维浮点数组,单位为元
unit_prices_cny  # 展示数组保留输入顺序且使用统一数值类型
array([12. , 15.5,  9.8, 20. ])
second_price_cny = unit_prices_cny[1]  # 用位置1取得第二项单价,索引从0开始
middle_prices_cny = unit_prices_cny[1:3]  # 用半开切片取得位置1和2的两项单价
print(second_price_cny, middle_prices_cny)  # 同时核对标量索引与一维切片的输出形状
15.5 [15.5  9.8]
tax_rates = np.array([0.03, 0.03, 0.05, 0.05])  # 为四项单价提供逐位置对应的题设税率
tax_inclusive_prices_cny = unit_prices_cny * (1.0 + tax_rates)  # 对应位置逐元素相乘得到含税单价
average_price_cny = tax_inclusive_prices_cny.mean()  # 沿唯一轴归约为一项平均含税单价
print(tax_inclusive_prices_cny, average_price_cny)  # 展示向量结果和标量聚合结果的区别
[12.36  15.965 10.29  21.   ] 14.90375

这个最小例子已经给出本章主线:数组保存同质元素,索引选择位置,逐元素运算保持形状,聚合把一条轴压缩为统计量。接下来再解释为什么这些表达式成立,以及 dtypeshapestrides 如何限定其行为。

3.3 机制深化:向量化与 ndarray 元数据

3.3.1 理论基础:向量化计算的数学原理

向量化计算的理论基础

向量化计算的核心思想是将标量操作(scalar operations)转换为向量操作(vector operations),从而避免显式循环。

从数学角度,向量化可以形式化地定义为:

给定: - 向量 \(\mathbf{x} \in \mathbb{R}^n\) - 标量 \(c \in \mathbb{R}\) - 一元函数 \(f: \mathbb{R} \to \mathbb{R}\)

逐元素操作定义为: \[ f(\mathbf{x})_i = f(x_i), \quad \forall i \in \{1, 2, ..., n\} \]

标量广播定义为: \[ (c \odot \mathbf{x})_i = c \cdot x_i, \quad \forall i \in \{1, 2, ..., n\} \]

其中 \(\odot\) 表示广播乘法运算。

计算复杂度对比

操作类型 Python 层实现 NumPy 实现 应如何判断性能
元素级加法 O(n),逐元素解释执行 O(n),编译内核,可能产生临时数组 在相同 dtype、形状与硬件上实测
元素级乘法 O(n),逐元素解释执行 O(n),编译内核 同时核对内存带宽与临时数组
矩阵乘法 朴素实现通常为 O(n³) 通常调用优化 BLAS,复杂度仍由算法决定 记录矩阵大小、线程和 BLAS 实现
两个向量点积 O(n),逐元素解释执行 O(n),可能使用 SIMD 小数组可能受调用开销主导

其中 SIMD(Single Instruction, Multiple Data)是现代 CPU 的并行计算特性。NumPy 的一些编译内核可能使用 SIMD;是否使用以及收益大小取决于算子、dtype、NumPy 构建方式和硬件。向量化更一般的含义,是把逐元素循环从 Python 解释器移交给编译内核,而不是承诺每个算子都使用 SIMD。

NumPy 之所以对数值计算如此重要,其主要原因之一是它为处理大型数据数组而设计的高效率。这种效率源于几个关键的架构决策:

  • 内存布局: 许多数值 ndarray 基数组使用同质、定宽的紧凑缓冲区,但切片、转置等操作产生的视图可以通过步长访问非连续内存。编译内核利用已经确定的 dtype,减少 Python 逐元素动态分派;具体操作仍可能检查类型和形状,也可能创建临时数组。在表示相同的同质数值数据时,NumPy 数组通常比保存 Python 对象引用的序列更紧凑,但这不是对对象 dtype 或任意布局的无条件结论。
  • 向量化操作: NumPy 允许把许多同质数组循环委托给预编译内核。它常比等价的 Python 逐元素循环更快,但不承诺固定倍数,也不意味着含分支、对象 dtype 或大量临时数组的任务都应强行向量化。

为了说明性能差异,我们来比较一个包含一百万个整数的 NumPy 数组和一个等效的 Python 列表。

import numpy as np                                          # 用NumPy构造百万级同质数组并统计重复计时的中位数
import timeit  # 对两种等价计算执行可重复的进程内计时

price_array = np.arange(1_000_000)                          # 构造0至999999的一百万项整数数组作为向量化乘法基准输入
price_list = list(range(1_000_000))                         # 一百万个Python整数作为逐项循环基准,与数组保持相同输入规模

numpy_timings = timeit.repeat('price_array * 2', globals=globals(), number=10, repeat=5)  # 重复测量百万级NumPy数组的向量化乘法
list_timings = timeit.repeat('[value * 2 for value in price_list]', globals=globals(), number=10, repeat=5)  # 用等价列表推导式建立对照计时
print(f'NumPy 中位单次耗时: {np.median(numpy_timings) / 10:.6f} 秒')  # 报告多次测量的稳健中心而非单次结果
print(f'Python 列表中位单次耗时: {np.median(list_timings) / 10:.6f} 秒')  # 使用相同重复次数报告列表基准
NumPy 中位单次耗时: 0.000302 秒
Python 列表中位单次耗时: 0.041419 秒

%timeit 给出的只是当前解释器、硬件和这一百万整数乘法的测量结果。应先核对列表与数组计算是否等价,再报告中位耗时与离散程度;不能把一次观测外推为所有算子的固定加速倍数。此后我们会在适合数组内核的任务中使用 NumPy,并在索引对齐、I/O 或对象类型占主导时重新测量。

3.3.2 NumPy ndarray:一个多维数组对象

NumPy 的核心对象是 N 维数组 ndarray。它保存同质元素,并允许你用接近标量运算的语法对整块数组执行逐元素计算。

3.3.2.1 ndarray 对象内部机制

内存布局的数学表示

从数学角度,NumPy数组可以定义为五元组: \[ \text{ndarray} = (B, S, T, \mathbf{s}, P) \]

其中: - \(B\) 是底层字节缓冲区,按 \(T\) 解释后可表示整数、浮点数、复数、布尔、日期时间、定长字符串、结构化记录或对象引用等 NumPy dtype,并不限于实数域 - \(S \in \mathbb{N}_0^k\) 是形状(\(k\) 维度,维长允许为 0) - \(T\) 是数据类型(dtype),规定每个元素的字节数和解释规则 - \(\mathbf{s} \in \mathbb{Z}^k\) 是以字节为单位的步长(strides)向量;反向切片会产生负步长 - \(P\) 是当前数组第一个逻辑元素所对应的内存地址

线性索引公式

给定一个数组 \(A\),其形状为 \(S = (n_0, n_1, ..., n_{k-1})\),步长为 \(\mathbf{s} = (s_0, s_1, ..., s_{k-1})\),索引 \(\mathbf{i} = (i_0, i_1, ..., i_{k-1})\) 对应的线性地址为:

\[ \text{address}(\mathbf{i}) = P + \sum_{j=0}^{k-1} i_j \cdot s_j \]

其中 \(0 \leq i_j < n_j\)

行优先序(Row-Major Order)

NumPy默认使用C风格的行优先序,即最后一个维度变化最快: \[ s_{k-1} = \text{sizeof}(T) \] \[ s_{j} = s_{j+1} \cdot n_{j+1}, \quad j = 0, 1, ..., k-2 \]

例如,一个 \((3, 4, 5)\) 形状的 float64 数组,其步长为: \[ \mathbf{s} = (160, 40, 8) \]

视图的数学定义

给定一个 ndarray \(A = (B, S, T, \mathbf{s}, P)\),一个视图(view)通过创建新的元数据元组 \((B, S', T', \mathbf{s}', P')\) 解释同一个底层缓冲区 \(B\),而不是复制缓冲区。

切片操作 \(A[\text{slice}]\) 返回一个视图,满足: \[ P' = P + \sum_{j=0}^{k-1} a_j s_j \] \[ B' = B \]

其中 \(a_j\) 是切片在第 \(j\) 轴选中的首个逻辑索引。例如一维连续数组的步长若为 (8,),反向视图 arr[::-1] 的步长通常为 (-8,),且 \(P'\) 指向原数组最后一个元素。因而“步长为正整数”只适用于某些正向连续布局,不是 ndarray 的一般定义。

为了判断视图、复制和轴移动的行为,需要进一步检查 ndarray 的元数据。NumPy ndarray 提供了一种将一块同质类型的数据(连续或跨步的内存块)解释为多维数组对象的方式。

ndarray 之所以如此灵活,部分原因在于每个数组对象都是数据块上的一个跨步视图(strided view)。例如,像 arr[::2, ::-1] 这样的数组视图不会复制任何数据。原因是 ndarray 不仅仅是一块内存;它还包含了跨步信息,使数组能够以不同的步长在内存中移动。更精确地说,ndarray 内部由以下部分组成:

  • 一个指向数据的指针——即指向内存(RAM)或内存映射文件中的数据块。
  • 数据类型dtype,描述数组中固定大小的值单元。
  • 一个表示数组形状的元组。
  • 一个跨度(strides)元组——整数元组,指示为了沿某个维度前进一个元素需要“跨过”的字节数。

图 3.1 展示本节讨论对象的可视化结果,读图时应结合正文给出的口径与限制。

import matplotlib.pyplot as plt                             # 用Matplotlib把ndarray元数据与底层缓冲区关系画成结构图
import matplotlib.patches as patches  # 提供圆角框、矩形和箭头图元以编码ndarray的四类元数据关系

fig, ax = plt.subplots(figsize=(10, 6))                     # 建立显式坐标轴,以分别控制数据映射和出版格式
ax.set_xlim(0, 10)                                          # 横向数据坐标覆盖结构框与三类元数据框的完整宽度
ax.set_ylim(0, 7)                                           # 纵向数据坐标为主对象、指针和元数据框保留布局空间
ax.axis('off')  # 关闭坐标轴刻度和边框,仅保留图形示意

# ndarray主对象:统摄指针与三类数组元数据
ax.add_patch(patches.FancyBboxPatch((1, 4.5), 8, 2, boxstyle='round,pad=0.1', fc='lightblue', ec='black'))  # 用顶部圆角框表示统摄指针、dtype、shape和strides的ndarray对象
ax.text(5, 5.5, 'NumPy ndarray Object', ha='center', va='center', fontsize=16, fontweight='bold')  # 标注ndarray对象标题

# 数据指针:连接逻辑数组与共享内存缓冲区
ax.add_patch(patches.Rectangle((1, 0.5), 2.5, 2.5, fc='lightgreen', ec='black'))  # 左下绿色框单独编码指向共享底层缓冲区的数据指针
ax.text(2.25, 1.75, 'Data Pointer\n(to memory block)', ha='center', va='center', fontsize=12)  # 标注数据指针说明
ax.add_patch(patches.Arrow(5, 4.5, 0, -1.2, width=0.3, color='gray'))  # 中央向下连线表明数据指针是ndarray对象持有的元数据之一

# dtype:规定缓冲区每个8字节单元的解释方式
ax.add_patch(patches.Rectangle((4, 2.5), 2, 1, fc='lightcoral', ec='black'))  # 中部红框表示每个缓冲区单元按float64解释的dtype合同
ax.text(5, 3, 'dtype: float64', ha='center', va='center', fontsize=12)  # float64说明每个缓冲区单元按双精度浮点数解释
ax.add_patch(patches.Arrow(3.5, 4.5, 1, -0.8, width=0.2, color='gray'))  # 左斜连线把dtype解释规则关联回ndarray主对象


# shape:规定两条逻辑轴分别含10项和5项
ax.add_patch(patches.Rectangle((7, 2.5), 2, 1, fc='lightgoldenrodyellow', ec='black'))  # 右中黄框记录示例数组两条轴长度为10和5
ax.text(8, 3, 'shape: (10, 5)', ha='center', va='center', fontsize=12)  # (10,5)把可用索引边界限定为轴0的0—9与轴1的0—4
ax.add_patch(patches.Arrow(6.5, 4.5, 1, -0.8, width=0.2, color='gray'))  # 右斜连线表明shape决定逻辑索引域而不复制缓冲区

# strides:规定沿两轴前进一步的字节偏移
ax.add_patch(patches.Rectangle((7, 0.5), 2, 1, fc='plum', ec='black'))  # 右下紫框给出沿两轴前进一步分别跨40和8字节
ax.text(8, 1, 'strides: (40, 8)', ha='center', va='center', fontsize=12)  # 40/8字节分别对应跨一整行与跨一个float64单元
ax.add_patch(patches.Arrow(5, 4.5, 2.5, -2.8, width=0.2, color='gray'))  # 最长连线强调strides把逻辑轴移动映射到底层字节偏移

ax.set_title('ndarray 内部结构示意图', fontsize=18, pad=20)        # 标明本图展示 shape、dtype、strides 与数据指针的关系
plt.show()                                                  # 输出shape、dtype、strides与数据指针之间的层级关系图
图 3.1: NumPy ndarray 对象的内部结构

例如,一个 10 × 5 的数组其形状为 (10, 5)。

np.ones((10, 5)).shape                                      # 返回示例数组两条轴的长度(10,5),不读取任何元素值
(10, 5)

一个典型的 C 顺序(行主序)的 3 × 4 × 5 float64(8 字节)数组,其跨度为 (160, 40, 8)。跨度说明相邻逻辑元素在缓冲区中的字节距离;它会影响缓存访问模式,但不能单独决定总耗时。

np.ones((3, 4, 5), dtype=np.float64).strides  # 查看3×4×5 float64数组的内存步长(单位:字节)
(160, 40, 8)

虽然大多数用户不会直接与跨度(strides)交互,但它们是“零拷贝”视图机制的基础。跨度甚至可以是负数,这使得数组能够“向后”移动通过内存,例如在 obj[::-1] 这样的切片中。

3.3.3 创建函数与形状控制

在最小任务已经展示 np.array 之后,本节补充嵌套序列、预分配函数和显式形状控制。

import numpy as np                                          # 用NumPy把嵌套收益与波动率序列转为同质ndarray

# 从一个列表的列表创建一个二维数组(模拟波动率数据)
volatility_data = np.array([[1.5, -0.1, 3], [0, -3, 6.5]])  # 用嵌套列表创建2×3的波动率数据矩阵
volatility_data  # 展示两条观测、三个指标的原始浮点矩阵以对照后续逐元素运算
array([[ 1.5, -0.1,  3. ],
       [ 0. , -3. ,  6.5]])

现在,我们可以直接在这个 data 对象上执行数学运算:

# 逐元素乘法
volatility_data * 10  # 将波动率数据的每个元素扩大10倍(向量化逐元素乘法)
array([[ 15.,  -1.,  30.],
       [  0., -30.,  65.]])
# 逐元素加法
volatility_data + volatility_data  # 波动率数据与自身逐元素相加(等价于乘以2)
array([[ 3. , -0.2,  6. ],
       [ 0. , -6. , 13. ]])

在第一个例子中,所有元素都乘以 10。在第二个例子中,数组中每个“单元格”的对应值被相加。

ndarray 是一个通用的多维同质数据容器;也就是说,所有元素必须是相同类型。每个数组都有一个 shape(一个表示各维度大小的元组)和一个 dtype(一个描述数组数据类型的对象):

# 读取volatility_data两条轴长度,确认嵌套输入没有被展平
volatility_data.shape                                       # 核对外层两组观测与每组三个指标形成(2,3)形状
(2, 3)
# 检查混合整数与小数输入经类型提升后的统一dtype
volatility_data.dtype  # 核对整数与小数混合输入被统一提升为浮点dtype
dtype('float64')

创建数组最简单的方法是使用 array 函数。它接受任何序列型对象(包括其他数组),并生成一个包含传入数据的新 NumPy 数组。

# 从列表创建(模拟收益率)
returns_data = [6, 7.5, 8, 0, 1]                            # 五期同质数值含整数与小数,用于观察数组统一提升dtype
returns_arr = np.array(returns_data)                        # 将Python列表转换为NumPy一维数组
returns_arr  # 展示列表转换后保持五项顺序的一维收益数组
array([6. , 7.5, 8. , 0. , 1. ])

嵌套序列,比如一个等长列表组成的列表,将被转换成一个多维数组:

portfolio_matrix = [[1, 2, 3, 4], [5, 6, 7, 8]]             # 定义包含两个等长列表的嵌套列表
portfolio_arr = np.array(portfolio_matrix)                  # 将嵌套列表转换为二维NumPy数组
portfolio_arr  # 展示两条组合记录各四项数值形成的二维布局
array([[1, 2, 3, 4],
       [5, 6, 7, 8]])

由于 portfolio_matrix 是一个列表的列表,NumPy 数组 portfolio_arr 有两个维度。我们可以通过检查 ndimshape 属性来确认这一点:

portfolio_arr.ndim  # 返回轴数量2,与嵌套列表的外层和内层对应
2
portfolio_arr.shape                                         # 核对两项组合记录、每项四个数值对应(2,4)形状
(2, 4)

除了 np.array,还有许多其他函数可以创建新数组。np.zerosnp.ones 分别创建全为 0 或全为 1 的数组。np.empty 创建一个数组,但不将其值初始化为任何特定值。要用这些方法创建更高维度的数组,只需传入一个元组作为形状参数。

# 标量形状10只建立一条轴,结果可作为十期观测的预分配容器
np.zeros(10)  # 返回shape为(10,)的一维浮点数组,不等同于10×1列向量
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
# 二元形状元组分别规定轴0含3行、轴1含6列
np.zeros((3, 6))  # 返回shape为(3,6)的二维浮点数组,供比较一维与二维形状
array([[0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0.]])
# 创建一个未初始化的2x3x2数组
np.empty((2, 3, 2))  # 创建2×3×2的三维未初始化数组(值为随机内存残留)
array([[[0., 0.],
        [1., 0.],
        [1., 1.]],

       [[0., 1.],
        [0., 0.],
        [0., 0.]]])

关于 np.empty 的重要警告

假设 numpy.empty 会返回一个全零数组是不安全的。此函数返回的是未初始化的内存,因此可能包含任意的“垃圾”值。在金融建模中,例如初始化投资组合权重矩阵时,你应当只在确定会立即用真实数据覆盖该数组的情况下使用此方法,以避免计算中引入随机不可控的数值误差。

np.arange 是内置 Python range 函数的数组版本:

np.arange(15)  # 生成0到14的等差整数数组(类似Python的range)
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])

表 3.1 列出了一些标准的数组创建函数。

表 3.1: 一些重要的 NumPy 数组创建函数。
函数 描述
array 将输入数据(列表、元组、数组或其他序列类型)转换为 ndarray,可推断数据类型或显式指定;默认复制输入数据。
asarray 将输入转换为 ndarray,但如果输入已经是 ndarray 则不复制。
arange 类似于内置的 range,但返回一个 ndarray 而不是列表。
ones, ones_like 生成一个给定形状和数据类型的全 1 数组;ones_like 接受另一个数组并生成一个相同形状和数据类型的全 1 数组。
zeros, zeros_like 类似于 onesones_like,但生成的是全 0 数组。
empty, empty_like 通过分配新内存创建新数组,但不像 ones 和 zeros 那样用任何值填充。
full, full_like 生成一个给定形状和数据类型的数组,所有值都设置为指定的“填充值”;full_like 接受另一个数组并生成一个相同形状和数据类型的填充数组。
eye, identity 创建一个 N × N 的单位矩阵(对角线上为 1,其他地方为 0)。

3.3.4 ndarray 的数据类型

数据类型dtype 是一个特殊的对象,它包含了 ndarray 需要将一块内存解释为特定类型数据的信息(或元数据)。除非显式指定,numpy.array 会尝试推断一个合适的数据类型。

prices_float = np.array([6, 7.5, 8, 0, 1])                  # 含浮点数的数组,NumPy将自动推断为float64类型
prices_int = np.array([[1, 2, 3], [4, 5, 6]])               # 纯整数的二维数组,NumPy将自动推断为int64类型
print(f'prices_float 的 dtype: {prices_float.dtype}')        # 报告混合整数/小数输入推断出的统一浮点dtype
print(f'prices_int 的 dtype: {prices_int.dtype}')            # 报告纯整数二维输入推断出的平台整数dtype
prices_float 的 dtype: float64
prices_int 的 dtype: int64

因为 arr1 包含浮点数,其 dtype 被推断为 float64arr2 只包含整数,被推断为 int64(在64位系统上)。

你可以在创建时显式指定 dtype

prices_float = np.array([1, 2, 3], dtype=np.float64)        # 显式指定dtype为64位浮点数
prices_int = np.array([1, 2, 3], dtype=np.int32)            # 显式指定dtype为32位整数
print(f'prices_float 的 dtype: {prices_float.dtype}')        # 核对显式声明已把三个整数存为float64
print(f'prices_int 的 dtype: {prices_int.dtype}')            # 核对显式声明把数组限制为32位整数
prices_float 的 dtype: float64
prices_int 的 dtype: int32

数据类型是 NumPy 与来自其他系统的数据交互灵活性的源泉。它们提供了到磁盘或内存底层表示的直接映射。数值型 dtype 的命名方式是类型名(intfloatcomplex)后跟一个表示每个元素位数的数字。一个标准的双精度浮点值(Python 的 float 对象所使用的)占用 8 字节或 64 位,因此得名 float64表 3.2 列出了主要的数据类型。

表 3.2: 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 标准单精度浮点数
float64 f8 标准双精度浮点数 (与 Python 的 float 兼容)
complex64, complex128 c8, c16 由两个 32 位或 64 位浮点数表示的复数
bool ? 存储 TrueFalse 值的布尔类型
object O Python 对象类型
bytes_ S 固定长度 ASCII 字节串类型 (例如 S10 代表 10 字节字符串)
unicode_ U 固定长度 Unicode 类型 (例如 U10 代表 10 字符字符串)

有符号与无符号整数的财务数据处理细节

在经济和金融数据处理中,负数经常出现(例如季度亏损、负利率或收益率回撤)。因此,必须清楚 signed(有符号)和 unsigned(无符号)整数类型的区别。

  • signed 整数: 可以表示正数和负数(如 int8 范围为 -128 到 127)。
  • unsigned 整数: 只能表示非负数(如 uint8 范围为 0 到 255)。

在处理会计分录或盈亏数据时,错误地使用 unsigned 类型可能导致负数被错误地解释为一个非常大的正数,从而引发严重的财务计算逻辑泄露。

你可以使用 astype 方法显式地将一个数组从一个 dtype 转换或强制转换(cast)为另一个。

int_prices = np.array([1, 2, 3, 4, 5])                      # 创建整数数组(默认推断为int64)
print(f'原始 dtype: {int_prices.dtype}')                      # 先报告整数解释规则,作为显式转型后的类型对照

float_prices = int_prices.astype(np.float64)                # 将整数数组转换为浮点数数组
print(f'新的 dtype: {float_prices.dtype}')                    # 转型后应为float64,说明元素值与底层表示一并改变
print(f'新的数组: {float_prices}')  # 查看转换后的数组内容(整数变为带小数点的浮点数)
原始 dtype: int64
新的 dtype: float64
新的数组: [1. 2. 3. 4. 5.]

如果你将浮点数转换为整数 dtype,小数部分将被截断:

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.bytes_)  # 创建字节串数组(模拟从文本文件读入的数字)
numeric_strings.astype(float)                               # 将字节串解析为浮点数(自动处理编码转换)
array([ 1.25, -9.6 , 42.  ])

astype 的复制条件

astype 默认使用 copy=True,因此通常返回新数组。若显式指定 copy=False,并且目标 dtype、内存顺序和子类要求已经由输入满足,NumPy 可以直接返回原数组;只要发生 dtype 转换、内存顺序调整或其他不兼容要求,仍会分配新缓冲区。因此 copy=False 是允许避免复制的请求,不是零拷贝保证。处理数千万行财务数据时,应在转换前估算峰值内存,并用对象身份或 np.shares_memory 检查实际结果。

same_dtype = float_prices.astype(float_prices.dtype, copy=False)  # dtype与布局已满足时请求复用原数组
converted_dtype = float_prices.astype(np.float32, copy=False)  # dtype改变时即使copy=False也需要新的数据表示
print(f'同dtype返回原对象: {same_dtype is float_prices}')  # 核对当前NumPy版本是否直接返回输入对象
print(f'转换dtype共享内存: {np.shares_memory(converted_dtype, float_prices)}')  # 验证实际类型转换没有共享元素缓冲区
同dtype返回原对象: True
转换dtype共享内存: False

3.3.5 NumPy 数据类型层次结构

有时需要检查一个数组的 dtype 是否属于某个通用类别,比如整数或浮点数。数据类型有诸如 np.integernp.floating 这样的超类,可以与 np.issubdtype 函数一起使用:

ints = np.ones(10, dtype=np.uint16)                         # 创建10个元素的无符号16位整数全一数组
floats = np.ones(10, dtype=np.float32)                      # 创建10个元素的单精度浮点数全一数组

print(f'ints.dtype 是 np.integer 的子类型吗? {np.issubdtype(ints.dtype, np.integer)}')  # 检查uint16是否属于整数大类
print(f'floats.dtype 是 np.floating 的子类型吗? {np.issubdtype(floats.dtype, np.floating)}')  # 检查float32是否属于浮点数大类
ints.dtype 是 np.integer 的子类型吗? True
floats.dtype 是 np.floating 的子类型吗? True

你可以通过调用特定数据类型的 mro 方法查看它的所有父类:

np.float64.mro()  # 查看float64类型的完整继承链(方法解析顺序)
[numpy.float64,
 numpy.floating,
 numpy.inexact,
 numpy.number,
 numpy.generic,
 float,
 object]
import matplotlib.pyplot as plt                             # 用Matplotlib递归绘制NumPy标量类型的父子继承层级

def plot_hierarchy(ax, hierarchy, x, y, dx, dy):            # 接收当前子树及坐标间距,把嵌套类型关系递归展开到同一坐标轴
    for key, value in hierarchy.items():                    # 当前字典的每个键是父节点,值决定继续展开子树还是落到叶节点
        ax.text(x, y, key, ha='center', va='center', bbox=dict(boxstyle='round,pad=0.5', fc='skyblue'))  # 在(x,y)处绘制带蓝色圆角背景框的节点标签
        if isinstance(value, dict):  # 判断子节点是否为字典(即还有更深层级)
            child_keys = list(value.keys())                 # 保留当前映射键顺序,使同层节点的视觉次序稳定
            num_children = len(child_keys)  # 子节点数决定横向间距,避免同层类型标签重叠
            child_x_start = x - (num_children - 1) * dx / 2  # 计算第一个子节点的起始x坐标(居中分布)
            for i, child_key in enumerate(child_keys):      # 遍历每个子节点及其索引
                child_x = child_x_start + i * dx  # 根据索引序号计算该子节点的x坐标
                child_y = y - dy  # 子节点y坐标向下偏移dy(下一层)
                ax.plot([x, child_x], [y - 0.3, child_y + 0.3], 'k-')  # 绘制从父节点到子节点的连接线
                plot_hierarchy(ax, {child_key: value[child_key]}, child_x, child_y, dx / num_children, dy)  # 递归绘制子节点及其下层结构
        elif isinstance(value, list):                       # 判断子节点是否为列表(即叶子节点集合)
            num_children = len(value)  # 计算叶子节点数量
            child_x_start = x - (num_children - 1) * dx / 2  # 计算叶子节点的起始x坐标
            for i, child_item in enumerate(value):          # 遍历每个叶子节点
                 child_x = child_x_start + i * dx  # 计算该叶子节点的x坐标
                 child_y = y - dy  # 叶子节点y坐标向下偏移
                 ax.plot([x, child_x], [y - 0.3, child_y + 0.3], 'k-')  # 绘制到叶子节点的连接线
                 ax.text(child_x, child_y, child_item, ha='center', va='center', bbox=dict(boxstyle='round,pad=0.3', fc='lightgreen'))  # 绘制带绿色背景的叶子节点标签

图 3.2 展示本节讨论对象的可视化结果,读图时应结合正文给出的口径与限制。

hierarchy = {                                               # 定义NumPy数据类型的层次结构字典
    'np.generic': {  # 根节点:所有NumPy类型的基类
        'np.number': {  # 数值类型分支
            'np.integer': {  # 整数类型分支
                'np.signedinteger': ['int8', 'int16', 'int32', 'int64'],  # 有符号整数(可表示负数)
                'np.unsignedinteger': ['uint8', 'uint16', 'uint32', 'uint64']  # 无符号整数(仅非负数)
            },  # 整数分支结束
            'np.floating': ['float16', 'float32', 'float64'],  # 浮点数类型(常用float64)
            'np.complexfloating': ['complex64', 'complex128']  # 复数类型
        },  # 数值类型分支结束
        'np.bool_': {},  # 布尔类型(无子节点)
        'np.character': {  # 字符类型分支
            'np.bytes_': [],  # 字节串类型
            'np.unicode_': []  # Unicode字符串类型
        }  # 字符类型分支结束
    }  # np.generic结束
}  # 层次结构字典定义完成

fig, ax = plt.subplots(figsize=(16, 10))                    # 创建16×10英寸的画布和坐标轴
ax.axis('off')  # 隐藏坐标轴(层次图不需要轴线)
plot_hierarchy(ax, hierarchy, 0.5, 1.0, 0.4, 0.15)  # 从中心点(0.5,1.0)开始递归绘制整个类型层次图
ax.set_title('NumPy Data Type Hierarchy', fontsize=20)      # 概括从 np.generic 到具体标量类型的继承层级
plt.show()                                                  # 输出从np.generic到具体数值与字符dtype的继承树
图 3.2: NumPy 数据类型类层次结构

3.4 NumPy 数组的运算

数组允许你表达对同质数据的批量操作,而不必在 Python 层逐元素编写 for 循环。这种写法称为向量化(vectorization)。等形状数组之间的算术运算按对应位置执行。

# 模拟资产收益率
asset_returns = np.array([[1., 2., 3.], [4., 5., 6.]])      # 创建2×3的浮点数数组(模拟两只资产三期收益率)
# 逐元素相乘
asset_returns * asset_returns  # 逐元素相乘:每个元素与自身相乘(即各元素的平方)
array([[ 1.,  4.,  9.],
       [16., 25., 36.]])
# 逐元素相减
asset_returns - asset_returns  # 逐元素相减:每个元素减去自身,结果为全零矩阵
array([[0., 0., 0.],
       [0., 0., 0.]])

与标量的算术运算会将标量参数传播到数组的每个元素。这是广播(broadcasting)的一种简单形式。

核心概念:广播 (Broadcasting)

广播的数学形式化定义

给定两个数组 \(A\)\(B\),其形状分别为 \(S_A \in \mathbb{N}_0^{n_A}\)\(S_B \in \mathbb{N}_0^{n_B}\),其中 \(\mathbb{N}_0\) 包含 0。令 \(n=\max(n_A,n_B)\),从左侧补 1 得到等长形状 \(\widetilde S_A,\widetilde S_B\in\mathbb{N}_0^{n}\)。广播能否成立以及输出形状由三条规则确定:

  1. 维度对齐:从尾部(右侧)开始对齐维度;缺失维度在左侧补 1。

  2. 兼容性检查:对每个维度 \(i\),必须满足 \(\widetilde S_{A,i}=\widetilde S_{B,i}\),或其中至少一个等于 1;只要有一个维度不满足,广播就不存在。

  3. 输出形状:仅在所有维度兼容时,若两个维长相等就取共同值;若不等,则其中一个必为 1,输出取另一个维长。逐维写为

    \[ S_{C,i}= \begin{cases} \widetilde S_{A,i}, & \widetilde S_{A,i}=\widetilde S_{B,i},\\ \widetilde S_{B,i}, & \widetilde S_{A,i}=1,\\ \widetilde S_{A,i}, & \widetilde S_{B,i}=1. \end{cases} \]

这一区分对零长度轴很重要:形状 (0, 3)(1, 3) 可以广播,第一维由“1 对应另一维”规则得到 0,因此输出形状是 (0, 3),而不是逐维取最大值得到的 (1, 3)

广播描述了 NumPy 在算术运算期间如何处理不同形状的数组。例如,在收益率归一化计算 returns - returns.mean(axis=0) 中,均值向量会被广播到原始矩阵的每一行。这是 NumPy 向量化运算兼顾效率和语义简洁的核心。

# 标量除法
1 / asset_returns  # 标量1除以数组的每个元素(广播机制:标量扩展到整个数组)
array([[1.        , 0.5       , 0.33333333],
       [0.25      , 0.2       , 0.16666667]])
# 标量指数运算
asset_returns ** 2  # 数组每个元素求平方(广播机制:标量2应用于每个元素)
array([[ 1.,  4.,  9.],
       [16., 25., 36.]])

同样大小数组之间的比较也会产生布尔数组:

benchmark_returns = np.array([[0., 4., 1.], [7., 2., 12.]])  # 创建基准收益率数组用于逐元素比较
benchmark_returns > asset_returns  # 逐元素比较:返回布尔数组,True表示基准收益高于资产收益
array([[False,  True, False],
       [ True, False,  True]])

3.5 基本索引和切片

在数据分析的过程中,我们经常需要从庞大的数据集中提取特定的子集——例如过去 5 个交易日的价格数据,或者涨幅超过 5% 的股票。NumPy 数组的索引(Indexing)切片(Slicing)是完成这类任务的基础工具。NumPy 数组索引包含多种子集选择方法;一维数组的基本操作与 Python 列表相似。

# 模拟日成交量
daily_volume = np.arange(10)                                # 生成0到9的等差整数数组(模拟10个交易日成交量索引)
daily_volume  # 展示0—9与位置一一对应的一维数组,便于追踪切片回写
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# 获取索引为 5 的元素
daily_volume[5]  # 获取索引为5的元素(即第6个元素,值为5)
5
# 获取从索引 5 到 8 (不包括) 的切片
daily_volume[5:8]  # 切片获取索引5到8(不含)的元素,返回[5,6,7]
array([5, 6, 7])

如果你将一个标量值赋给一个切片,比如 arr[5:8] = 12,这个值会被传播(或广播)到整个选区。

# 广播赋值
daily_volume[5:8] = 12  # 将索引5到8的元素全部赋值为12(广播赋值)
daily_volume  # 查看修改后的数组,索引5-7的元素变为12
array([ 0,  1,  2,  3,  4, 12, 12, 12,  8,  9])

核心规则:视图(Views) vs. 副本(Copies)

与 Python 内置列表不同,NumPy 中的数组切片(Slicing)返回的是底层数据的视图

  • 优势: 基本切片通常只创建新的数组元数据并共享原缓冲区,因此附加内存较小;耗时仍受索引方式、后续访问模式和缓存影响,且原始数组在视图存续期间不能释放。
  • 风险: 对视图的修改会直接“回溯”原数组数据。
  • 解决方案: 若需独立处理子集,必须显式调用 .copy()(例如 latest_prices = all_prices[:100].copy())。

为了演示这一点,我们首先创建 arr 的一个切片:

# 半开边界5:8命中位置5、6、7,并保留与原数组共享缓冲区的关系
volume_slice = daily_volume[5:8]  # 三项视图没有复制数据,后续原地赋值会回写对应原位置
volume_slice  # 展示位置5—7的三个值,并保留其与原数组共享缓冲区的语义
array([12, 12, 12])

现在,当我们改变 arr_slice 中的值时,这些改动会反映在原始数组 arr 中:

volume_slice[:] = 12345  # 修改切片中所有元素为12345
daily_volume  # 查看原数组:索引5-7也变为12345(视图特性)
array([    0,     1,     2,     3,     4, 12345, 12345, 12345,     8,
           9])

“裸”切片 [:] 会赋值给数组中的所有值:

volume_slice[:] = 64  # 再次修改切片中所有元素为64
daily_volume  # 原数组索引5-7同步变为64(进一步验证视图特性)
array([ 0,  1,  2,  3,  4, 64, 64, 64,  8,  9])

对于更高维度的数组,你有更多的选择。在一个二维数组中,每个索引处的元素是一维数组:

# 模拟价格矩阵 (3x3)
price_matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])  # 创建3×3的二维整数数组(模拟价格矩阵)

# 访问第三行 (索引 2)
price_matrix[2]  # 获取第三行(索引2),返回一维数组[7,8,9]
array([7, 8, 9])

单个元素可以通过递归方式访问,但更方便的是传递一个逗号分隔的索引列表:

# 这两种形式是等价的
print(f'递归访问: {price_matrix[0][2]}')  # 先取第0行再取第2个元素,结果为3
print(f'元组访问: {price_matrix[0, 2]}')  # 直接用逗号分隔行列索引,更简洁高效,结果同样为3
递归访问: 3
元组访问: 3

图 3.3 展示了在一个二维数组上的索引。将轴 0 想象为“行”,轴 1 想象为“列”,会很有帮助。

import matplotlib.pyplot as plt                             # 用Matplotlib按矩阵惯例绘制二维数组的行列索引坐标

def plot_indexing_diagram():                                # 把二维索引元组映射到3×3网格,直观区分轴0行位置与轴1列位置
    fig, ax = plt.subplots(figsize=(4, 4))                  # 为3×3索引网格建立等宽高坐标轴,保证单元格视觉尺寸一致
    ax.set_xlim(-0.5, 2.5)                                  # 横轴恰好包围列索引0—2对应的三个单元格
    ax.set_ylim(-0.5, 2.5)                                  # 纵轴恰好包围行索引0—2对应的三个单元格
    ax.invert_yaxis()  # 翻转y轴使行号从上到下递增(符合矩阵惯例)
    
    for i in range(3):  # 外层位置写入索引元组首项,对应从上到下的轴0行号
        for j in range(3):  # 内层位置写入索引元组次项,对应从左到右的轴1列号
            ax.add_patch(plt.Rectangle((j-0.5, i-0.5), 1, 1, fill=True, color='skyblue', alpha=0.6, ec='black'))  # 为每个单元格绘制浅蓝色矩形
            ax.text(j, i, f'({i}, {j})', ha='center', va='center', fontsize=14)  # 在单元格中心标注(行,列)索引

    ax.set_xticks(range(3))                                 # 横轴只标出二维数组的列位置0、1、2
    ax.set_yticks(range(3))                                 # 纵轴只标出二维数组的行位置0、1、2
    ax.xaxis.tick_top()  # 把列位置放到矩阵上缘,使读图顺序与表格列标题一致
    ax.xaxis.set_label_position('top')  # 同步上移“轴1(列)”标签,避免与下缘空白混淆
    ax.set_xlabel('轴 1 (列)', fontsize=12)                   # 横向刻度对应二维数组的列索引(axis=1)
    ax.set_ylabel('轴 0 (行)', fontsize=12)                   # 纵向刻度对应二维数组的行索引(axis=0)
    
    # 用红色覆盖(0,2),把“先行后列”的索引顺序落到具体单元格
    ax.add_patch(plt.Rectangle((2-0.5, 0-0.5), 1, 1, fill=True, color='tomato', alpha=0.8, ec='black'))  # 用红色高亮(0,2)位置的单元格
    ax.text(1.5, -1.0, "arr2d[0, 2]", color='tomato', fontsize=14, ha='center')  # 在图表上方标注该索引表达式

    ax.set_title('二维数组索引', loc='center', pad=40, fontsize=16)  # 提示读者按“先行后列”解释单元格坐标

    plt.show()                                              # 输出高亮(0,2)并标注axis=0/1含义的索引示意图

plot_indexing_diagram()  # 调用函数绘制二维数组索引示意图
图 3.3: 在一个二维 NumPy 数组中索引元素

在多维数组中,如果你省略了后面的索引,返回的对象将是一个较低维度的 ndarray 切片。对于一个 2 × 2 × 3 的数组 arr3d

# 模拟三维张量 (2x2x3)
tensor_data = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])  # 创建2×2×3的三维数组(模拟张量数据)
# 在轴 0 上切片返回一个 2x3 数组
tensor_data[0]  # 沿轴0取第0个切片,返回一个2×3的子数组
array([[1, 2, 3],
       [4, 5, 6]])

3.5.1 使用切片进行索引

和一维数组一样,多维数组也可以被切片。

# 选择 price_matrix 的前两行
price_matrix[:2]  # 选择price_matrix的前两行(沿轴0切片)
array([[1, 2, 3],
       [4, 5, 6]])

这是沿着轴 0 进行切片。你可以将 price_matrix[:2] 理解为“选择 price_matrix 的前两行”。你可以传递多个切片:

# 选择前两行以及从索引 1 开始的所有列
price_matrix[:2, 1:]  # 选择前两行,且每行只取从第1列开始到末尾的列
array([[2, 3],
       [5, 6]])

当你这样切片时,你总是得到相同维度的数组视图。通过混合使用整数索引和切片,你可以得到较低维度的切片。例如,要选择第二行但只取前两列:

lower_dim_slice = price_matrix[1, :2]  # 选择第1行的前两列(整数索引+切片,降维为一维)
lower_dim_slice  # 展示第1行前两列,整数行索引使结果由二维降为长度2
array([4, 5])

这里,lower_dim_slice 是一维的。 单独一个冒号 : 表示取整个轴:

# 选择所有行中的第一列
price_matrix[:, :1]  # 选择所有行的第0列(保持二维形状)
array([[1],
       [4],
       [7]])

图 3.4 提供了切片的可视化总结。

import numpy as np                                          # 用NumPy构造5×5编号矩阵及各切片区域的布尔掩码
import matplotlib.pyplot as plt                             # 用Matplotlib并排比较四种二维切片的命中单元格

def show_slice(ax, slice_str, highlight):                   # 将行列边界转成5×5布尔掩码,使切片表达式与保留形状可逐格核对
    arr = np.arange(25).reshape((5, 5))                     # 创建5×5的示例数组(0~24)
    ax.imshow(np.zeros_like(arr), cmap='Pastel2', vmin=-1, vmax=1)  # 用浅色底图绘制数组背景
    
    # 根据传入的行列边界生成布尔掩码,只给当前切片覆盖区着色
    rows, cols = highlight  # 解包高亮区域的行列切片参数
    if isinstance(rows, int): rows = slice(rows, rows + 1)  # 若行索引为整数则转为切片以统一处理
    if isinstance(cols, int): cols = slice(cols, cols + 1)  # 若列索引为整数则转为切片以统一处理
    
    highlight_mask = np.zeros_like(arr, dtype=bool)         # 初始化全False的布尔掩码矩阵
    highlight_mask[rows, cols] = True  # 将指定切片区域标记为True
    ax.imshow(np.where(highlight_mask, 1, 0), cmap='Set1', alpha=0.5, vmin=-1, vmax=1)  # 用红色半透明覆盖高亮区域

    for i in range(5):  # 轴0位置决定编号矩阵中每个文本标记的纵坐标
        for j in range(5):  # 轴1位置决定横坐标,并与掩码中的同一单元格配对
            ax.text(j, i, str(arr[i, j]), ha='center', va='center', color='black')  # 在每个单元格中心显示数值
    
    ax.set_title(f"arr{slice_str}", fontsize=14)            # 直接把当前切片表达式写入面板题名,便于与红色区域配对
    ax.set_xticks([])                                       # 移除列刻度,避免与单元格内0—24数值发生双重编码
    ax.set_yticks([])                                       # 移除行刻度,仅以高亮位置传达切片覆盖范围

fig, axes = plt.subplots(2, 2, figsize=(8, 8))              # 分配四个等大面板,一一承载两种二维切片和两种降维选择
fig.suptitle('二维数组切片示例', fontsize=20)  # 用总题名说明四面板共享同一5×5编号矩阵

show_slice(axes[0, 0], '[:2, 1:]', (slice(0, 2), slice(1, None)))  # 演示取前两行、从第1列到末尾
show_slice(axes[0, 1], '[2]', (2, slice(None, None)))  # 演示取第2行的所有列
show_slice(axes[1, 0], '[:, :1]', (slice(None, None), slice(0, 1)))  # 演示取所有行的第0列
show_slice(axes[1, 1], '[1, :2]', (1, slice(0, 2)))  # 演示取第1行的前两列

plt.tight_layout(rect=[0, 0.03, 1, 0.95])                   # 在保留总题名顶部空间的同时消除四个切片面板重叠
plt.show()                                                  # 输出四个切片表达式与其行列覆盖范围的对照图
图 3.4: 二维数组切片

3.6 布尔索引

布尔索引(Boolean Indexing)用与数组同形状的真假条件选择元素,例如筛选“值大于阈值”的位置。它避免手工枚举位置,但仍需核对掩码形状、缺失值语义和结果是否为副本。布尔索引从数学角度可视为示性函数(indicator function)的应用。

布尔索引的集合论基础

布尔索引从数学角度是特征函数(Characteristic Function)的应用。

给定: - 源集合 \(A = \{a_1, a_2, ..., a_n\}\) - 谓词 \(P: A \to \{\text{True}, \text{False}\}\)

过滤定义为: \[ S_P = \{a \in A \mid P(a) = \text{True}\} \]

在NumPy中,布尔数组 \(\mathbf{b} \in \{\text{True}, \text{False}\}^n\) 本质上是一个特征函数向量: \[ \mathbf{b}(i) = \begin{cases} 1 & \text{if } A[i] \text{ satisfies } P \\ 0 & \text{otherwise} \end{cases} \]

布尔运算的集合对应

布尔运算 NumPy语法 集合论表示 数学含义
与 (AND) & \(A \cap B\) 交集
或 (OR) \| \(A \cup B\) 并集
非 (NOT) ~ \(A^c\) 补集
异或 (XOR) ^ \((A \cup B) \setminus (A \cap B)\) 对称差

让我们考虑一个例子,我们有一个数组中的数据,以及一个包含重复名字的数组。

# 下面的公司标签均来自长三角,但数值是专为布尔索引机制设计的假设标准化输入
stock_names = np.array(['恒瑞医药', '宁波银行', '宁波港', '恒瑞医药', '宁波港', '宁波银行', '宁波银行'])  # 让公司标签重复出现以演示掩码选行
# 四列依次代表假设的市盈率、市净率、总市值和市销率标准分数,不作任何经验性解释
valuation_data = np.array([  # 构造与七个公司标签逐行对应的假设指标矩阵
    [ 1.39,  1.41,  1.34,  1.41],  # 首行给恒瑞医药四项正标准分,用于命中公司等值掩码
    [-0.39, -0.29, -0.59, -0.35],  # 次行给宁波银行四项负标准分,用于首次命中银行掩码
    [-0.12, -0.60, -0.54, -0.36],  # 第三行对应宁波港,提供与银行不同的负值组合
    [ 1.58,  1.70,  1.79,  1.65],  # 恒瑞医药再次出现,验证布尔索引可返回多行而非首个命中
    [-0.33, -0.75, -0.55, -0.46],  # 宁波港再次出现,使取反掩码也保留重复公司观测
    [-0.96, -0.66, -0.69, -0.82],  # 宁波银行第二次出现,扩充同名标签命中集合
    [-1.18, -0.81, -0.76, -1.07]])  # 末行仍属宁波银行,用于验证三行选择结果保持原顺序

每个公司标签对应 valuation_data 的一行。若要选择所有“恒瑞医药”行,将 stock_names 与字符串比较即可得到布尔数组:

stock_names == '恒瑞医药'  # 生成布尔数组,标记恒瑞医药对应的假设记录
array([ True, False, False,  True, False, False, False])

这个布尔数组可以在索引 valuation_data 数组时传入:

valuation_data[stock_names == '恒瑞医药']  # 用布尔掩码选择恒瑞医药对应的假设记录
array([[1.39, 1.41, 1.34, 1.41],
       [1.58, 1.7 , 1.79, 1.65]])

布尔数组的长度必须与其索引的数组轴的长度相同。你可以将布尔数组与切片或整数混合使用:

# 选择恒瑞医药对应的行,并显示从第2列开始的假设指标
valuation_data[stock_names == '恒瑞医药', 2:]  # 行掩码保留两条公司记录,列边界2:只返回市净率和杠杆率两项指标
array([[1.34, 1.41],
       [1.79, 1.65]])

要选择除了“恒瑞医药”以外的所有行,可以使用 != 或者用 ~ 对条件取反:

valuation_data[~(stock_names == '恒瑞医药')]  # 取反后选择另外两家长三角公司的假设记录
array([[-0.39, -0.29, -0.59, -0.35],
       [-0.12, -0.6 , -0.54, -0.36],
       [-0.33, -0.75, -0.55, -0.46],
       [-0.96, -0.66, -0.69, -0.82],
       [-1.18, -0.81, -0.76, -1.07]])

要组合多个布尔条件,使用布尔算术运算符 & (与) 和 | (或)。Python 关键字 andor 不适用于布尔数组。

mask = (stock_names == '恒瑞医药') | (stock_names == '宁波银行')  # 组合两家公司标签的布尔条件
valuation_data[mask]  # 用组合掩码选择满足任一条件的行
array([[ 1.39,  1.41,  1.34,  1.41],
       [-0.39, -0.29, -0.59, -0.35],
       [ 1.58,  1.7 ,  1.79,  1.65],
       [-0.96, -0.66, -0.69, -0.82],
       [-1.18, -0.81, -0.76, -1.07]])

使用布尔数组设置值时,会将右侧的值或多个值赋到布尔数组为 True 的位置。

# 将估值数据中的所有负值(标准化后低于均值的部分)设置为0
valuation_data[valuation_data < 0] = 0  # 用布尔索引批量修改满足条件的元素
valuation_data  # 核对对恒瑞医药行的原地赋值只改变掩码命中的两行四列
array([[1.39, 1.41, 1.34, 1.41],
       [0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  ],
       [1.58, 1.7 , 1.79, 1.65],
       [0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  ]])

布尔索引必创副本

通过布尔掩码(Boolean Mask)选择数据时,NumPy 总是会创建数据的独立副本。这与切片视图机制截然不同。在进行复杂的数据过滤和清洗时,务必注意这种动态差异对程序总内存占用的影响。

3.7 花式索引 (Fancy Indexing)

花式索引的数学表示

花式索引从数学角度是索引映射(Indexing Map)的应用。

给定:

  • 源数组 \(A \in \mathbb{R}^n\)
  • 有序整数索引序列 \(I=(i_1,i_2,\ldots,i_k)\in\{0,1,\ldots,n-1\}^k\)

\(I\) 不是集合:索引的先后次序决定输出次序,同一索引也可以重复出现。例如,\(I=(2,0,2)\) 会依次取得 \(A[2]\)\(A[0]\)\(A[2]\)

索引映射定义为: \[ \operatorname{fancy}_I(A) = (A[i_1], A[i_2], ..., A[i_k]) \in \mathbb{R}^k. \]

当索引序列 \(I\) 固定时,\(\operatorname{fancy}_I\) 可写成一个由 0—1 选择矩阵表示的线性映射;重复索引对应选择矩阵中的重复行。实际 NumPy 表达式还允许负整数索引,解释前会按轴长度规范化,并始终按给定索引数组的顺序产生结果。花式索引返回副本,而不是基本切片产生的视图。

多维花式索引

对于多维数组 \(A \in \mathbb{R}^{n_0 \times n_1 \times \cdots \times n_{d-1}}\),若每个轴都由整数数组 \(I_0,I_1,\ldots,I_{d-1}\) 索引,NumPy 会先把这些索引数组广播到共同形状 \(B\)。结果的形状就是 \(B\),且对每个广播位置 \(\boldsymbol{b}\)

\[ \operatorname{fancy\_multi}(A,I_0,\ldots,I_{d-1})[\boldsymbol{b}] = A[I_0^{(B)}[\boldsymbol{b}],\ldots,I_{d-1}^{(B)}[\boldsymbol{b}]], \]

其中 \(I_j^{(B)}\) 表示广播到形状 \(B\) 的第 \(j\) 个索引数组。因此,结果是否一维取决于索引数组的广播形状,而不是取决于“每个轴恰有一个整数数组”。

组合索引

花式索引可以与其他索引类型组合。设 \(S\) 为切片对象,\(\mathbf{b}\) 为布尔数组,\(I\) 为整数索引数组:

\[ A[I, S, \mathbf{b}] \]

这种组合提供了极大的灵活性,但需要注意返回的是副本还是视图。

花式索引是 NumPy 采纳的一个术语,用来描述使用整数数组进行索引。假设我们有一个 8 × 4 的数组:

simulated_factors = np.empty((8, 4))                        # 预分配八行四列缓冲区,随后用行号覆盖每个未初始化单元
for i in range(8):  # 每次覆盖一整行,行号决定该行四项共同写入的值
    simulated_factors[i] = i  # 将每行所有元素设为行索引值
simulated_factors  # 展示第i行四项均为i,便于核对后续行花式索引顺序
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 整数:

# 选择第 4, 3, 0, 和 6 行
simulated_factors[[4, 3, 0, 6]]  # 花式索引:按指定顺序选取第4、3、0、6行
array([[4., 4., 4., 4.],
       [3., 3., 3., 3.],
       [0., 0., 0., 0.],
       [6., 6., 6., 6.]])

使用负数索引会从末尾选择行:

# 选择倒数第 3, 5, 和 7 行
simulated_factors[[-3, -5, -7]]  # 负数花式索引:从末尾倒数选取第3、5、7行
array([[5., 5., 5., 5.],
       [3., 3., 3., 3.],
       [1., 1., 1., 1.]])

传递多个一维且形状相同的索引数组时,它们逐位置配对,所以这个特例会返回一维数组;更一般地,多个整数索引数组会先广播,广播形状进入结果形状。

reshaped_factors = np.arange(32).reshape((8, 4))            # 创建0~31的等差序列并重塑为8×4数组
# 选择元素 (1,0), (5,3), (7,1), 和 (2,2)
reshaped_factors[[1, 5, 7, 2], [0, 3, 1, 2]]  # 多索引数组花式索引:选取四个单独元素
array([ 4, 23, 29, 10])

在这个特例中,两个形状均为 (4,) 的索引数组配对选择元素 (1, 0)(5, 3)(7, 1)(2, 2),所以结果形状为 (4,)。下面用可广播的二维索引验证一般规则:行索引形状为 (2, 1),列索引形状为 (1, 3),共同广播形状为 (2, 3)

row_index = np.array([[1], [5]])  # 两个目标行按(2,1)排列,等待沿列轴广播
column_index = np.array([[0, 2, 3]])  # 三个目标列按(1,3)排列,等待沿行轴广播
broadcast_selection = reshaped_factors[row_index, column_index]  # 逐广播位置配对行列坐标,得到2×3元素网格
assert np.broadcast_shapes(row_index.shape, column_index.shape) == (2, 3)  # 核对两个索引数组的共同广播形状
assert broadcast_selection.shape == (2, 3)  # 验证广播形状完整进入高级索引结果形状
np.testing.assert_array_equal(broadcast_selection, reshaped_factors[np.ix_([1, 5], [0, 2, 3])])  # 与显式笛卡尔积选择交叉核验数值
broadcast_selection  # 展示第1、5行与第0、2、3列交叉形成的矩形区域
array([[ 4,  6,  7],
       [20, 22, 23]])

要用一维行列列表选择一个矩形区域,也可以使用以下方法:

# 选择第 1, 5, 7, 2 行,并从这些行中选择第 0, 3, 1, 2 列
reshaped_factors[[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]])

3.8 数组转置和轴交换

转置是一种特殊的重塑形式,它同样返回底层数据的视图而不复制任何内容。数组有 transpose 方法和特殊的 T 属性:

factor_data = np.arange(15).reshape((3, 5))                 # 以3个观测×5个变量组织0—14,便于核对转置后两轴互换为5×3
factor_data.T                                               # 转置矩阵:行变列、列变行,返回视图
array([[ 0,  5, 10],
       [ 1,  6, 11],
       [ 2,  7, 12],
       [ 3,  8, 13],
       [ 4,  9, 14]])

在进行矩阵计算时,这非常常见。例如,当使用 np.dot 计算 \(X^{\mathsf T}X\) 时,结果首先是列向量两两内积组成的 Gram 矩阵(交叉乘积矩阵),不能在未中心化、未缩放时称为协方差矩阵:

# 真实世界数据背景下的例子
# 让我们模拟一个有两个因子的简单回归,所以我们的 factors 矩阵是 n x 2
# 我们将使用随机数据进行演示
rng = np.random.default_rng(seed=12345)                     # 初始化随机数生成器(固定种子保证可复现)
factors = rng.standard_normal((100, 2)) # 生成100×2的标准正态随机因子矩阵
gram_matrix = np.dot(factors.T, factors)  # 计算未中心化的X'X Gram矩阵,保留各列均值对交叉乘积的贡献
gram_matrix  # 展示2×2因子Gram矩阵,行列均对应两项模拟因子
array([[89.36134837, 11.49065471],
       [11.49065471, 93.25390931]])

\(n\) 行观测、\(p\) 列变量的数据矩阵为 \(X\),列均值向量为 \(\bar{x}\),中心化矩阵为 \(X_c=X-\boldsymbol{1}\bar{x}^{\mathsf T}\)。按通常的无偏样本口径,样本协方差矩阵是

\[ S=\frac{1}{n-1}X_c^{\mathsf T}X_c, \tag{3.1}\]

因此只有先按列中心化,再除以 \(n-1\),才得到 式 3.1@ 中缀运算符是进行矩阵乘法的另一种方式:

centered_factors = factors - factors.mean(axis=0)  # 按列减去样本均值,移除均值项对交叉乘积的影响
sample_covariance = centered_factors.T @ centered_factors / (factors.shape[0] - 1)  # 按n-1缩放中心化Gram矩阵得到样本协方差
assert np.allclose(sample_covariance, np.cov(factors, rowvar=False, ddof=1))  # 用NumPy协方差接口核验中心化与自由度口径
sample_covariance  # 展示2×2无偏样本协方差矩阵
array([[0.90262378, 0.1160565 ],
       [0.1160565 , 0.94195154]])

使用 .T 进行简单转置是轴交换的一个特例。ndarrayswapaxes 方法,它接受一对轴号并交换指定的轴来重新排列数据:

tensor_3d = np.arange(16).reshape((2, 2, 4))                # 创建0~15的等差序列并重塑为2×2×4的三维数组
tensor_3d.swapaxes(1, 2)  # 交换轴1和轴2,形状变2×4×2,返回视图
array([[[ 0,  4],
        [ 1,  5],
        [ 2,  6],
        [ 3,  7]],

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

swapaxes 同样返回数据的视图而不进行复制。

3.9 通用函数:逐元素数组函数

通用函数的数学定义

通用函数(ufunc)从数学角度是元素级映射(Element-wise Mapping)的形式化。

给定: - 源数组 \(\mathbf{x} \in \mathbb{R}^n\) - 一元函数 \(f: \mathbb{R} \to \mathbb{R}\)

一元ufunc定义为: \[ \text{ufunc}_f(\mathbf{x})_i = f(x_i), \quad \forall i \in \{0, 1, ..., n-1\} \]

对于二元函数 \(g: \mathbb{R} \times \mathbb{R} \to \mathbb{R}\) 和两个数组 \(\mathbf{x}, \mathbf{y} \in \mathbb{R}^n\)

二元ufunc定义为: \[ \text{ufunc}_g(\mathbf{x}, \mathbf{y})_i = g(x_i, y_i), \quad \forall i \in \{0, 1, ..., n-1\} \]

特殊方法

  1. reduce(归约)\[ \text{reduce}(g, \mathbf{x}) = g(g(...g(g(x_0, x_1), x_2), ..., x_{n-1}) \]

  2. accumulate(累积)\[ \text{accumulate}(g, \mathbf{x})_i = g(x_0, x_1, ..., x_i), \quad \forall i \]

  3. outer(外积)\[ \text{outer}(g, \mathbf{x}, \mathbf{y})_{i,j} = g(x_i, y_j), \quad \forall i, j \]

通用函数(universal function),或称 ufunc,是一种对 ndarray 中的数据执行逐元素操作的函数。你可以将它们看作是简单函数的快速向量化包装器,这些函数接受一个或多个标量值并产生一个或多个标量结果。

许多 ufuncs 是简单的逐元素变换,如 np.sqrtnp.exp

values = np.arange(10)                                      # 构造0—9作为where三种调用方式共享的待筛选值域
np.sqrt(values)                                             # 对每个元素计算平方根(一元ufunc)
array([0.        , 1.        , 1.41421356, 1.73205081, 2.        ,
       2.23606798, 2.44948974, 2.64575131, 2.82842712, 3.        ])

这些被称为一元 ufuncs。其他的,如 np.addnp.maximum,接受两个数组(因此是二元 ufuncs)并返回一个单一数组作为结果:

rng = np.random.default_rng(seed=12345)                     # 初始化随机数生成器(固定种子)
vec_a = rng.standard_normal(8)                              # 生成8个标准正态随机数作为向量a
vec_b = rng.standard_normal(8)                              # 生成8个标准正态随机数作为向量b
np.maximum(vec_a, vec_b)  # 逐元素取两个数组的较大值(二元ufunc)
array([ 0.36105811,  1.26372846,  2.34740965,  0.96849691, -0.07534331,
        0.90219827, -0.46695317,  0.6488928 ])

表 3.3表 3.4 列出了一些 NumPy 的 ufuncs。

表 3.3: 一些一元通用函数。
函数 描述
abs, fabs 计算整数、浮点数或复数的绝对值。对于非复数,fabs 速度更快。
sqrt 计算每个元素的平方根。相当于 arr ** 0.5
square 计算每个元素的平方。相当于 arr ** 2
exp 计算每个元素的指数 \(e^x\)
log, log10, log2, log1p 分别为自然对数(底e)、底10对数、底2对数和 log(1 + x)
sign 计算每个元素的符号:1(正数),0(零),-1(负数)。
ceil 计算每个元素的上取整(大于等于该值的最小整数)。
floor 计算每个元素的下取整(小于等于该值的最大整数)。
rint 将每个元素四舍五入到最近的整数,保留 dtype。
modf 将数组的小数和整数部分作为两个独立的数组返回。
isnan 返回一个布尔数组,指示每个值是否为 NaN(非数字)。
isfinite, isinf 分别返回布尔数组,指示每个元素是有限的(非 inf,非 NaN)还是无限的。
cos, cosh, sin, sinh, tan, tanh 常规和双曲三角函数。
arccos, arccosh, … 反三角函数。
logical_not 对数组中的每个元素计算 not x。相当于 ~arr
表 3.4: 一些二元通用函数。
函数 描述
add 将两个数组中对应的元素相加。
subtract 从第一个数组的元素中减去第二个数组的元素。
multiply 将数组元素相乘。
divide, floor_divide 除法或向下取整的除法(丢弃余数)。
power 将第一个数组中的元素作为底,第二个数组中的元素作为指数,进行幂运算。
maximum, fmax 逐元素计算最大值。fmax 忽略 NaN
minimum, fmin 逐元素计算最小值。fmin 忽略 NaN
mod 逐元素计算模数(余数)。
copysign 将第二个数组中值的符号复制到第一个数组的值上。
greater, ..., not_equal 逐元素比较,产生布尔数组。
logical_and, logical_or, logical_xor 逐元素计算逻辑与、或、异或。

3.9.1 ufunc 的高级用法

NumPy 的每个二元 ufunc 都有特殊的方法来执行特殊的向量化操作。

reduce 接受一个数组,并通过执行一系列二元操作来聚合其值。例如,np.add.reduce 等同于 arr.sum()

values = np.arange(10)                                      # 构造0—9以演示把大于5的命中元素写入预分配输出
np.add.reduce(values)  # 用reduce累加所有元素(等价于sum)
45

accumulatereduce 相关,但它会产生一个与输入大小相同的数组,其中包含中间的“累积”值,等同于 cumsum

matrix_data = np.arange(15).reshape((3, 5))                 # 每行放置5个连续整数,以观察axis=1累积只在行内进行
np.add.accumulate(matrix_data, axis=1)  # 沿轴1计算累积和(每行从左到右累加)
array([[ 0,  1,  3,  6, 10],
       [ 5, 11, 18, 26, 35],
       [10, 21, 33, 46, 60]])

outer 对两个数组执行逐对的叉积。输出数组的形状是输入数组形状的拼接。

vec = np.arange(4)                                          # 创建0~3的整数向量
np.multiply.outer(vec, np.arange(5))  # 计算两个向量的外积(乘法表)
array([[ 0,  0,  0,  0,  0],
       [ 0,  1,  2,  3,  4],
       [ 0,  2,  4,  6,  8],
       [ 0,  3,  6,  9, 12]])

reduceat 执行“局部 reduce”或“分组”操作。它接受一个“分箱边界”序列,指示如何分割和聚合值。

values = np.arange(10)                                      # 构造0—9以比较ufunc.reduce与普通sum的聚合结果
# 对 values[0:5], values[5:8], 和 values[8:] 求和
np.add.reduceat(values, [0, 5, 8])  # 按分箱边界局部求和(分段归约)
array([10, 18, 17])

3.10 使用数组进行面向数组编程

使用 NumPy 数组,你可以把许多同质数值任务表达为数组表达式。这种用数组表达式替代 Python 层显式循环的做法称为向量化

3.10.1 将条件逻辑表达为数组操作

numpy.where 函数是三元表达式 x if condition else y 的向量化版本。

bull_return_scenarios = np.array([1.1, 1.2, 1.3, 1.4, 1.5])  # 创建牛市情景下的收益率数组
bear_return_scenarios = np.array([2.1, 2.2, 2.3, 2.4, 2.5])  # 创建熊市情景下的收益率数组
bull_market_flag = np.array([True, False, True, True, False])  # 创建布尔数组表示各期是否为牛市

假设我们想在 bull_market_flagTrue 时从 bull_return_scenarios 中取值,否则从 bear_return_scenarios 中取值。 …

trading_signal = np.where(bull_market_flag, bull_return_scenarios, bear_return_scenarios)  # 根据牛熊市标志选择对应情景的收益率
trading_signal  # 展示五期结果逐位置来自牛市或熊市场景数组
array([1.1, 2.2, 1.3, 1.4, 2.5])

在数据分析中,where 的一个典型用途是根据另一个数组生成一个新值的数组。例如,假设你有一个随机生成的数据矩阵,你想将所有正值替换为 2,所有负值替换为 –2。

# 模拟收益率矩阵
rng = np.random.default_rng(seed=12345)                     # 初始化随机数生成器
returns_matrix = rng.standard_normal((4, 4))                # 生成4×4的随机收益率矩阵
# 生成方向标记:收益率为正时为 1,为负时为 -1
np.where(returns_matrix > 0, 1, -1)                         # 根据正负号生成二值信号矩阵
array([[-1,  1, -1, -1],
       [-1, -1, -1,  1],
       [ 1, -1,  1,  1],
       [-1,  1, -1, -1]])

3.10.2 数学和统计方法

计算整个数组或沿轴统计量的函数可以作为数组方法访问。你可以调用 summeanstd 等聚合(也称为归约),既可使用实例方法(如 arr.sum()),也可使用顶层函数(如 np.sum(arr))。

# 计算均值和总和
returns_matrix = rng.standard_normal((5, 4))                # 生成5×4的随机收益率矩阵
print(f'returns_matrix 的均值: {returns_matrix.mean()}')  # 计算并输出全部元素的算术平均值
print(f'returns_matrix 的总和: {returns_matrix.sum()}')  # 计算并输出全部元素的总和
returns_matrix 的均值: 0.17933634979615845
returns_matrix 的总和: 3.586726995923169

meansum 这样的函数接受一个可选的 axis 参数,该参数会计算给定轴上的统计量,结果是一个维度减少一的数组:

# 沿着列计算均值 (axis=1)
returns_matrix.mean(axis=1)  # 沿轴1计算均值(每行的平均收益率)
array([ 0.37675318,  0.07598404, -0.2834985 ,  1.48722347, -0.75978045])
# 沿着行计算总和 (axis=0)
returns_matrix.sum(axis=0)  # 沿轴0计算总和(每列的收益率之和)
array([ 2.71870476,  0.30188842, -0.49975489,  1.0658887 ])

其他方法如 cumsumcumprod 不进行聚合,而是产生一个包含中间结果的数组:

matrix_data = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])   # 创建3×3的示例矩阵
# 沿列的累积和
matrix_data.cumsum(axis=1)                                  # 沿轴1计算累积和(每行从左到右逐步累加)
array([[ 0,  1,  3],
       [ 3,  7, 12],
       [ 6, 13, 21]])

表 3.5 列出了一些基本的数组统计方法。

表 3.5: 基本数组统计方法。
方法 描述
sum 数组中所有元素的总和,或沿一个轴的总和;零长度数组的总和为 0。
mean 算术平均值;在零长度数组上无效(返回 NaN)。
std, var 分别为标准差和方差。
min, max 最小值和最大值。
argmin, argmax 分别为最小值和最大值元素的索引。
cumsum 从 0 开始的元素的累积和。
cumprod 从 1 开始的元素的累积积。

3.10.3 布尔数组的方法

在前面的方法中,布尔值被强制转换为 1 (True) 和 0 (False)。因此,sum 经常被用作计算布尔数组中 True 值数量的手段:

pnl_vector = rng.standard_normal(100)                       # 生成100个模拟盈亏值(标准正态分布)
(pnl_vector > 0).sum() # 布尔值按1/0聚合,把逐日盈亏条件压缩为正收益观测数
51

另外两个方法 anyall 对布尔数组很有用。any 测试数组中是否有一个或多个值为 True,而 all 检查是否每个值都为 True

# 模拟每日盈亏状态(已盈利为 True)
profitable_days_mask = np.array([False, False, True, False])  # 创建布尔数组模拟4天的盈亏状态
print(f'是否有盈余天? {profitable_days_mask.any()}')  # 只要一个True即为真,用于回答“样本内是否曾盈利”
print(f'是否所有天都盈利? {profitable_days_mask.all()}')  # 只有四项全为True才为真,检验条件比any更严格
是否有盈余天? True
是否所有天都盈利? False

3.10.4 排序

和 Python 内置的 list 类型一样,NumPy 数组可以使用 sort 方法进行原地排序:

unsorted_returns = rng.standard_normal(6)                   # 生成6个随机收益率值
unsorted_returns.sort()  # 原地升序排序(修改原数组)
unsorted_returns  # 核对原数组已被就地升序排列,未保留排序前顺序
array([-0.7758961 , -0.62361213,  0.28208603,  0.41071644,  0.84122103,
        1.12182226])

你可以通过将轴号传递给 sort,对多维数组中每个一维部分的值沿一个轴进行原地排序:

matrix_data = rng.standard_normal((5, 3))                   # 生成5×3的随机矩阵
matrix_data.sort(axis=1) # 沿轴1原地排序(对每行元素升序排列)
matrix_data  # 展示每一行内部升序而三行之间不交换位置的axis=1结果
array([[-2.7224161 , -0.6733048 ,  1.24622153],
       [-0.0292946 ,  0.17534089,  0.79020803],
       [-1.41951426, -1.35996632,  0.22341156],
       [-2.17088985,  0.62848817,  1.76177943],
       [-0.86924667,  0.60119653,  0.95075786]])

顶层方法 numpy.sort 返回一个排序后的数组副本,而不是在原地修改数组。

3.10.5 唯一值和其他集合逻辑

NumPy 对一维 ndarray 有一些基本的集合操作。一个常用的是 np.unique,它返回数组中排序后的唯一值:

sectors = np.array(['Tech', 'Energy', 'Finance', 'Tech', 'Energy', 'Finance', 'Finance'])  # 创建含重复行业名称的字符串数组
np.unique(sectors)                                          # 返回去重并排序后的唯一值
array(['Energy', 'Finance', 'Tech'], dtype='<U7')

另一个函数 np.in1d,测试一个数组中的值是否在另一个数组中,返回一个布尔数组:

stock_ids = np.array([600000, 0, 0, 600926, 601009, 5, 600000])  # 用长三角银行代码和无效值构造成员测试数组
np.isin(stock_ids, [600926, 601009, 600000])  # 检查每个元素是否属于目标长三角证券池
array([ True, False, False,  True,  True, False,  True])

3.11 使用数组进行文件输入和输出

numpy.savenumpy.load 用于在磁盘上保存和加载数组数据。默认情况下,数组以未压缩的二进制格式保存,文件扩展名为 .npy。本节所有教学写出都放进同一个临时目录;numpy_io_tmpdir 对象会一直存活到本章内核结束,所以后续代码块能够复用路径,同时不会在仓库根目录留下文件:

from pathlib import Path  # 使用路径对象构造同一临时目录中的数组文件名
from tempfile import TemporaryDirectory  # 创建内核结束后自动清理的教学目录
numpy_io_tmpdir = TemporaryDirectory()  # 保持对象存活,使跨代码块回读路径持续有效
numpy_io_root = Path(numpy_io_tmpdir.name)  # 把临时目录包装为跨平台路径对象
npy_array_path = numpy_io_root / 'some_array.npy'  # 固定本节单数组临时文件路径
npz_archive_path = numpy_io_root / 'array_archive.npz'  # 固定本节多数组归档临时路径
saved_data = np.arange(10)                                  # 创建0~9的整数数组作为待保存数据
np.save(npy_array_path, saved_data)  # 把.npy教学文件写入临时目录而非仓库根目录

磁盘上的数组随后可以用 np.load 加载:

loaded_saved_data = np.load(npy_array_path)  # 从上一块的同一临时路径加载数组
np.testing.assert_array_equal(loaded_saved_data, saved_data)  # 核对往返没有改变形状、dtype或元素值
loaded_saved_data  # 展示已通过逐元素断言的回读数组
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

你可以使用 numpy.savez 将多个数组保存在一个未压缩的归档文件中,并将数组作为关键字参数传递:

np.savez(npz_archive_path, a=saved_data, b=saved_data)  # 在临时NPZ中以键a、b保存两份0—9数组

加载 .npz 文件时,你会得到一个类似字典的对象,它会延迟加载各个数组:

with np.load(npz_archive_path) as array_archive:  # 从同一临时路径打开延迟加载归档并确保文件句柄关闭
    archived_b = array_archive['b'].copy()  # 在关闭归档前复制成员b,避免依赖已释放句柄
np.testing.assert_array_equal(archived_b, saved_data)  # 核对归档成员与原始数组完全一致
archived_b  # 展示已通过往返断言的归档成员
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

如果你的数据压缩效果好,你可能希望使用 numpy.savez_compressed

3.12 示例:随机游走

随机游走的模拟为利用数组操作提供了一个说明性的应用。让我们首先考虑一个从 0 开始的简单随机游走,步长为 1 和 –1,且以相等的概率发生。

纯 Python 的实现可能如下所示:

import random                                               # 用标准库逐步抽取二元方向,作为NumPy向量化随机游走的对照
python_walk_rng = random.Random(12345)                      # 用局部固定状态隔离本例,避免依赖全局随机流
position = 0                                                # 初始位置设为原点
walk = [position]                                           # 用列表记录每步位置
nsteps = 1000                                               # 为纯Python基准路径约定1000个独立±1增量
for i in range(nsteps):  # 循环执行1000步随机游走
    step = 1 if python_walk_rng.randint(0, 1) else -1       # 从本例局部随机流等概率选择向前或向后
    position += step                                        # 更新当前位置
    walk.append(position)  # 将新位置追加到轨迹列表

图 3.5 展示了前 100 步的图。

import matplotlib.pyplot as plt                             # 用Matplotlib把纯Python路径的前100个累计位置映射到步序号
plt.plot(walk[:100])                                        # 横轴取列表位置0—99,纵轴取对应累计±1净步数
plt.title('随机游走的前100步')                                     # 限定图中仅展示模拟轨迹的前100次增量
plt.xlabel('步数')                                            # 横轴是离散模拟步序号,而非自然日期
plt.ylabel('位置')                                            # 纵轴是从零点累计的净步数(无量纲)
plt.grid(True)                                              # 用整数步对应的参考线辅助读取累计位置
plt.show()                                                  # 输出前100个离散步的累计位置路径及零点附近波动
图 3.5: 一个简单的随机游走

这个游走是随机步长的累积和。NumPy 可以用一次数组归约表达相同计算;实际耗时仍取决于问题规模和运行环境。

nsteps = 1000                                               # 为NumPy向量化对照使用同样的1000步样本长度
rng = np.random.default_rng(seed=12345)                     # 创建带种子的随机数生成器
draws = rng.integers(0, 2, size=nsteps)                     # 生成1000个0或1的随机整数
steps = np.where(draws == 0, 1, -1)                         # 将0映射为+1步、1映射为-1步
price_path = steps.cumsum()                                 # 对步长求累积和得到位置序列

由此,我们可以提取统计数据,比如游走轨迹中的最小值和最大值:

print(f'最小值: {price_path.min()}')  # 输出游走轨迹中的最小位置
print(f'最大值: {price_path.max()}')  # 输出游走轨迹中的最大位置
最小值: -8
最大值: 50

一个更复杂的统计量是首次穿越时间,即随机游走达到特定值的步数。这里我们可能想知道游走需要多长时间才能距离原点至少 10 步。(np.abs(price_path) >= 10).argmax() 给了我们这个条件首次为真的索引。

(np.abs(price_path) >= 10).argmax()                         # 找到距原点首次超过10步的时刻索引
155

3.12.1 一次模拟多次随机游走

为了模拟多次随机游走,比如 5000 次,我们可以一次性生成所有的随机游走。rng.integers 函数可以生成一个二维的抽样数组,然后我们可以计算每一行的累积和(axis=1)。

nwalks = 5000                                               # 模拟5000条随机游走路径
nsteps = 1000                                               # 每条路径模拟1000步
draws = rng.integers(0, 2, size=(nwalks, nsteps)) # 0 或 1
steps = np.where(draws > 0, 1, -1)                          # 将1映射为+1步、0映射为-1步
simulated_paths = steps.cumsum(axis=1)                      # 沿行累积和得到每条路径的位置序列
simulated_paths  # 展示5000条路径各1000步的累计位置矩阵及轴含义
array([[  1,   2,   3, ...,  22,  23,  22],
       [  1,   0,  -1, ..., -50, -49, -48],
       [  1,   2,   3, ...,  50,  49,  48],
       ...,
       [ -1,  -2,  -1, ..., -10,  -9, -10],
       [ -1,  -2,  -3, ...,   8,   9,   8],
       [ -1,   0,   1, ...,  -4,  -3,  -2]])

现在,我们可以计算所有游走中获得的最大值和最小值:

print(f'所有游走中的最大值: {simulated_paths.max()}')  # 输出5000次游走中的全局最大位置
print(f'所有游走中的最小值: {simulated_paths.min()}')  # 输出5000次游走中的全局最小位置
所有游走中的最大值: 114
所有游走中的最小值: -120

让我们计算达到 30 或 –30 的最小穿越时间。这有点棘手,因为不是所有的游走都会达到 30。我们可以使用 any 方法来检查这一点:

hits30 = (np.abs(simulated_paths) >= 30).any(axis=1)        # 标记每条路径是否曾距原点超过30步
print(f'达到 30 或 -30 的游走次数: {hits30.sum()}')  # 统计满足穿越条件的路径数量
达到 30 或 -30 的游走次数: 3395

我们可以使用这个布尔数组来选择那些实际穿越了绝对 30 水平的游走的行,并在 axis=1 上调用 argmax 来获得穿越时间:

crossing_times = (np.abs(simulated_paths[hits30]) >= 30).argmax(axis=1)  # 对已穿越的路径找出首次距原点超过30步的时刻
print(f'平均最小穿越时间: {crossing_times.mean()}')  # 输出平均首次穿越步数
平均最小穿越时间: 500.5699558173785

3.12.2 一个真实世界的市场路径:上证综指日收益率

随机游走模型是有效市场研究中的基准模型之一,但单条真实市场路径与随机游走“看起来相似”并不能证明收益不可预测。下面绘制本地上证综指的累积对数收益率,只演示如何把NumPy数组工具用于真实数据;正式结论还需要收益自相关、稳健性和样本外检验。

教师提供的数据准备(只运行核对,不评分)

下一代码块只负责把本地 HDF5/Pandas 对象转换为两个本章可直接使用的 NumPy 数组:index_trade_datesindex_close_prices。学生考核从这两个 ndarray 开始,Pandas 读取、筛选、索引和日期解析不在本章评分范围。

import pandas as pd                                         # 用pandas读取、排序并按日期切取本地上证综指行情
import matplotlib.pyplot as plt                             # 用Matplotlib展示样本期累计对数收益路径
import numpy as np                                          # 用NumPy对正收盘价取对数并计算相邻差分
from pathlib import Path                                    # 导入路径处理工具

# 从本地数据获取上证指数数据,并复用首块已验证的跨平台入口
index_path = f'{DATA_ROOT}/index/indexes.h5'               # 复用本章首块的跨平台数据入口定位指数文件

index_data = pd.read_hdf(index_path)                    # 读取本地HDF5格式的指数行情数据
sse_index = index_data[index_data['symbol'] == '000001.XSHG'].copy()  # 筛选上证综指并创建副本
sse_index['datetime'] = pd.to_datetime(sse_index['datetime'].astype(str), format='%Y%m%d%H%M%S')  # 将日期列转换为datetime格式
sse_index.set_index('datetime', inplace=True)               # 将日期设为DataFrame索引
sse_index = sse_index.sort_index()                          # 按日期升序排列

# 选取时间范围
sse_index = sse_index.loc['2014-01-01':'2024-01-01']        # 截取2014年至2024年的数据

# 将真实Pandas对象交付为本章已经讲授的一维NumPy数组
index_trade_dates = sse_index.index.to_numpy(dtype='datetime64[D]')  # 教师把Pandas日期索引交付为一维NumPy日期数组
index_close_prices = sse_index['close'].to_numpy(dtype=np.float64)  # 教师把正收盘点位交付为同长度浮点数组
assert index_trade_dates.shape == index_close_prices.shape  # 核对日期和价格逐项一一对应
assert np.all(index_close_prices > 0)  # 对数变换前核对全部真实指数点位为正

从这里开始只使用本章已经讲授的 NumPy 运算,计算内容属于正文学习范围。

log_returns = np.diff(np.log(index_close_prices))  # 从这里开始仅用已讲授NumPy操作计算相邻日对数收益

# 通过累积和创建“随机游走”
random_walk_sse = np.cumsum(log_returns)  # 累加日对数收益得到区间对数收益;未减基准或无风险收益,不能称为超额收益
# 将2014—2024共同交易日上的累计对数收益映射为日期—收益路径
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC']   # 让上证综指图题和日期轴标签使用项目统一中文字体
plt.rcParams['axes.unicode_minus'] = False                  # 解决负号显示为方块的问题
plt.figure(figsize=(12, 6))                                 # 为十年真实交易日路径分配宽画布,减少日期刻度拥挤
plt.plot(index_trade_dates[1:], random_walk_sse)  # 按真实交易日连接从2014年起累加的上证综指日对数收益
plt.title('作为随机游走的上证指数累积对数收益率')                             # 明确曲线是2014—2024样本内日对数收益的累加
plt.xlabel('日期')                                            # 横轴沿上证指数实际交易日推进
plt.ylabel('累积对数收益率')                                       # 纵轴为无量纲的累计对数收益,不是价格点位
plt.grid(True)                                              # 参考线便于对照不同日期的累计收益水平
plt.show()                                                  # 输出2014—2024实际交易日上的上证综指累计对数收益路径
图 3.6: 上证指数的累积对数收益率 (2014-2024)

图 3.6 展示的是这一真实样本区间的累计对数收益路径。视觉上的相似只说明两条路径都随时间累积波动,不能据此接受随机游走假说或断言收益不可预测。

3.13 习题

请先完成同一难度层级中的题面,再展开完整解答。网页中的解答默认折叠;浏览器打印或导出 PDF 时,下面的打印逻辑会临时展开全部解答,完成后恢复原有开合状态。

3.13.1 基础习题

习题 3.1: 数组创建与基本运算

定位:核心;先修:一维数组创建与聚合;预计时间:15 分钟;评分产出:数值结果、数组元数据与字节偏移解释。

创建两个数组:

import numpy as np  # 使用NumPy构造题设价格数组并完成向量化统计
prices_a = np.array([35.5, 36.2, 35.9, 36.8, 37.1])  # 题设宁波港连续五日收盘价,单位为元
prices_b = np.array([12.3, 12.5, 12.4, 12.8, 12.7])  # 题设宁波银行同一五日口径的收盘价,单位为元
  1. 计算两个数组的价格之和
  2. 计算两个数组价格的均值和标准差
  3. 找出两个数组的最大值和最小值
  4. 对数组进行排序
  5. 计算价格的相关系数
  6. 报告 prices_ashapedtypestrides,并计算索引 3 相对首元素的字节偏移
展开习题 3.1 完整解答

代码清单 列表 3.1 给出本题的可复核解答。

列表 3.1: 习题3.1解答
import numpy as np                                          # 用NumPy完成两只股票五日价格的逐元素统计与相关计算

prices_a = np.array([35.5, 36.2, 35.9, 36.8, 37.1])         # 重建宁波港五日价格,顺序与题设交易日一致
prices_b = np.array([12.3, 12.5, 12.4, 12.8, 12.7])         # 重建宁波银行同五日价格以支持逐位置比较

# (a) 价格之和
total = prices_a + prices_b  # 对应位置价格相加(向量化操作)
print('(a) 价格之和:')  # 标明下方五项数值是两只股票同日价格的逐元素和
print(total)  # 五个位置分别是两条价格序列同日相加值,不是跨日聚合总额

# (b) 均值和标准差
mean_a = np.mean(prices_a)                                  # 计算宁波港均价
std_a = np.std(prices_a)                                    # 计算宁波港价格波动率
mean_b = np.mean(prices_b)                                  # 计算宁波银行均价
std_b = np.std(prices_b)                                    # 计算宁波银行价格波动率
print(f'\n(b) 统计量:')                                        # 区分两只股票各自五日价格的总体均值与总体标准差
print(f'宁波港: 均值={mean_a:.2f}, 标准差={std_a:.2f}')  # 输出宁波港的描述性统计量
print(f'宁波银行: 均值={mean_b:.2f}, 标准差={std_b:.2f}')  # 输出宁波银行的描述性统计量

# (c) 最大值和最小值
print(f'\n(c) 极值:')                                         # 引出每只股票在五个观测日内的最高价与最低价
print(f'宁波港: 最大={np.max(prices_a):.2f}, 最小={np.min(prices_a):.2f}')  # 输出宁波港价格极值
print(f'宁波银行: 最大={np.max(prices_b):.2f}, 最小={np.min(prices_b):.2f}')  # 输出宁波银行价格极值

# (d) 排序
print(f'\n(d) 排序后:')                                        # 标明下方数组按价格升序排列而非按交易日排列
print(f'宁波港: {np.sort(prices_a)}')                          # 宁波港价格升序排列
print(f'宁波银行: {np.sort(prices_b)}')                         # 宁波银行价格升序排列

# (e) 相关系数
corr = np.corrcoef(prices_a, prices_b)[0, 1]                # 计算两只股票价格的Pearson相关系数
print(f'\n(e) 相关系数: {corr:.4f}')                            # 输出相关系数,接近1表示强正相关
(a) 价格之和:
[47.8 48.7 48.3 49.6 49.8]

(b) 统计量:
宁波港: 均值=36.30, 标准差=0.58
宁波银行: 均值=12.54, 标准差=0.19

(c) 极值:
宁波港: 最大=37.10, 最小=35.50
宁波银行: 最大=12.80, 最小=12.30

(d) 排序后:
宁波港: [35.5 35.9 36.2 36.8 37.1]
宁波银行: [12.3 12.4 12.5 12.7 12.8]

(e) 相关系数: 0.9432

代码清单 列表 3.2 给出本题的可复核解答。

列表 3.2: 习题3.1数组元数据解答
element_index = 3  # 选择第四个价格元素用于核对线性地址规则
byte_offset = element_index * prices_a.strides[0]  # 用索引乘字节步长计算相对首元素的偏移
assert byte_offset == prices_a.itemsize * element_index  # 连续一维数组中步长应等于单个元素字节数
print(f'(f) shape={prices_a.shape}, dtype={prices_a.dtype}, strides={prices_a.strides}')  # 报告数组逻辑结构与内存步长
print(f'索引 {element_index} 的相对字节偏移: {byte_offset}')  # 输出可与itemsize手工核对的地址偏移
(f) shape=(5,), dtype=float64, strides=(8,)
索引 3 的相对字节偏移: 24

习题 3.2: 数组索引与切片

定位:核心;先修:基本索引、切片与花式索引;预计时间:20 分钟;评分产出:子矩阵、索引结果及视图/副本原地验证。

给定一个10x10的矩阵:

rng_index_prompt = np.random.default_rng(42)  # 为题面建立独享的固定种子生成器
matrix = rng_index_prompt.integers(1, 100, size=(10, 10))  # 生成待切片和索引的10×10整数矩阵
  1. 提取前3行前3列的子矩阵
  2. 提取所有行第2、4、6列
  3. 提取对角线元素
  4. 将矩阵转置
  5. 找出每行的最大值及其索引
  6. 预测基本切片、花式索引和布尔索引是否共享内存;分别原地修改结果,并验证哪些修改会回写原数组
展开习题 3.2 完整解答

代码清单 列表 3.3 给出本题的可复核解答。

列表 3.3: 习题3.2解答
import numpy as np                                          # 用NumPy复现10×10整数矩阵并比较切片与花式索引

rng_index_solution = np.random.default_rng(42)              # 为答案建立同种子独享生成器,使输入与题面一致
matrix = rng_index_solution.integers(1, 100, size=(10, 10))  # 生成10×10的随机整数矩阵(值域1~99)

# (a) 前3行前3列
sub_matrix = matrix[:3, :3]  # 用切片提取左上角3×3子矩阵
print('(a) 前3行前3列:')  # 标明输出是原矩阵左上角、保留二维形状的3×3视图
print(sub_matrix)  # 输出保持二维3×3形状,并与原矩阵左上角共享缓冲区

# (b) 所有行的第2、4、6列
selected_cols = matrix[:, [1, 3, 5]]  # 用花式索引选取第2、4、6列(索引从0开始)
print(f'\n(b) 第2、4、6列:')                                    # 标明所有10行均保留、列顺序按位置1/3/5抽取
print(selected_cols)  # 十行全部保留,三列按位置1、3、5的给定次序组成副本

# (c) 对角线元素
diagonal = np.diag(matrix)                                  # 提取主对角线元素(即matrix[i,i])
print(f'\n(c) 对角线元素:')                                      # 标明输出依次对应matrix[i,i]的十个主对角元素
print(diagonal)  # 十项依次对应行列位置相同的元素,结果降为一维

# (d) 转置
transposed = matrix.T                                       # 矩阵转置(行变列、列变行)
print(f'\n(d) 转置后形状: {transposed.shape}')                   # 确认形状从(10,10)变为(10,10)

# (e) 每行最大值及索引
max_values = np.max(matrix, axis=1)                         # 沿轴1求每行的最大值
max_indices = np.argmax(matrix, axis=1)                     # 沿轴1求每行最大值所在列索引
print(f'\n(e) 每行最大值:')                                      # 引出每一行最大值及其在该行中的列位置
for i in range(10):  # 用同一行位置配对该行最大值与argmax返回的列位置,防止两组结果错位
    print(f'  行{i}: 值={max_values[i]}, 列={max_indices[i]}')  # 输出每行的最大值及其列位置
(a) 前3行前3列:
[[ 9 77 65]
 [53 97 73]
 [50 37 19]]

(b) 第2、4、6列:
[[77 44 86]
 [97 76 78]
 [37 92 64]
 [23 55  7]
 [76 36 97]
 [20 47  5]
 [74 96 33]
 [19 13 48]
 [44 83 70]
 [80 39 29]]

(c) 对角线元素:
[ 9 97 19 55  7  5 90 23 77 14]

(d) 转置后形状: (10, 10)

(e) 每行最大值:
  行0: 值=86, 列=5
  行1: 值=97, 列=1
  行2: 值=92, 列=3
  行3: 值=88, 列=4
  行4: 值=97, 列=5
  行5: 值=76, 列=0
  行6: 值=96, 列=3
  行7: 值=79, 列=0
  行8: 值=94, 列=0
  行9: 值=89, 列=4

代码清单 列表 3.4 给出本题的可复核解答。

列表 3.4: 习题3.2视图与副本的原地修改核验
memory_test = np.arange(10)  # 构造值与位置一一对应的数组以观察回写效果
slice_view = memory_test[2:5]  # 基本切片通常返回共享底层缓冲区的视图
fancy_copy = memory_test[[5, 6, 7]]  # 整数数组花式索引返回独立副本
boolean_copy = memory_test[memory_test >= 8]  # 布尔索引把命中元素收集到独立副本
slice_view[0] = -2  # 修改视图的首元素,应回写memory_test的索引2
fancy_copy[0] = -5  # 修改花式索引结果,不应改变memory_test的索引5
boolean_copy[0] = -8  # 修改布尔索引结果,不应改变memory_test的索引8
assert np.shares_memory(slice_view, memory_test)  # 验证基本切片与原数组共享内存
assert not np.shares_memory(fancy_copy, memory_test)  # 验证花式索引结果不共享内存
assert not np.shares_memory(boolean_copy, memory_test)  # 验证布尔索引结果不共享内存
assert (memory_test[2], memory_test[5], memory_test[8]) == (-2, 5, 8)  # 同时核对只有视图修改回写原数组
print(memory_test)  # 展示三类修改后的原数组作为可观察证据
[ 0  1 -2  3  4  5  6  7  8  9]

习题 3.3: 布尔索引与条件过滤

定位:核心;先修:布尔掩码与逐元素运算;预计时间:15 分钟;评分产出:筛选结果、计数和复合累计收益路径。

给定一个包含股票收益率的数组:

returns = np.array([0.02, -0.01, 0.03, -0.02, 0.01, 0.04, -0.03, 0.02])  # 构造用于布尔筛选的题设日收益数组
  1. 找出所有正收益
  2. 找出收益率大于2%的数据
  3. 计算正收益和负收益的个数
  4. 将所有负收益替换为0
  5. 计算累计收益率
展开习题 3.3 完整解答

代码清单 列表 3.5 给出本题的可复核解答。

列表 3.5: 习题3.3解答
import numpy as np                                          # 用NumPy对八个题设日收益执行布尔筛选与复利累计

returns = np.array([0.02, -0.01, 0.03, -0.02, 0.01, 0.04, -0.03, 0.02])  # 创建含正负值的股票日收益率数组

# (a) 正收益
positive_mask = returns > 0  # 生成布尔掩码标记收益为正的交易日
positive_returns = returns[positive_mask]  # 用布尔索引提取所有正收益
print('(a) 正收益:')  # 标明输出仅含严格大于零的日收益观测
print(positive_returns)  # 只保留严格大于零的观测,结果长度等于掩码中True的数量

# (b) 收益率大于2%
high_returns = returns[returns > 0.02]  # 筛选收益率超过2%的交易日
print(f'\n(b) 收益率大于2%:')                                    # 标明阈值为严格大于0.02,不包含恰好2%的观测
print(high_returns)  # 严格阈值排除恰好等于2%的观测,并保持原时间顺序

# (c) 正负收益个数
n_positive = np.sum(returns > 0)                            # 统计盈利天数(True计为1)
n_negative = np.sum(returns < 0)                            # 统计亏损天数
print(f'\n(c) 正收益个数: {n_positive}, 负收益个数: {n_negative}')    # 输出盈亏天数统计

# (d) 负收益替换为0
returns_clean = np.where(returns > 0, returns, 0)           # 正收益保留,负收益替换为0
print(f'\n(d) 负收益替换为0后:')                                   # 标明输出保留正收益并把非正观测按题意截为零
print(returns_clean)  # 显示将负收益截断后的结果

# (e) 累计收益率
cumulative = np.cumprod(1 + returns) - 1                    # 计算(1+r)累积乘积再减1得到累计收益率
print(f'\n(e) 累计收益率:')                                      # 标明每项对应截至该日按(1+r)复合的净累计收益
print(cumulative)  # 展示八个观测日逐期复合后的累计收益路径
(a) 正收益:
[0.02 0.03 0.01 0.04 0.02]

(b) 收益率大于2%:
[0.03 0.04]

(c) 正收益个数: 5, 负收益个数: 3

(d) 负收益替换为0后:
[0.02 0.   0.03 0.   0.01 0.04 0.   0.02]

(e) 累计收益率:
[0.02       0.0098     0.040094   0.01929212 0.02948504 0.07066444
 0.03854451 0.0593154 ]

3.13.2 进阶习题

习题 3.4: 矩阵运算与线性代数

定位:核心综合题;先修:矩阵乘法、方差与特征分解;预计时间:25 分钟;评分产出:二次型、标准差、特征对和有边界的经济解释。

在投资组合理论中,矩阵运算是核心工具。给定:

# 3个资产的协方差矩阵
cov_matrix = np.array([
    [0.04, 0.02, 0.01],
    [0.02, 0.09, 0.03],
    [0.01, 0.03, 0.16]
])

# 等权重组合
weights = np.array([1/3, 1/3, 1/3])
  1. 计算投资组合方差:σ²p = w^T Σ w
  2. 计算投资组合标准差
  3. 找出协方差矩阵的特征值和特征向量
  4. 解释特征值的经济学含义
展开习题 3.4 完整解答

代码清单 列表 3.6 给出本题的可复核解答。

列表 3.6: 习题3.4解答
import numpy as np                                          # 用NumPy完成三资产协方差矩阵的二次型与特征分解

cov_matrix = np.array([                                     # 构造三资产收益的对称协方差矩阵
    [0.04, 0.02, 0.01],  # 第一资产方差0.04及其与第二、三资产的协方差
    [0.02, 0.09, 0.03],  # 第二资产方差0.09及其与其余资产的协方差
    [0.01, 0.03, 0.16]  # 第三资产方差0.16并补全对称矩阵
])  # 行列顺序均对应资产1、2、3,数值单位为收益率平方
weights = np.array([1/3, 1/3, 1/3])                         # 按同一资产顺序设置和为1的等权向量

# (a) 投资组合方差
portfolio_variance = np.dot(weights.T, np.dot(cov_matrix, weights))  # 按w转置×协方差矩阵×w计算等权三资产组合方差
print(f'(a) 投资组合方差: {portfolio_variance:.6f}')  # 方差按收益率平方单位保留6位小数

# (b) 投资组合标准差
portfolio_std = np.sqrt(portfolio_variance)                 # 将收益率平方单位的组合方差还原为同频率标准差
print(f'(b) 投资组合标准差: {portfolio_std:.6f} ({portfolio_std*100:.2f}%)')  # 同时报出小数与百分数口径的组合波动率

# (c) 特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)       # numpy线性代数运算
print(f'\n(c) 特征值:')                                        # 引出三个正交风险方向各自解释的方差量
for i, ev in enumerate(eigenvalues):                        # 按NumPy返回次序逐一编号三个特征值
    print(f'  特征值 {i+1}: {ev:.6f}')  # 逐个报告协方差矩阵主方向所解释的方差量

print(f'\n特征向量:')                                           # 引出每个风险方向在三项资产上的单位载荷向量
for i, ev in enumerate(eigenvectors.T):                     # 按列转置后逐一报告与特征值对应的载荷向量
    print(f'  特征向量 {i+1}: {ev}')  # 每行显示相应主方向在三项资产上的载荷

# (d) 经济学含义解释
print(f'\n(d) 经济学含义:')                                      # 将特征分解结果转述为共同风险方向及方差贡献
print(f'特征值代表主成分(Principal Components)解释的方差量。')  # 解释特征值与风险方差分解的对应关系
print(f'最大特征值 {eigenvalues.max():.6f} 对应第一主成分,')  # 标出样本协方差中方差最大的共同方向
print(f'解释了 {eigenvalues.max()/eigenvalues.sum()*100:.1f}% 的总方差。')  # 用最大特征值占迹的比例量化第一主成分贡献
print('该分解可用于描述题设协方差矩阵中的主要共同变化方向。')  # 限定解释对象为题设矩阵而不外推市场事实
(a) 投资组合方差: 0.045556
(b) 投资组合标准差: 0.213437 (21.34%)

(c) 特征值:
  特征值 1: 0.173136
  特征值 2: 0.032982
  特征值 3: 0.083882

特征向量:
  特征向量 1: [-0.12390393 -0.3630554  -0.92349261]
  特征向量 2: [-0.94289213  0.33306872 -0.00443355]
  特征向量 3: [-0.30919612 -0.87020458  0.3835906 ]

(d) 经济学含义:
特征值代表主成分(Principal Components)解释的方差量。
最大特征值 0.173136 对应第一主成分,
解释了 59.7% 的总方差。
该分解可用于描述题设协方差矩阵中的主要共同变化方向。

投资组合风险的矩阵代数表示

在马科维茨现代投资组合理论中,投资组合方差的计算是核心:

\[ \sigma_p^2 = \mathbf{w}^T \Sigma \mathbf{w} = \sum_{i=1}^{n}\sum_{j=1}^{n} w_i w_j \sigma_{ij} \]

其中 \(\mathbf{w}\) 是权重向量,\(\Sigma\) 是协方差矩阵。

谱分解的经济解释: 对实对称协方差矩阵作特征值分解(\(\Sigma = Q \Lambda Q^T\)),可以把题设总方差写成若干正交方向的方差贡献。较大的特征值只表示对应方向在该协方差矩阵中解释的方差较多;若要把方向命名为“市场因子”或据此制定对冲方案,还必须结合真实样本、资产载荷稳定性和样本外检验。


习题 3.5: 广播机制的应用

定位:核心综合题;先修:广播、矩阵乘法与复合收益;预计时间:25 分钟;评分产出:形状推断、组合收益、资产终值及不兼容形状修正。

NumPy的广播机制使得不同形状的数组可以进行运算。给定:

# 每日收益率 (5天 x 3个资产)
daily_returns = np.array([
    [0.01, 0.02, 0.015],
    [-0.01, 0.01, 0.005],
    [0.02, -0.01, 0.01],
    [0.015, 0.005, -0.005],
    [-0.005, 0.01, 0.02]
])

# 3个资产的权重
weights = np.array([0.4, 0.3, 0.3])

在 (a)—(b) 中,weights 表示每日期初的目标权重:假设组合每天在当期收益发生前恢复到 40%、30%、30%,并忽略交易成本。若只在第一天按该权重买入并持有,权重会随各资产收益漂移,累计结果通常不同。(c) 另设每项资产初始投资 100000 元,是独立的买入持有终值计算。

  1. 先使用广播计算每天各资产的加权收益贡献,再沿资产轴求和得到每日投资组合收益率
  2. 在上述每日再平衡假设下计算累计收益率
  3. 独立假设每个资产初始投资100000元并买入持有,计算各资产及合计最终价值
  4. 添加一个新的资产(权重为0),使投资组合包含4个资产
  5. 解释为什么形状为 (5, 3) 的收益矩阵不能与长度为 2 的权重向量逐元素相乘,并给出形状兼容的修正
展开习题 3.5 完整解答

代码清单 列表 3.7 给出本题的可复核解答。

列表 3.7: 习题3.5解答
import numpy as np                                          # 用NumPy按五日×三资产矩阵与标签顺序一致的权重计算组合收益

daily_returns = np.array([                                  # 构建5天×3资产的日收益率矩阵
    [0.01, 0.02, 0.015],  # 第1天:三个资产分别上涨1%、2%、1.5%
    [-0.01, 0.01, 0.005],  # 第2天:资产1下跌,资产2、3上涨
    [0.02, -0.01, 0.01],  # 第3天:资产1大涨,资产2下跌
    [0.015, 0.005, -0.005],  # 第4天:资产3小幅下跌
    [-0.005, 0.01, 0.02]  # 第5天:资产3表现最佳
])  # 每行是一天,每列是一个资产
weights = np.array([0.4, 0.3, 0.3])                         # 每日期初再平衡目标权重依次为40%、30%、30%

# (a) 广播得到逐资产贡献,再沿资产轴聚合为每日组合收益率
weighted_asset_returns = daily_returns * weights  # 形状(5,3)与(3,)广播,保留每天×资产的逐项加权贡献
assert weighted_asset_returns.shape == daily_returns.shape  # 广播本身不收缩资产轴,输出仍为(5,3)
portfolio_returns = weighted_asset_returns.sum(axis=1)  # 沿资产轴求和后,才得到五个交易日的组合收益
matmul_returns = daily_returns @ weights  # 矩阵乘法直接完成内积,作为数值等价但机制不同的交叉核验
np.testing.assert_allclose(portfolio_returns, matmul_returns)  # 核对“广播后求和”与矩阵乘法在本题中数值相同
print('(a) 每日投资组合收益率:')  # 标明输出五项依次对应五个交易日的加权收益
print(portfolio_returns)  # 结果长度为五日,每项由当日三资产行与权重向量内积得到

# (b) 累计收益率
cumulative_returns = np.cumprod(1 + portfolio_returns) - 1  # (1+r)累积乘积再减1得累计收益
print(f'\n(b) 累计收益率:')                                      # 标明输出为组合每日收益逐期复合后的净累计收益
print(cumulative_returns)  # 展示五个交易日结束时组合财富相对起点的净变化
(a) 每日投资组合收益率:
[0.0145 0.0005 0.008  0.006  0.007 ]

(b) 累计收益率:
[0.0145     0.01500725 0.02312731 0.02926607 0.03647093]

这里反复使用同一 weights,所以 cumprod(1 + portfolio_returns) 对应忽略交易成本的每日再平衡策略。若只在期初配置后买入持有,应计算 weights @ np.prod(1 + daily_returns, axis=0) - 1;由于持仓权重随相对价格变化而漂移,该结果通常不同。

代码清单 列表 3.8 继续核对资产终值与零权重资产边界。

列表 3.8: 习题3.5资产终值与零权重边界解答
# (c) 最终价值
initial_value = np.array([100000, 100000, 100000])          # 每个资产初始投资10万元
# 每个资产的价值变化
asset_growth_factors = np.prod(1.0 + daily_returns, axis=0)  # 沿日期轴分别复合三项资产的五日收益
final_values = initial_value * asset_growth_factors  # 用各资产自身增长因子计算三项终值
total_final = np.sum(final_values)  # 汇总三笔各10万元初始投资的期末总价值
print(f'\n(c) 各资产最终价值: {final_values}')  # 输出三项资产分别按自身收益路径得到的终值
print(f'(c) 最终总价值: ¥{total_final:,.2f}')  # 输出30万元初始总额对应的期末总价值

# (d) 添加第4个资产
weights_4 = np.append(weights, 0.0)  # 按题意追加权重为0的第四项资产,原三项权重不变
returns_4 = np.column_stack([daily_returns, np.zeros(5)])  # 拼接第4资产(收益率全为0)
portfolio_returns_4 = returns_4 @ weights_4  # 计算4资产组合的每日加权收益
np.testing.assert_allclose(portfolio_returns_4, portfolio_returns)  # 核对零权重资产不改变原组合收益
print(f'\n(d) 4资产组合每日收益率:')  # 标明零权重资产加入后结果应与三资产组合一致
print(portfolio_returns_4)  # 输出形状仍为(5,)且数值与原组合收益相同

(c) 各资产最终价值: [103002.048765  103524.74649   104562.6858675]
(c) 最终总价值: ¥311,089.48

(d) 4资产组合每日收益率:
[0.0145 0.0005 0.008  0.006  0.007 ]

代码清单 列表 3.9 给出本题的可复核解答。

列表 3.9: 习题3.5广播形状检查解答
incompatible_weights = np.array([0.6, 0.4])  # 构造无法匹配三资产列轴的错误权重
assert daily_returns.shape[-1] != incompatible_weights.shape[0]  # 在运算前识别末轴长度3与权重长度2不兼容
compatible_weights = np.array([0.4, 0.3, 0.3])  # 按三资产列顺序补齐并归一化权重
assert daily_returns.shape[-1] == compatible_weights.shape[0]  # 验证权重长度与收益矩阵资产轴相容
assert np.isclose(compatible_weights.sum(), 1.0)  # 验证组合权重满足全额投资约束
corrected_weighted_returns = daily_returns * compatible_weights  # 广播得到形状(5,3)的逐资产加权贡献
assert corrected_weighted_returns.shape == (5, 3)  # 核对逐元素广播保留日期轴与资产轴
corrected_portfolio_returns = corrected_weighted_returns.sum(axis=1)  # 沿资产轴求和才得到一维五日序列
matmul_check = daily_returns @ compatible_weights  # 矩阵乘法收缩资产轴,作为独立机制的交叉核验
np.testing.assert_allclose(corrected_portfolio_returns, matmul_check)  # 核对两条计算路径在本题中数值一致
assert corrected_portfolio_returns.shape == (5,)  # 核对聚合后的目标结果是一维五日序列

习题 3.6: 本地上证综指收益口径核验

定位:核心综合题;先修:相邻切片、diffcumprodcumsum预计时间:30 分钟;评分产出:三类收益对象、形状断言、恒等式验证和口径解释。

使用 小节 3.12.2 中教师提供且已核对形状的 index_trade_datesindex_close_prices 两个 ndarray 完成下列任务:(a) 仅用 NumPy 计算相邻交易日简单收益与对数收益;(b) 构造起点为 1 的财富指数、累计简单收益和累计对数收益;(c) 验证“累计简单收益 = 财富指数 − 1 = exp(累计对数收益) − 1”;(d) 报告数组首尾日期并解释为什么收益数组比价格数组少一个元素。不得在学生答案中调用 Pandas,也不得把任何一种未经基准扣除的收益称为超额收益。

展开习题 3.6 完整解答

代码清单 列表 3.10 给出本题的可复核解答。

列表 3.10: 习题3.6本地上证综指收益口径解答
assert index_trade_dates.shape == index_close_prices.shape  # 从教师交付合同确认每个日期恰对应一个价格
assert index_close_prices.ndim == 1 and np.all(index_close_prices > 0)  # 核对一维正价格满足相邻比值与对数定义域
simple_index_returns = index_close_prices[1:] / index_close_prices[:-1] - 1.0  # 用相邻切片计算简单收益数组
log_index_returns = np.diff(np.log(index_close_prices))  # 用对数价格的一阶差分计算同区间对数收益
assert simple_index_returns.shape == (index_close_prices.size - 1,)  # 收益区间数应比价格观测数少一
np.testing.assert_allclose(log_index_returns, np.log1p(simple_index_returns))  # 核对两种单期收益口径的精确转换
wealth_index = np.concatenate(([1.0], np.cumprod(1.0 + simple_index_returns)))  # 从1元起点复合得到财富指数
cumulative_simple_return = wealth_index - 1.0  # 财富指数减1得到累计简单收益率
cumulative_log_return = np.concatenate(([0.0], np.cumsum(log_index_returns)))  # 从0起点累加对数收益
np.testing.assert_allclose(cumulative_simple_return, np.exp(cumulative_log_return) - 1.0)  # 核对两条累计收益路径恒等
sample_start_date = np.datetime_as_string(index_trade_dates[0], unit='D')  # 把首个NumPy日期转为ISO日字符串
sample_end_date = np.datetime_as_string(index_trade_dates[-1], unit='D')  # 把末个NumPy日期转为同口径ISO字符串
print(f'样本区间: {sample_start_date}{sample_end_date}')  # 报告教师数组实际首尾日期
print(f'期末财富指数: {wealth_index[-1]:.6f}')  # 展示1元初始财富在样本末的相对价值
print(f'期末累计简单收益: {cumulative_simple_return[-1]:.6%}')  # 展示财富指数减1后的净累计收益
print(f'期末累计对数收益: {cumulative_log_return[-1]:.6f}')  # 展示可跨期相加但未经基准调整的对数收益
样本区间: 2014-01-02 至 2023-12-29
期末财富指数: 1.410331
期末累计简单收益: 41.033144%
期末累计对数收益: 0.343825

收益需要一对相邻价格,因此长度为 \(n\) 的价格数组只能产生 \(n-1\) 个收益区间;这里直接用相邻切片表达该结构,不伪造首期收益。财富指数、累计简单收益和累计对数收益是三个不同对象;它们可以按上述恒等式转换,但都没有扣除市场基准或无风险收益,因此不是超额收益。


3.13.3 可选路线图:随机过程、期权定价与 VaR(不纳入本章考核)

定位:选修路线图;先修:概率分布、随机过程、无套利定价与统计推断;预计时间:20 分钟阅读;产出:能够指出模型测度、参数和损失方向,而不是在本章实现完整定价器。

本章的 cumsumexp、广播和批量抽样可以继续用于金融工程,但数组语法不能替代模型条件。以下三条路线只说明 NumPy 操作与后续课程的接口,题设过程均为假设模型,不是上市公司真实行情,也不构成经验结论。

  • 几何布朗运动:在真实世界测度下,常写为 \(dS_t=\mu S_tdt+\sigma S_tdW_t\)。离散实现会用固定种子生成标准正态增量,再用 cumsumexp 构造路径。\(\mu\)\(\sigma\) 恒定和连续路径只是模型假设,不能由一条模拟路径验证。
  • 欧式期权蒙特卡洛:在 Black–Scholes–Merton 的无套利、无摩擦、连续交易、常数利率与波动率条件下,若股息率 \(q=0\),风险中性测度 \(Q\) 下有 \(dS_t=rS_tdt+\sigma S_tdW_t^Q\),并按 \(C_0=e^{-rT}E^Q[(S_T-K)^+]\) 贴现。这里的 \(E^Q\) 不是按真实世界预期收益率 \(\mu\) 计算的期望。经典推导可追溯到 Black 与 Scholes(1973)Merton(1973)
  • 损失分位数与 VaR:先定义组合损失 \(L=-VR_p\),再讨论尾部。置信水平 \(\alpha\) 下采用损失分布分位数

\[ \operatorname{VaR}_{\alpha}(L)=\inf\{\ell:\Pr(L\leq\ell)\geq\alpha\}. \tag{3.2}\]

式 3.2 明确了正损失方向。若进一步假定 \(R_p\sim N(\mu_p,\sigma_p^2)\),则

\[ \operatorname{VaR}_{\alpha}=V\left(z_{\alpha}\sigma_p-\mu_p\right). \tag{3.3}\]

式 3.3 表明常见的 \(Vz_{\alpha}\sigma_p\) 还采用了 \(\mu_p=0\) 近似。VaR 不描述超过分位点后的平均损失;厚尾、波动聚集、持有期和参数误差都可能改变其解释。

如果后续课程要求实现这些模型,应把终值分布、贴现、标准误、解析解校准和模型失效情形分别验证,而不是把它们扩展成本章第二条主线。

3.14 核心内容回顾

本章从一个四元素数组出发,逐步建立了以下可检验能力:

  • shape 描述轴长度,dtype 描述元素解释规则,strides 把逻辑索引映射为字节偏移;
  • 基本切片通常返回共享缓冲区的视图,整数数组花式索引和布尔索引返回副本;
  • 逐元素运算、广播、矩阵乘法和沿轴聚合解决的问题不同,运行前应先写出输入与输出形状;
  • 简单收益、对数收益、财富指数和累计收益各有明确口径,数组恒等式可以用于核对实现。

随机游走只用于演示批量抽样与累计运算。期权定价和 VaR 还需要测度、无套利条件、损失方向与误差分析,不能由数组语法本身保证正确。后续章节将把这些数组操作接入表格数据、文件格式和时间索引。

章节导航:上一章为内置数据结构、函数与文件,下一章为数据读取与存储:从本地到大数据系统