在单位格网上的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}")

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(它有分歧)

这篇关于在单位格网上的2D随机游动中命中给定线的平均时间的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，WP2