在单位网格上的 2D 随机游走中击中给定线的平均时间

Question

我正在尝试模拟以下问题：

给定从原点开始的二维随机游走（在格子网格中），到达直线 y=1-x 的平均等待时间是多少

import numpy as np
from tqdm import tqdm
N=5*10**3
results=[]
for _ in tqdm(range(N)):
  current = [0,0]
  step=0
  while (current[1]+current[0] != 1):
    step += 1
    a = np.random.randint(0,4)
    if (a==0):
      current[0] += 1
    elif (a==1):
      current[0] -= 1
    elif (a==2):
      current[1] += 1
    elif (a==3):
      current[1] -= 1
  results.append(step)

此代码即使对于 N<10**4 也很慢我不确定如何优化或更改它以正确模拟问题。

Answer 1

与其按顺序模拟一堆随机游走，不如尝试同时模拟多条路径并跟踪这些路径发生的概率，例如我们从位置 0 开始，概率为 1:

states = {0+0j: 1}

可能的动作及其相关概率如下所示：

moves = {1+0j: 0.25,  0+1j: 0.25, -1+0j: 0.25, 0-1j: 0.25}
# moves = {1: 0.5, -1:0.5} # this would basically be equivelent

使用此构造，我们可以通过遍历每个状态和每个移动的组合并相应地更新概率来更新到新状态

def simulate_one_step(current_states):
    newStates = {}
    for cur_pos, prob_of_being_here in current_states.items():
        for movement_dist,prob_of_moving_this_way in moves.items():
            newStates.setdefault(cur_pos+movement_dist, 0)
            newStates[cur_pos+movement_dist] += prob_of_being_here*prob_of_moving_this_way
    return newStates

然后我们只是迭代这个在每一步弹出所有获胜状态：

for stepIdx in range(1, 100):
    states = simulate_one_step(states)
    winning_chances = 0
    # use set(keys) to make copy so we can delete cases out of states as we go.
    for pos, prob in set(states.items()):
        # if y = 1-x
        if pos.imag == 1 - pos.real:
            winning_chances += prob
            # we no longer consider this a state that propogated because the path stops here.
            del states[pos]
    
    print(f"probability of winning after {stepIdx} moves is: {winning_chances}")

您还可以查看 states 以了解可能位置的分布，尽管根据与线的距离对其进行总计可以简化数据。无论如何，最后一步是将采取那么多步骤的概率平均所采取的步骤，看看它是否收敛：

total_average_num_moves += stepIdx * winning_chances

但我们也许可以通过使用符号变量来收集更多信息！（注意我将其简化为一维问题，我在底部描述了如何）

import sympy
x = sympy.Symbol("x") # will sub in 1/2 later
moves = {
    1: x, # assume x is the chances for us to move towards the target
   -1: 1-x # and therefore 1-x is the chance of moving away
}

这与上面编写的确切代码给了我们这个序列：

probability of winning after 1 moves is: x
probability of winning after 2 moves is: 0
probability of winning after 3 moves is: x**2*(1 - x)
probability of winning after 4 moves is: 0
probability of winning after 5 moves is: 2*x**3*(1 - x)**2
probability of winning after 6 moves is: 0
probability of winning after 7 moves is: 5*x**4*(1 - x)**3
probability of winning after 8 moves is: 0
probability of winning after 9 moves is: 14*x**5*(1 - x)**4
probability of winning after 10 moves is: 0
probability of winning after 11 moves is: 42*x**6*(1 - x)**5
probability of winning after 12 moves is: 0
probability of winning after 13 moves is: 132*x**7*(1 - x)**6

如果我们用 (2n)!/(n!(n+1)!) 的公式询问 the OEIS what the sequence 1,2,5,14,42,132... means it tells us those are Catalan numbers 那么我们可以为该系列中的非零项编写一个函数：

f(n,x) = (2n)! / (n! * (n+1)!) * x^(n+1) * (1-x)^n

或在实际代码中：

import math
def probability_of_winning_after_2n_plus_1_steps(n, prob_of_moving_forward = 0.5):
    return (math.factorial(2*n)/math.factorial(n)/math.factorial(n+1) 
           * prob_of_moving_forward**(n+1) * (1-prob_of_moving_forward)**n)

现在为我们提供了一种相对即时的方法来计算任何长度的相关参数，或者更有用 ask wolfram alpha what the average would be（发散）

请注意，我们可以通过将 y-x 视为一个变量来将其简化为一维问题：“我们从 y-x = 0 开始移动，使得 y-x 增加或减少 1每次移动的机会均等，我们对 y-x = 1 感兴趣。这意味着我们可以通过代入 z=y-x.

来考虑一维情况

Answer 2

矢量化会产生更快的代码，大约快 90K 倍。这是一个函数，它将 return 从 (0,0) 开始点击 y=1-x 线，并在二维网格上以单位步长生成轨迹。

import numpy as np

def _random_walk_2D(sim_steps):
    """ Walk on 2D unit steps
        return  x_sim, y_sim, trajectory, number_of_steps_first_hit to y=1-x """
    random_moves_x = np.insert(np.random.choice([1,0,-1], sim_steps), 0, 0)
    random_moves_y = np.insert(np.random.choice([1,0,-1], sim_steps), 0, 0)
    x_sim = np.cumsum(random_moves_x)
    y_sim = np.cumsum(random_moves_y)
    trajectory = np.array((x_sim,y_sim)).T
    y_hat = 1-x_sim # checking if hit y=1-x
    y_hit = y_hat-y_sim
    hit_steps = np.where(y_hit == 0)
    number_of_steps_first_hit = -1
    if hit_steps[0].shape[0] > 0:
        number_of_steps_first_hit = hit_steps[0][0]
    return x_sim, y_sim, trajectory, number_of_steps_first_hit

如果number_of_steps_first_hit为-1表示弹道未命中直线

更长的模拟和重复可能会给出平均行为，但下面的行为表明它是否没有逃脱到英菲尼迪，它平均达到约 84 步的直线。

sim_steps= 5*10**3 # 5K steps
#Repeat 
nrepeat = 40000
hit_step = [_random_walk_2D(sim_steps)[3] for _ in range(nrepeat)]
hit_step = [h for h in hit_step if h > -1]
np.mean(hit_step) #  ~84 step

更长的时间 sim_steps 会改变结果。

PS: 很好的练习，希望这不是作业，如果是作业，请引用此答案。

编辑正如评论中所讨论的，当前 _random_walk_2D 适用于 8 个方向。要将其限制在基本方向，我们可以进行以下过滤：

cardinal_x_y = [(t[0], t[1]) for t in zip(random_moves_x, random_moves_y)  
                if np.abs(t[0]) != np.abs(t[1])]
random_moves_x = [t[0] for t in cardinal_x_y]
random_moves_y = [t[1] for t in cardinal_x_y]

虽然这会稍微降低函数的速度，但与 for loop 解决方案相比仍然非常快。

在单位网格上的 2D 随机游走中击中给定线的平均时间

Average time to hit a given line on 2D random walk on a unit grid

python

simulation

probability

discrete-mathematics