初识数据分析之NumPy 笔记七 随机漫步

来源:《利用Python进行数据分析·第2版》

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

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

In [146]: position = 0

In [147]: walk = [position]

In [148]: steps = 10000

In [149]: for _ in range(steps):
     ...:     step = 1 if random.randint(0, 1) else -1
     ...:     position += step
     ...:     walk.append(position)
     ...:

In [150]: plt.plot(walk)

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

In [164]: plt.plot(walk[:100])
Out[164]: [<matplotlib.lines.Line2D at 0x1b114b83c08>]

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

In [165]: nsteps = 1000

In [166]: draws = np.random.randint(0, 2, size=nsteps)

In [168]: steps = np.where(draws > 0, 1, -1)

In [169]: walk = steps.cumsum()

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

In [173]: walk.min()
Out[173]: -9

In [174]: walk.max()
Out[174]: 60

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

In [178]: (np.abs(walk) >= 10).argmax()
Out[178]: 297

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

一次模拟多个随机漫步

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

In [185]: nwalks = 5000

In [186]: nsteps = 1000

In [187]: draws = np.random.randint(0, 2, size=(nwalks, nsteps))

In [189]: steps = np.where(draws > 0, 1, -1)

In [191]: walks = steps.cumsum(1)

In [192]: walks
Out[192]:
array([[  1,   2,   3, ...,  46,  47,  46],
       [  1,   0,   1, ...,  40,  41,  42],
       [  1,   2,   3, ..., -26, -27, -28],
       ...,
       [  1,   0,   1, ...,  64,  65,  66],
       [  1,   2,   1, ...,   2,   1,   0],
       [ -1,  -2,  -3, ...,  32,  33,  34]], dtype=int32)

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

In [193]: walks.max()
Out[193]: 122

In [195]: walks.min()
Out[195]: -128

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

In [202]: hits30 = (np.abs(walks) >=30).any(1)

In [203]: hits30
Out[203]: array([ True,  True,  True, ...,  True, False,  True])

In [204]: hits30.sum()
Out[204]: 3368

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

In [206]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)

In [207]: crossing_times.mean()
Out[207]: 509.99762470308787

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

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

In [209]: steps
Out[209]:
array([[-0.20472197, -0.03407671, -0.24781626, ...,  0.22195392,
         0.21250584,  0.43994277],
       [-0.09066595,  0.04741116,  0.09872352, ...,  0.10684413,
         0.2540548 , -0.22255203],
       [ 0.00942324,  0.03418277, -0.42305453, ..., -0.05711872,
        -0.11820873, -0.71117832],
       ...,
       [ 0.1230118 , -0.18155722,  0.12782267, ...,  0.0730254 ,
         0.06546779,  0.47920034],
       [-0.20550269,  0.16459572, -0.03138668, ...,  0.06915243,
         0.07831634,  0.27554056],
       [-0.23056172,  0.28924403,  0.08444807, ..., -0.27018412,
        -0.02140732,  0.06343304]])

 

  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值