Average Time To Hit A Given Line On 2d Random Walk On A Unit Grid

I am trying to simulate the following problem: Given a 2D random walk (in a lattice grid) starting from the origin what is the average waiting time to hit the line y=1-x import num

Solution 1:

Instead of simulating a bunch of random walks sequentially, lets try simulating multiple paths at the same time and tracking the probabilities of those happening, for instance we start at position 0 with probability 1:

states = {0+0j: 1}

and the possible moves along with their associated probabilities would be something like this:

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

With this construct we can update to new states by going over the combination of each state and each move and update probabilities accordingly

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

Then we just iterate this popping out all winning states at each step:

for stepIdx inrange(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 inset(states.items()):
        # if y = 1-xif 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}")

you would also be able to look at states for an idea of the distribution of possible positions, although totalling it in terms of distance from the line simplifies the data. Anyway, the final step would be to average the steps taken by the probability of taking that many steps and see if it converges:

total_average_num_moves += stepIdx * winning_chances

But we might be able to gather more insight by using symbolic variables! (note I'm simplifying this to a 1D problem which I describe how at the bottom)

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
}

This with the exact code as written above gives us this sequence:

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

And if we ask the OEIS what the sequence 1,2,5,14,42,132... means it tells us those are Catalan numbers with the formula of (2n)!/(n!(n+1)!) so we can write a function for the non-zero terms in that series as:

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

or in actual code:

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)

which now gives us a relatively instant way of calculating relevant parameters for any length, or more usefully ask wolfram alpha what the average would be (it diverges)

---

Note that we can simplify this to a 1D problem by considering y-x as one variable: "we start at y-x = 0 and move such that y-x either increases or decreases by 1 each move with equal chance and we are interested when y-x = 1. This means we can consider the 1D case by subbing in z=y-x.

Solution 2:

Vectorisation would result in much faster code, approximately ~90K times faster. Here is the function that would return step to hit y=1-x line starting from (0,0) and trajectory generation on the 2D grid with unit steps .

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 = -1if 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

if number_of_steps_first_hit is -1 it means trajectory does not hit the line.

A longer simulation and repeating might give the average behaviour, but the following one tells if it does not escape to Infiniti it hits line on average ~84 steps.

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

Much longer sim_steps will change the result though.

PS: Good exercise, hope that this wasn't a homework, if it was homework, please cite this answer if it is used.

Edit As discussed in the comments current _random_walk_2D works for 8-directions. To restrict it to cardinal direction we could do the following filtering:

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]

though this would slow it down the function a bit but still will be super fast compare to for loop solutions.