3  NumPy 基础

NumPy，全称 Numerical Python，是 Python 中用于数值计算最重要的基础包之一。对于经济学和金融学的学生来说，掌握 NumPy 不仅仅是一项技术练习；它是开启现代经济分析中强大定量工具的钥匙，从计量经济建模到金融工程。许多提供科学计算功能的包，如 pandas 和 statsmodels，都使用 NumPy 的数组对象作为数据交换的标准“通用语言”之一。我们在 sec-numpy-basics 中涵盖的关于 NumPy 的大部分知识也同样适用于这些更高级的库。

以下是您将在 NumPy 中找到的一些核心功能：

虽然 NumPy 本身不提供高级的数据建模功能，但对其数组和面向数组的计算范式的深刻理解，是高效使用像 pandas 这类工具的先决条件。

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

作为经济学专业的学生，你可能会问：为什么我们要从如此底层的库学起？答案很简单：NumPy是整个Python数据科学生态系统的基石。你将来要用到的几乎所有计量经济学、金融建模和机器学习库（例如 pandas, statsmodels, scikit-learn, tensorflow）的底层数据结构都是NumPy数组。理解NumPy的运行机制——特别是其对内存的有效管理和向量化计算的理念——能够让你写出运行速度快几个数量级的代码。在处理大规模经济数据（如面板数据、高频交易数据）时，这种效率的差异是决定性的。

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

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

列表 3.1 
import numpy as np

my_arr = np.arange(1_000_000)
my_list = list(range(1_000_000))

print('NumPy 数组乘以 2:')
%timeit my_arr2 = my_arr * 2

print('\nPython 列表乘以 2:')
%timeit my_list2 = [x * 2 for x in my_list]

正如 列表 lst-numpy-performance 中的计时结果所示，基于 NumPy 的算法通常比纯 Python 的实现快 10 到 100 倍（甚至更多），并且使用的内存也显著减少。

3.1 NumPy ndarray：一个多维数组对象

NumPy 的一个关键特性是其 N 维数组对象，即 ndarray。它是一个快速、灵活的容器，适用于 Python 中的大型数据集，允许您使用类似于标量元素之间等效运算的语法，对整块数据执行数学运算。

3.1.1 ndarray 对象内部机制

要真正理解 ndarray 为何如此高效，我们必须超越其 Python 接口，深入其内部。NumPy ndarray 提供了一种将一块同质类型的数据（连续或跨步的内存块）解释为多维数组对象的方式。

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

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

fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10)
ax.set_ylim(0, 7)
ax.axis('off')

# Main ndarray object
ax.add_patch(patches.FancyBboxPatch((1, 4.5), 8, 2, boxstyle='round,pad=0.1', fc='lightblue', ec='black'))
ax.text(5, 5.5, 'NumPy ndarray Object', ha='center', va='center', fontsize=16, fontweight='bold')

# Data Pointer
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'))

# Dtype
ax.add_patch(patches.Rectangle((4, 2.5), 2, 1, fc='lightcoral', ec='black'))
ax.text(5, 3, 'dtype: float64', ha='center', va='center', fontsize=12)
ax.add_patch(patches.Arrow(3.5, 4.5, 1, -0.8, width=0.2, color='gray'))


# Shape
ax.add_patch(patches.Rectangle((7, 2.5), 2, 1, fc='lightgoldenrodyellow', ec='black'))
ax.text(8, 3, 'shape: (10, 5)', ha='center', va='center', fontsize=12)
ax.add_patch(patches.Arrow(6.5, 4.5, 1, -0.8, width=0.2, color='gray'))

# Strides
ax.add_patch(patches.Rectangle((7, 0.5), 2, 1, fc='plum', ec='black'))
ax.text(8, 1, 'strides: (40, 8)', ha='center', va='center', fontsize=12)
ax.add_patch(patches.Arrow(5, 4.5, 2.5, -2.8, width=0.2, color='gray'))

ax.set_title('ndarray 内部结构示意图', fontsize=18, pad=20)
plt.show()
图 3.1

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

np.ones((10, 5)).shape

一个典型的 C 顺序（行主序）的 3 × 4 × 5 的 float64（8字节）数组，其跨度为 (160, 40, 8)。这个信息对性能至关重要，因为特定轴上的跨度越大，意味着沿着该轴进行计算的成本越高。

np.ones((3, 4, 5), dtype=np.float64).strides

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

3.1.2 创建 ndarray

为了让您体验 NumPy 如何实现批量计算，我们首先导入 NumPy 并创建一个小数组。我们将遵循标准约定 import numpy as np。

import numpy as np

# 从一个列表的列表创建一个二维数组
data = np.array([[1.5, -0.1, 3], [0, -3, 6.5]])
data

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

# 逐元素乘法
data * 10
# 逐元素加法
data + data

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

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

# 获取数组的形状
data.shape
# 获取数组的数据类型
data.dtype

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

# 从列表创建
data1 = [6, 7.5, 8, 0, 1]
arr1 = np.array(data1)
arr1

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

data2 = [[1, 2, 3, 4], [5, 6, 7, 8]]
arr2 = np.array(data2)
arr2

由于 data2 是一个列表的列表，NumPy 数组 arr2 有两个维度。我们可以通过检查 ndim 和 shape 属性来确认这一点：

arr2.ndim
arr2.shape

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

# 创建一个长度为10的全零一维数组
np.zeros(10)
# 创建一个3x6的全零数组
np.zeros((3, 6))
# 创建一个未初始化的2x3x2数组
np.empty((2, 3, 2))
关于 np.empty 的警告

假设 numpy.empty 会返回一个全零数组是不安全的。此函数返回的是未初始化的内存，因此可能包含任意的“垃圾”值。你应当只在确定会立即用新数据填充该数组的情况下使用此函数，以避免潜在的数值错误。

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

np.arange(15)

表 tbl-array-creation 列出了一些标准的数组创建函数。

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

3.1.3 ndarray 的数据类型

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

arr1 = np.array([6, 7.5, 8, 0, 1])
arr2 = np.array([[1, 2, 3], [4, 5, 6]])
print(f'arr1 的 dtype: {arr1.dtype}')
print(f'arr2 的 dtype: {arr2.dtype}')

因为 arr1 包含浮点数，其 dtype 被推断为 float64。arr2 只包含整数，被推断为 int64（在64位系统上）。

您可以在创建时显式指定 dtype：

arr1 = np.array([1, 2, 3], dtype=np.float64)
arr2 = np.array([1, 2, 3], dtype=np.int32)
print(f'arr1 的 dtype: {arr1.dtype}')
print(f'arr2 的 dtype: {arr2.dtype}')

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

表 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 ? 存储 True 和 False 值的布尔类型
object O Python 对象类型
string_ S 固定长度 ASCII 字符串类型 (例如 S10 代表 10 字节字符串)
unicode_ U 固定长度 Unicode 类型 (例如 U10 代表 10 字符字符串)
有符号 vs. 无符号整数 (Signed vs. Unsigned Integers)

在经济和金融数据中，负数是很常见的（例如，亏损、负增长率）。因此，你需要了解signed和unsigned整数类型的区别。signed整数可以表示正数和负数，而unsigned整数只能表示非负数。例如，int8（有符号8位整数）的范围是-128到127，而uint8（无符号8位整数）的范围是0到255。错误地使用unsigned类型可能会导致负数被错误地解释为一个非常大的正数，从而引发严重的计算错误。

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

arr = np.array([1, 2, 3, 4, 5])
print(f'原始 dtype: {arr.dtype}')

float_arr = arr.astype(np.float64)
print(f'新的 dtype: {float_arr.dtype}')
print(f'新的数组: {float_arr}')

如果您将浮点数转换为整数 dtype，小数部分将被截断：

arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1])
arr.astype(np.int32)

如果您有一个由表示数字的字符串组成的数组，您可以使用 astype 将它们转换为数值形式：

numeric_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_)
numeric_strings.astype(float)
astype 总是创建一个新数组

调用 astype 总是会创建一个新的数组（数据的副本），即使新的 dtype 与旧的 dtype 相同。这确保了原始数组的数据不会被意外修改，但在处理非常大的数组时需要注意内存使用。

3.1.4 NumPy 数据类型层次结构

您可能偶尔需要检查一个数组的 dtype 是否属于某个通用类别，比如整数或浮点数。数据类型有诸如 np.integer 和 np.floating 这样的超类，可以与 np.issubdtype 函数一起使用：

ints = np.ones(10, dtype=np.uint16)
floats = np.ones(10, dtype=np.float32)

print(f'ints.dtype 是 np.integer 的子类型吗? {np.issubdtype(ints.dtype, np.integer)}')
print(f'floats.dtype 是 np.floating 的子类型吗? {np.issubdtype(floats.dtype, np.floating)}')

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

np.float64.mro()
import matplotlib.pyplot as plt

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'))
        if isinstance(value, dict):
            child_keys = list(value.keys())
            num_children = len(child_keys)
            child_x_start = x - (num_children - 1) * dx / 2
            for i, child_key in enumerate(child_keys):
                child_x = child_x_start + i * dx
                child_y = 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
            for i, child_item in enumerate(value):
                 child_x = child_x_start + i * dx
                 child_y = y - dy
                 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'))


hierarchy = {
    'np.generic': {
        'np.number': {
            'np.integer': {
                'np.signedinteger': ['int8', 'int16', 'int32', 'int64'],
                'np.unsignedinteger': ['uint8', 'uint16', 'uint32', 'uint64']
            },
            'np.floating': ['float16', 'float32', 'float64'],
            'np.complexfloating': ['complex64', 'complex128']
        },
        'np.bool_': {},
        'np.character': {
            'np.string_': [],
            'np.unicode_': []
        }
    }
}

fig, ax = plt.subplots(figsize=(16, 10))
ax.axis('off')
plot_hierarchy(ax, hierarchy, 0.5, 1.0, 0.4, 0.15)
ax.set_title('NumPy Data Type Hierarchy', fontsize=20)
plt.show()
图 3.2

3.2 NumPy 数组的运算

数组的重要性在于，它们使您能够表达对数据的批量操作，而无需编写任何 for 循环。NumPy 用户称之为向量化（vectorization）。任何在等大小数组之间的算术运算都将应用该运算于逐个元素。

arr = np.array([[1., 2., 3.], [4., 5., 6.]])
# 逐元素相乘
arr * arr
# 逐元素相减
arr - arr

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

核心概念：广播 (Broadcasting)

广播描述了 NumPy 在算术运算期间如何处理不同形状的数组。在满足特定规则的情况下，较小的数组会被“广播”到较大数组的形状，以便它们兼容。最简单的情况是将一个标量（可以看作是零维数组）与一个数组相结合。在这种情况下，NumPy 会将标量值“拉伸”或“复制”，使其形状与另一个数组匹配。例如，在 arr * 2 中，2 被广播到 arr 的形状，使得每个元素都乘以2。这个机制是 NumPy 向量化运算效率和灵活性的关键。

# 标量除法
1 / arr
# 标量指数运算
arr ** 2

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

arr2 = np.array([[0., 4., 1.], [7., 2., 12.]])
arr2 > arr

3.3 基本索引和切片

NumPy 数组索引是一个深度话题，有很多方法可以选择数据的子集。一维数组很简单；它们的操作类似于 Python 列表。

arr = np.arange(10)
arr
# 获取索引为 5 的元素
arr[5]
# 获取从索引 5 到 8 (不包括) 的切片
arr[5:8]

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

arr[5:8] = 12
arr
关键区别：视图（Views） vs. 副本（Copies）

与Python的内置列表不同，NumPy中的数组切片是原始数组的视图。这意味着数据没有被复制，对视图的任何修改都会直接反映在源数组上。这个设计是为了处理非常大的数组，避免因不必要的复制而导致的性能和内存问题。

如果你需要一个切片的副本而不是视图，你必须显式地进行复制，例如 arr[5:8].copy()。

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

arr_slice = arr[5:8]
arr_slice

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

arr_slice[:] = 12345
arr

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

arr_slice[:] = 64
arr

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

arr2d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# 访问第三行 (索引 2)
arr2d[2]

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

# 这两种形式是等价的
print(f'递归访问: {arr2d[0][2]}')
print(f'元组访问: {arr2d[0, 2]}')

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

import matplotlib.pyplot as plt

def plot_indexing_diagram():
    fig, ax = plt.subplots(figsize=(4, 4))
    ax.set_xlim(-0.5, 2.5)
    ax.set_ylim(-0.5, 2.5)
    ax.invert_yaxis()
    
    for i in range(3):
        for j in range(3):
            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))
    ax.set_yticks(range(3))
    ax.xaxis.tick_top()
    ax.xaxis.set_label_position('top')
    ax.set_xlabel('轴 1 (列)', fontsize=12)
    ax.set_ylabel('轴 0 (行)', fontsize=12)
    
    # Highlight one element
    ax.add_patch(plt.Rectangle((2-0.5, 0-0.5), 1, 1, fill=True, color='tomato', alpha=0.8, ec='black'))
    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()

plot_indexing_diagram()
图 3.3

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

arr3d = np.array([[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]])
# 在轴 0 上切片返回一个 2x3 数组
arr3d[0]

3.3.1 使用切片进行索引

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

# 选择 arr2d 的前两行
arr2d[:2]

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

# 选择前两行以及从索引 1 开始的所有列
arr2d[:2, 1:]

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

lower_dim_slice = arr2d[1, :2]
lower_dim_slice

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

# 选择所有行中的第一列
arr2d[:, :1]

图 fig-2d-slicing 提供了切片的可视化总结。

import numpy as np
import matplotlib.pyplot as plt

def show_slice(ax, slice_str, highlight):
    arr = np.arange(25).reshape((5, 5))
    ax.imshow(np.zeros_like(arr), cmap='Pastel2', vmin=-1, vmax=1)
    
    # Highlight the slice
    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)
    highlight_mask[rows, cols] = True
    ax.imshow(np.where(highlight_mask, 1, 0), cmap='Set1', alpha=0.5, vmin=-1, vmax=1)

    for i in range(5):
        for j in range(5):
            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([])
    ax.set_yticks([])

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
fig.suptitle('二维数组切片示例', fontsize=20)

show_slice(axes[0, 0], '[:2, 1:]', (slice(0, 2), slice(1, None)))
show_slice(axes[0, 1], '[2]', (2, slice(None, None)))
show_slice(axes[1, 0], '[:, :1]', (slice(None, None), slice(0, 1)))
show_slice(axes[1, 1], '[1, :2]', (1, slice(0, 2)))

plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
图 3.4

3.4 布尔索引

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

names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
data = np.random.standard_normal((7, 4)) # 使用随机数据进行说明

假设每个名字对应 data 数组中的一行，我们想选择所有名字为 ‘Bob’ 的行。像算术运算一样，比较也是向量化的。将 names 与字符串 ‘Bob’ 比较会产生一个布尔数组：

names == 'Bob'

这个布尔数组可以在索引 data 数组时传入：

data[names == 'Bob']

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

# 选择名字为 'Bob' 的行，并显示从第 2 列开始的列
data[names == 'Bob', 2:]

要选择除了 ‘Bob’ 以外的所有内容，你可以使用 != 或者用 ~ 对条件取反：

data[~(names == 'Bob')]

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

mask = (names == 'Bob') | (names == 'Will')
data[mask]

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

# 将 data 中的所有负值设置为 0
data[data < 0] = 0
data
布尔索引总是创建副本

通过布尔索引选择数据总是会创建数据的副本，即使返回的数组与原数组没有任何变化。这与切片操作返回视图的行为形成鲜明对比，是NumPy中一个需要特别注意的重要细节。

3.5 花式索引 (Fancy Indexing)

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

arr = np.empty((8, 4))
for i in range(8):
    arr[i] = i
arr

要以特定顺序选择行的子集，你可以传递一个指定所需顺序的列表或 ndarray 整数：

# 选择第 4, 3, 0, 和 6 行
arr[[4, 3, 0, 6]]

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

# 选择倒数第 3, 5, 和 7 行
arr[[-3, -5, -7]]

传递多个索引数组会做一些稍微不同的事情；它会选择一个一维数组，其元素对应于每个索引元组。

arr = np.arange(32).reshape((8, 4))
# 选择元素 (1,0), (5,3), (7,1), 和 (2,2)
arr[[1, 5, 7, 2], [0, 3, 1, 2]]

在这种情况下，花式索引的行为是选择元素 (1, 0)、(5, 3)、(7, 1) 和 (2, 2)。当传递的整数数组数量与轴数相同时，花式索引的结果总是一维的。要选择一个矩形区域，你可以使用以下方法：

# 选择第 1, 5, 7, 2 行，并从这些行中选择第 0, 3, 1, 2 列
arr[[1, 5, 7, 2]][:, [0, 3, 1, 2]]

3.6 数组转置和轴交换

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

arr = np.arange(15).reshape((3, 5))
arr.T

在进行矩阵计算时，这非常常见。例如，当使用 np.dot 计算内积 \(\boldsymbol{X'X}\) 时：

# 真实世界数据背景下的例子
# 让我们模拟一个有两个变量的简单回归，所以我们的 X 矩阵是 n x 2
# 我们将使用随机数据进行演示
rng = np.random.default_rng(seed=12345)
X = rng.standard_normal((100, 2)) # 100 个观测值, 2 个变量
X_prime_X = np.dot(X.T, X)
X_prime_X

@ 中缀运算符是进行矩阵乘法的另一种方式：

X.T @ X

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

arr = np.arange(16).reshape((2, 2, 4))
arr.swapaxes(1, 2)

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

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

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

许多 ufuncs 是简单的逐元素变换，如 np.sqrt 或 np.exp：

arr = np.arange(10)
np.sqrt(arr)

这些被称为一元 ufuncs。其他的，如 np.add 或 np.maximum，接受两个数组（因此是二元 ufuncs）并返回一个单一数组作为结果：

rng = np.random.default_rng(seed=12345)
x = rng.standard_normal(8)
y = rng.standard_normal(8)
np.maximum(x, y)

表 tbl-unary-ufuncs 和 表 tbl-binary-ufuncs 列出了一些 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.7.1 ufunc 的高级用法

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

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

arr = np.arange(10)
np.add.reduce(arr)

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

arr = np.arange(15).reshape((3, 5))
np.add.accumulate(arr, axis=1)

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

arr = np.arange(4)
np.multiply.outer(arr, np.arange(5))

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

arr = np.arange(10)
# 对 arr[0:5], arr[5:8], 和 arr[8:] 求和
np.add.reduceat(arr, [0, 5, 8])

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

使用 NumPy 数组，您可以将许多类型的数据处理任务表达为简洁的数组表达式，而这些任务可能需要编写循环。这种用数组表达式替代显式循环的做法被称为向量化。

3.8.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 为 True 时从 xarr 中取值，否则从 yarr 中取值。列表推导式会是： result = [(x if c else y) for x, y, c in zip(xarr, yarr, cond)]

对于大数组来说，这很慢，并且不适用于多维数组。使用 np.where 你可以简洁地完成这个任务：

result = np.where(cond, xarr, yarr)
result

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

rng = np.random.default_rng(seed=12345)
arr = rng.standard_normal((4, 4))
np.where(arr > 0, 2, -2)

3.8.2 数学和统计方法

一组计算整个数组或沿轴的数据统计量的数学函数，可以作为数组类的方法来访问。您可以使用聚合（有时称为归约）如 sum、mean 和 std（标准差），通过调用数组实例方法（例如 arr.sum()）或使用顶层 NumPy 函数（例如 np.sum(arr)）。

arr = rng.standard_normal((5, 4))
print(f'arr 的均值: {arr.mean()}')
print(f'arr 的总和: {arr.sum()}')

像 mean 和 sum 这样的函数接受一个可选的 axis 参数，该参数会计算给定轴上的统计量，结果是一个维度减少一的数组：

# 沿着列计算均值 (axis=1)
arr.mean(axis=1)
# 沿着行计算总和 (axis=0)
arr.sum(axis=0)

其他方法如 cumsum 和 cumprod 不进行聚合，而是产生一个包含中间结果的数组：

arr = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
# 沿列的累积和
arr.cumsum(axis=1)

表 tbl-stat-methods 列出了一些基本的数组统计方法。

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

3.8.3 布尔数组的方法

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

arr = rng.standard_normal(100)
(arr > 0).sum() # 正值的数量

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

bools = np.array([False, False, True, False])
print(f'有 True 吗? {bools.any()}')
print(f'都是 True 吗? {bools.all()}')

3.8.4 排序

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

arr = rng.standard_normal(6)
arr.sort()
arr

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

arr = rng.standard_normal((5, 3))
arr.sort(axis=1) # 对每一行进行排序
arr

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

3.8.5 唯一值和其他集合逻辑

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

names = np.array(['Bob', 'Will', 'Joe', 'Bob', 'Will', 'Joe', 'Joe'])
np.unique(names)

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

values = np.array([6, 0, 0, 3, 2, 5, 6])
np.in1d(values, [2, 3, 6])

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

numpy.save 和 numpy.load 是在磁盘上高效保存和加载数组数据的两个主力函数。默认情况下，数组以未压缩的原始二进制格式保存，文件扩展名为 .npy：

arr = np.arange(10)
np.save('some_array', arr)

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

np.load('some_array.npy')

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

np.savez('array_archive.npz', a=arr, b=arr)

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

arch = np.load('array_archive.npz')
arch['b']

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

3.10 示例：随机游走

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

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

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

图 fig-random-walk 展示了前 100 步的图。

import matplotlib.pyplot as plt
plt.plot(walk[:100])
plt.title('随机游走的前100步')
plt.xlabel('步数')
plt.ylabel('位置')
plt.grid(True)
plt.show()
图 3.5

这个游走是随机步长的累积和。我们可以使用 NumPy 更高效地计算它。

nsteps = 1000
rng = np.random.default_rng(seed=12345)
draws = rng.integers(0, 2, size=nsteps)
steps = np.where(draws == 0, 1, -1)
walk = steps.cumsum()

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

print(f'最小值: {walk.min()}')
print(f'最大值: {walk.max()}')

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

(np.abs(walk) >= 10).argmax()

3.10.1 一次模拟多次随机游走

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

nwalks = 5000
nsteps = 1000
draws = rng.integers(0, 2, size=(nwalks, nsteps)) # 0 或 1
steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(axis=1)
walks

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

print(f'所有游走中的最大值: {walks.max()}')
print(f'所有游走中的最小值: {walks.min()}')

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

hits30 = (np.abs(walks) >= 30).any(axis=1)
print(f'达到 30 或 -30 的游走次数: {hits30.sum()}')

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

crossing_times = (np.abs(walks[hits30]) >= 30).argmax(axis=1)
print(f'平均最小穿越时间: {crossing_times.mean()}')

3.10.2 一个真实世界的“随机游走”：标普500指数日收益率

虽然随机游走模型是一个理论构建，但它构成了金融学中有效市场假说的基础。该假说认为，资产价格反映了所有可用信息，因此价格变化应该遵循“随机游走”，是不可预测的。让我们通过绘制标普500指数的累积对数收益率，来看看这如何应用于真实数据。

import pandas as pd
from fredapi import Fred
import matplotlib.pyplot as plt
import numpy as np

# 使用在提示中提供的 API 密钥
fred_key = 'f2a2c60b6dc82682031f4ce84bf6da18'
fred = Fred(api_key=fred_key)

# 获取标普500数据
sp500 = fred.get_series('SP500', start_date='2014-01-01', end_date='2024-01-01')

# 计算对数收益率
log_returns = np.log(sp500 / sp500.shift(1)).dropna()

# 通过累积和创建“随机游走”
random_walk_sp500 = log_returns.cumsum()

# 绘图
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用于显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
plt.figure(figsize=(12, 6))
plt.plot(random_walk_sp500)
plt.title('作为随机游走的标普500累积对数收益率')
plt.xlabel('日期')
plt.ylabel('累积对数收益率')
plt.grid(True)
plt.show()
图 3.6

正如您在 图 fig-sp500-walk 中所看到的，标普500指数累积收益率的路径与我们模拟的随机游走相似。这展示了在 NumPy 中开发的强大、抽象的工具如何能被直接应用于分析和理解复杂的经济和金融现象。

3.11 结论

本章介绍了 NumPy 中的基本对象 ndarray，以及它所支持的面向数组的编程风格。虽然本课程的其余大部分内容将侧重于使用 pandas 构建数据整理技能，但我们将继续以类似的基于数组的风格工作。理解 NumPy 的核心概念对于在 Python 中进行有效的定量分析至关重要。