一、用于数组的文件输入和输出

NumPy能够读写磁盘上的文本数据或二进制数据。这里只讨论NumPy的内置二进制格式。

np.save和np.load是读写磁盘数组数据的两个主要函数。默认情况下，数组是以未压缩的原始二进制格式保存在扩展名为**.npy**的文件中的：

1
2
3
4
5
6
import numpy as np
arr = np.arange(10)
np.save('some_array', arr) # 如果文件路径末尾没有扩展名.npy，则该扩展名会被自动加上
print(np.load('some_array.npy'))
# 结果
[0 1 2 3 4 5 6 7 8 9]

通过np.savez可以将多个数组保存到一个未压缩文件中，将数组以关键字参数的形式传入即可，然后加载.npz文件时，你会得到一个类似字典的对象，该对象会对各个数组进行延迟加载：：

1
2
3
4
5
6
7
8
···
arr1 = np.arange(11)
np.savez('array_archive.npz', a=arr, b=arr1)
arch = np.load("array_archive.npz")
print(arch['b'])

# 结果
[0 1 2 3 4 5 6 7 8 9 10]

如果要将数据压缩，可以使用numpy.savez_compressed：

1
2
3
4
5
6
7
···
np.savez_compressed('arrays_compressed.npz', a=arr, b=arr1)
arch = np.load("array_archive.npz")
print(arch['b'])

# 结果
[0 1 2 3 4 5 6 7 8 9 10]

二、线性代数

2.1 矩阵乘法

线性代数（如矩阵乘法、矩阵分解、行列式以及其他方阵数学等）是任何数组库的重要组成部分。不像某些语言（如MATLAB），通过*对两个二维数组相乘得到的是一个元素级的积，而不是一个矩阵点积（内积）。因此，NumPy提供了一个用于矩阵乘法的dot函数（既是一个数组方法也是numpy命名空间中的一个函数）：

1
2
3
4
5
6
7
8
9
import numpy as np
x = np.array([[1., 2., 3.], [4., 5., 6.]]) # 2,3
y = np.array([[6., 23.], [-1, 7], [8, 9]]) # 3,2
print(x.dot(y))
print(np.dot(x, y)) # x.dot(y)等价于np.dot(x, y)

# 结果
[[ 28.  64.]
 [ 67. 181.]]

一个二维数组跟一个大小合适的一维数组的矩阵点积运算之后将会得到一个一维数组：

1
2
3
4
5
···
print(np.dot(x, np.ones(3))) # np.ones(3) -> 3，1

# 结果
[ 6. 15.]

@符（类似Python 3.5）也可以用作中缀运算符，进行矩阵乘法：

1
2
3
4
5
···
print(x @ np.ones(3))

# 结果
[ 6. 15.]

2.2 Numpy.linalg

numpy.linalg中有一组标准的矩阵分解运算以及诸如求逆和行列式之类的东西。它们跟MATLAB和R等语言所使用的是相同的行业标准线性代数库，如BLAS、LAPACK、Intel MKL（Math Kernel Library，可能有，取决于你的NumPy版本）等：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
from numpy.linalg import inv, qr
import numpy as np
X = np.random.randn(5, 5)
mat = X.T.dot(X)

print("-----inv(mat)-----")
print(inv(mat))
print('-----mat.dot(inv(mat))')
print(mat.dot(inv(mat)))
print("-----qr(mat)-----")
q,r = qr(mat)
print(r)

# 结果
-----inv(mat)-----
[[23.57378102 25.59614691 -1.14247588 -6.16387588  2.59455459]
 [25.59614691 28.03184927 -1.2771725  -6.70868113  2.82896575]
 [-1.14247588 -1.2771725   0.22263419  0.41629033 -0.19599961]
 [-6.16387588 -6.70868113  0.41629033  2.07429897 -0.54562351]
 [ 2.59455459  2.82896575 -0.19599961 -0.54562351  0.58413686]]
-----mat.dot(inv(mat))
[[ 1.00000000e+00  4.25060197e-15  5.44449844e-16  1.46585436e-15
  -3.17887458e-16]
 [-1.01348550e-14  1.00000000e+00 -4.93962174e-16 -2.42741991e-15
   1.21280219e-15]
 [-9.15105290e-15  2.67323017e-15  1.00000000e+00 -6.67063364e-16
  -7.14920612e-16]
 [ 2.11353448e-15  8.14671769e-16  7.99875283e-17  1.00000000e+00
  -3.79620488e-17]
 [-3.70507907e-15  3.36817506e-15 -6.43457141e-17 -2.90779435e-18
   1.00000000e+00]]
-----qr(mat)-----
[[-7.84722433  6.33808317  6.73324921 -3.30152419  3.62042646]
 [ 0.         -1.49681976  8.5187475  -5.19351382  6.23948695]
 [ 0.          0.         -6.0917072   0.84541523 -0.74963641]
 [ 0.          0.          0.         -0.87160841 -2.06438679]
 [ 0.          0.          0.          0.          0.25472357]]

表达式X.T.dot(X)计算X和它的转置X.T的点积。

三、伪随机数生成

numpy.random模块对Python内置的random进行了补充，增加了一些用于高效生成多种概率分布的样本值的函数。例如，你可以用normal来得到一个标准正态分布的4×4样本数组：

1
2
3
4
5
6
7
8
9
import numpy as np
samples = np.random.normal(size=(4,4)) # 这里必须加size
print(samples)

# 结果
[[ 0.86019211 -1.69311951  0.3569995   0.01053486]
 [ 0.73102713  0.35746827 -1.42414119  0.76492792]
 [-0.18037867 -0.58179004  3.34013328  0.26265435]
 [-2.37648991  1.02255001 -1.32951733 -2.02195647]]

这些都是伪随机数，是因为它们都是通过算法基于随机数生成器种子，在确定性的条件下生成的。你可以用NumPy的np.random.seed更改随机数生成种子：

1
2
···
np.random.seed(1234)

numpy.random的数据生成函数使用了全局的随机种子。要避免全局状态，你可以使用numpy.random.RandomState，创建一个与其它隔离的随机数生成器：

1
2
3
4
5
6
rng = np.random.RandomState(1234)
print(rng.randn(10))

# 结果
[ 0.47143516 -1.19097569  1.43270697 -0.3126519  -0.72058873  0.88716294
  0.85958841 -0.6365235   0.01569637 -2.24268495]

四、随机漫步

4.1 简单随机漫步

通过模拟随机漫步来说明如何运用数组运算。先来看一个简单的随机漫步的例子：从0开始，步长1和－1出现的概率相等。

下面是使用Python内置的random包进行的1000步的随机漫步：

1
2
3
4
5
6
7
8
9
10
11
12
import random
import matplotlib.pyplot as plt
position = 0
walk = [position]
steps = 1000
for i in range(steps):
    step = 1 if random.randint(0, 1) else -1
    position += step
    walk.append(position)

plt.plot(walk[:100])
plt.show()

根据前100个随机漫步值生成的折线图为：

不难看出，这其实就是随机漫步中各步的累计和，可以用一个数组运算来实现。因此，我用np.random模块一次性随机产生1000个“掷硬币”结果（即两个数中任选一个），将其分别设置为1或－1，然后计算累计和：

1
2
3
4
5
···
nsteps = 1000
draws = np.random.randint(0, 2, size=nsteps)
steps = np.where(draws > 0, 0, -1, 1)
walk = steps.cunsum()

可以沿着漫步路径做一些统计工作了，比如求取最大值和最小值：

1
2
3
···
print(walk.min())
print(walk.max())

看一个复杂点的统计任务——首次穿越时间，即随机漫步过程中第一次到达某个特定值的时间。假设我们想要知道本次随机漫步需要多久才能距离初始0点至少10步远（任一方向均可）。np.abs(walk)>=10可以得到一个布尔型数组，它表示的是距离是否达到或超过10，而我们想要知道的是第一个10或－10的索引。可以用argmax来解决这个问题，它返回的是该布尔型数组第一个最大值的索引（True就是最大值）：

1
2
3
···
print((np.abs(walk) >= 10).argmax())
37

4.2 一次模拟多个随机漫步

模拟多个随机漫步过程（比如5000个），只需对上面的代码做一点点修改即可生成所有的随机漫步过程。只要给numpy.random的函数传入一个二元元组就可以产生一个二维数组，然后我们就可以一次性计算5000个随机漫步过程（一行一个）的累计和了：

1
2
3
4
5
6
7
nwalks = 5000
nsteps = 1000
draws = np.random.randint(0,2,size=(nwalks, nsteps)) # 0 或者 1, 5000个1000步的随机漫步过程

steps = np.where(draws > 0, -1, 1)
walks = steps.cumsum(1)
print(walks)

计算所有随机漫步过程的最大值和最小值：

1
2
walks.max()
walks.min()

计算30或－30的最小穿越时间。这里稍微复杂些，因为不是5000个过程都到达了30。可以用any方法来对此进行检查：

1
2
3
hits30 = (np.abs(walks) >= 30).any(1)
print(hits30)
print(hits30.sum())

利用这个布尔型数组选出那些穿越了30（绝对值）的随机漫步（行），并调用argmax在轴1上获取穿越时间：

1
2
crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1) # 选出行
print(crossing_times.mean())

尝试用其他分布方式得到漫步数据。只需使用不同的随机数生成函数即可，如normal用于生成指定均值和标准差的正态分布数据：

1
steps = np.random.normal(loc=0, scale=0.25, size=(nwalks, nsteps))