4.7 示例：随机漫步

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

下面是一个通过内置的random模块以纯Python的方式实现1000步的随机漫步：

  1. In [247]: import random
  2.    .....: position = 0
  3.    .....: walk = [position]
  4.    .....: steps = 1000
  5.    .....: for i in range(steps):
  6.    .....:     step = 1 if random.randint(0, 1) else -1
  7.    .....:     position += step
  8.    .....:     walk.append(position)
  9.    .....:

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

  1. In [249]: plt.plot(walk[:100])

图4-4 简单的随机漫步

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

  1. In [251]: nsteps = 1000
  3. In [252]: draws = np.random.randint(0, 2, size=nsteps)
  5. In [253]: steps = np.where(draws > 0, 1, -1)
  7. In [254]: walk = steps.cumsum()

有了这些数据之后，我们就可以沿着漫步路径做一些统计工作了，比如求取最大值和最小值：

  1. In [255]: walk.min()
  2. Out[255]: -3
  4. In [256]: walk.max()
  5. Out[256]: 31

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

  1. In [257]: (np.abs(walk) >= 10).argmax()
  2. Out[257]: 37

注意，这里使用argmax并不是很高效，因为它无论如何都会对数组进行完全扫描。在本例中，只要发现了一个True，那我们就知道它是个最大值了。

一次模拟多个随机漫步

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

  1. In [258]: nwalks = 5000
  3. In [259]: nsteps = 1000
  5. In [260]: draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1
  7. In [261]: steps = np.where(draws > 0, 1, -1)
  9. In [262]: walks = steps.cumsum(1)
  11. In [263]: walks
  12. Out[263]: 
  13. array([[  1,   0,   1, ...,   8,   7,   8],
  14.        [  1,   0,  -1, ...,  34,  33,  32],
  15.        [  1,   0,  -1, ...,   4,   5,   4],
  16.        ..., 
  17.        [  1,   2,   1, ...,  24,  25,  26],
  18.        [  1,   2,   3, ...,  14,  13,  14],
  19.        [ -1,  -2,  -3, ..., -24, -23, -22]])

现在，我们来计算所有随机漫步过程的最大值和最小值：

  1. In [264]: walks.max()
  2. Out[264]: 138
  4. In [265]: walks.min()
  5. Out[265]: -133

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

  1. In [266]: hits30 = (np.abs(walks) >= 30).any(1)
  3. In [267]: hits30
  4. Out[267]: array([False,  True, False, ..., False,  True, False], dtype=bool)
  6. In [268]: hits30.sum() # Number that hit 30 or -30
  7. Out[268]: 3410

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

  1. In [269]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)
  3. In [270]: crossing_times.mean()
  4. Out[270]: 498.88973607038122

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

  1. In [271]: steps = np.random.normal(loc=0, scale=0.25,
  2.    .....:                          size=(nwalks, nsteps))