随机漫步模拟及概率分析

《利用Python进行数据分析》一书中提到一个随机漫步。同时，在求职过程中，也遇到过要求计算落入两端的概率之比的问题。正好可以用此案例进行模拟。这里整理出来，以享网友。

基本样例：

在一个无限的数轴上，一个小球从0开始随机移动，每一步以相同的概率向左或者向右移动1个单位。请模拟一下其过程

模拟一次过程

导入相关包

from random import randint
import matplotlib.pyplot as plt
import numpy as np

从零开始依次模拟

postition = 0
walk = [postition]

steps = 1000

for i in range(steps):
    step = 1 if randint(0, 1) else -1
    postition += step
    walk.append(postition)

plt.plot(walk)

利用Numpy库进行模拟

nsteps = 1000
# 左闭右开的区间
# 首先生成要每步的随机运动量
draws = np.random.randint(0, 2, size=nsteps)
# 如果是不同概率左右移动，则通过控制范围及比例，即可，不再描述
steps = np.where(draws > 0, 1, -1)
# 计算每步时到达的位置，用累计求和的方式
walk = steps.cumsum()

一次性模拟多个过程

通过传入参数，可以一次性模拟多个过程

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.shape

(5000, 1000)

# 计算下所有过程的最大值及最小值
walks.min()

-114

walks.max()

125

分析是否达到某个点或者某个距离上

# 计算下每个过程是否达到了某个点，比如30(每个过程中，只要有一个点达到即算)
hists30 = (np.abs(walks) >= 30).any(1)
hists30

array([ True, False,  True, ...,  True,  True,  True])

# 计算达到30这个距离的数量
hists30.sum()

3354

# 什么时候穿越了30这个点呢？用argmax在轴1上获取， argmax返回第一个最大值的索引
# 掌握这种切片方式

crossing_times = (np.abs(walks[hists30]) >= 30).argmax(1)
crossing_times

array([647, 327, 663, ..., 855, 709, 913], dtype=int64)

无限运行，最终会在哪个位置

直觉一下，应该会回到原点。
看一下程序执行的结果，即求最终位置的平均值

# 直觉应该是回到原点
walks[:,-1].mean()

0.0584

与直觉相符

如果两侧有坑，计算落入坑的概率比

两侧坑等距

先模拟一下过程

nwalks = 5000
nsteps = 1000
# 左闭右开的区间
draws = np.random.randint(0, 2, size=(nwalks, nsteps))
# 如果是不同概率左右移动，则通过控制范围及比例，即可，不再描述
steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(1)

筛选出落入坑的过程

# 计算落入坑的
a = 20
all_in = (np.abs(walks) >= a).any(1)
all_in.sum()

4704

# 在落入坑中的里面判断是左侧还是右侧
# 注意是要找第一个进坑的
# 判断进左右坑的穿越时间，如果为0 ，则表明是落到另一侧；如果都不为0，则落到索引较小的那一侧
fall_in_hole = walks[all_in]
fall_in_hole.shape

(4704, 1000)

判断落入左侧还是右侧

基本思路：对任何一个过程，以通过以上的实现方式过程中，累计和有三种情况：只存在小于左侧距离的值，只存在大于右侧距离的值，同时存大小于左侧或右侧距离的值。
分别求其落入左侧坑和落入右侧坑的索引，用argmax求解。如果是情况1、2，则必定落入索引不为0的那一侧（因考虑实际情况，第一个动作不会入坑-或者更完善的情况，可以在所有数据添加一个点，从为开始），如果是情况3，则取索引较小的那个

# 判断落入左侧的索引
cross_left = (fall_in_hole <= - a).argmax(1)
cross_left_flag = np.where(cross_left == 0, 0, 1)
cross_left_flag.sum()

2687

# 判断落入右侧的索引
cross_right = (fall_in_hole >= a).argmax(1)
cross_right_flag = np.where(cross_right == 0, 0, 1)
cross_right_flag.sum()

2591

# 输出到excel中看看计算是否准确
def get_left_right_num(cross_left, cross_right):
    # left right 分别表示落入左右侧的个数 res为索引（用于输出实际序列用）
    left = 0
    right = 0
    res = []
    # 最终的序列索引
    for left_index, right_index in zip(cross_left, cross_right):
        if left_index == 0:
            right += 1
            res.append(right_index)
        elif right_index == 0:
            left += 1
            res.append(left_index)
        else:
            if left_index > right_index:
                right += 1
                res.append(right_index)
            else:
                left += 1
                res.append(left_index)
    print(left, right)
    return left, right, res
left, right, res = get_left_right_num(cross_left, cross_right)
print(left / right)   

两侧坑不等距

nwalks = 5000
nsteps = 1000
# 左闭右开的区间
draws = np.random.randint(0, 2, size=(nwalks, nsteps))
# 如果是不同概率左右移动，则通过控制范围及比例，即可，不再描述
steps = np.where(draws > 0, 1, -1)
walks = steps.cumsum(1)

# 控制比例，这里是1：9
a = 5
left_dis = 1 * a
right_dis = 9 * a

# 计算落入坑的
# 小心：不能用or
all_in = ((walks <= - left_dis) | (walks >= right_dis)).any(1)
all_in.sum()

4733

fall_in_hole = walks[all_in]

cross_left = (fall_in_hole <= - left_dis).argmax(1)
cross_left_flag = np.where(cross_left == 0, 0, 1)
cross_left_flag.sum()

4360

cross_right = (fall_in_hole >= right_dis).argmax(1)
cross_right_flag = np.where(cross_right == 0, 0, 1)
cross_right_flag.sum()

798

left, right, res = get_left_right_num(cross_left, cross_right)
print(left / right)

4357 376
11.587765957446809

# 根据最终的序列索引，生成想要的序列（即一旦入坑，后面的值不再发生变动）
out = []
row, col = fall_in_hole.shape
fall_in_hole_2 = fall_in_hole.copy()
for index, fall_index in enumerate(res):
    fall_in_hole_2[index, fall_index:] = fall_in_hole_2[index, fall_index]
fall_in_hole_2

array([[ 1,  2,  3, ..., -5, -5, -5],
       [-1, -2, -1, ..., -5, -5, -5],
       [-1, -2, -1, ..., -5, -5, -5],
       ...,
       [ 1,  0, -1, ..., -5, -5, -5],
       [ 1,  2,  3, ..., -5, -5, -5],
       [-1,  0,  1, ..., -5, -5, -5]], dtype=int32)

# 保存为文件，查看不同
np.savetxt("fall_in_hole.txt",fall_in_hole)
np.savetxt("fall_in_hole_2.txt", fall_in_hole_2)