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

问题描述

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

给定从原点开始的2D随机游动(在网格中)，到达直线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，此代码也很慢。我不确定如何对其进行优化或更改以正确模拟该问题。

---

解决方案

不是按顺序模拟一组随机游动，而是尝试同时模拟多条路径并跟踪这些路径发生的概率，例如，我们从位置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中的子句来考虑一维情况。