使用 NumPy 来模拟随机游走（Random Walk）

考虑一种最简单的情形，即从位置 0 开始，每次前进一步（step=1）或者后退一步（step=-1），前进与后退的概率相同。

不使用 NumPy，直接使用 Python 内置的 random 函数：

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()
1
2
3
4
5
6
7
8
9
10
11
12
13

其实，walk 的值其实就是 step 的累加值。我们考虑使用 NumPy 来实现：

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

我们可以注意到两个 randint 函数的区别：

• Python 的内置 random 模块每次只会采样一个值，而 numpy.random 可以同时采样一系列值。这意味着在产生大量样本时，numpy.random 的速度会比 random 模块快上几个量级
• random.randint 取值范围：[low, high]；numpy.random.randint 取值范围：[low, high)

我们可以很轻易的得到 walk 的最大值或者最小值：

walk.min()
"""
-7
"""
1
2
3
4

如果我们想要找到 the first crossing time，即第一次穿过某个位置的 step，例如，我们想找第一次离原点 0 距离为 10 的 step：

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

np.abs(walk) >= 10 返回一个布尔数组，我们用 argmax 方法找到第一个值为 True 的位置（True 即为最大值，且argmax 会返回第一个最大值出现的位置）。

---

上面我们考虑了模拟一次随机游走的流程，如果我们想模拟很多次，例如，5000次，那么使用 NumPy 也可以很容易实现：

nwalks = 5000
nsteps = 1000

draws = np.random.randint(0, 2, size=(nwalks, nsteps))

steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(1) # 计算每行，也就是每次模拟的累加值

walks
"""
array([[  1,   2,   3, ...,  52,  51,  50],
       [  1,   2,   1, ...,  20,  19,  20],
       [  1,   2,   3, ..., -34, -35, -36],
       ...,
       [  1,   2,   3, ...,  44,  45,  46],
       [ -1,   0,  -1, ...,   0,  -1,   0],
       [  1,   2,   1, ...,  48,  49,  50]])
"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

我们计算这 5000 次模拟中有多少次模拟穿过了 30 或 -30 位置：

hits30 = (np.abs(walks) >= 30).any(1)
hits30.sum()
"""
3361
"""
1
2
3
4
5

进一步，我们可以得出这 3361 次模拟中第一次穿过 30 或 -30 位置的平均 step：

crossing_time = (np.abs(walks[hits30]) >= 30).argmax(1)
crossing_time.mean()
"""
500.99226420708123
"""
1
2
3
4
5

---

References

Python for Data Analysis, 2${}^{\mathrm{n}\mathrm{d}}$ edition. Wes McKinney.