How would I start this problem?

We consider again the one-dimensional version of the drunkard problem. This time he starts out at position x (e.g., x=10), and takes unit steps left or right with equal probability. Suppose at 0 and a (a>x>0, e.g., a=30) are two bars. When he reaches either one, he is trapped forever.

We want to determine the average time t (the mean number of steps) before the drunkard is trapped. This is a diffusion problem and the 'time' t is called the first passage time. Do a simulation to compute t. (Again note that only the average of t has a meaning and it is that average that we want to compute.)

Answer 1

[Image Analyst]

See my attached collection of random walk programs/demos I've written. Adapt as needed.

[Gaëtan Poirier]

Would this also help with the edited part? Namely, finding that average number of steps before the random walker becomes trapped?

[Image Analyst]

Yes. Random_walk2 is probably the one you want to start adapting.

[Gaëtan Poirier]

How would I have it so the computer generates a random number that between 0 and 1 and this provides the equal probability of going left or right. Your code seems to have the user input a certain number of steps and calculate how far from the origin the particle goes, were I am finding the number of steps before the particle reaches another "bar."

[Walter Roberson]

Convert the "for" to a "while true". When a blockage is detected, break from the loop.

[Gaëtan Poirier]

Could you walk me through that? I'm really not good at compsci; it doesn't come intrinsic...

[Image Analyst]

The code has the particle keep moving until it hits the boundary, then it records the number of steps it took as well as the distance. Look at random_walk2.m and see this code:

% Record the number of steps 
% it took to get out of bounds.
out_of_bounds(experiment) = stepNumber;

That, out_of_bounds, is the number of steps it took to hit the boundary for each experiment.

You need to actually RUN random_walk2.m and check out the figures it makes. You'll see it gives a bar chart of the distribution of the number of steps. Isn't that what you want? Of course you could also take the mean, median, standard deviation, etc. of that array if you want. All you have to do is to delete or ignore the distance calculation if you don't care about distance.

[Gaëtan Poirier]

So, how would I put in the "second bar" in my problem? And to make it a one-dimensional problem, how would this effect the code?

[Gaëtan Poirier]

I just don't understand how to initially make code a random position, a random position for the second bar, and a 50% chance of going left or right...

[Image Analyst]

You need to compute x, and don't have any y. Then

if x(stepNumber, experiment) < a || x(stepNumber, experiment) > b
      % It's out of bounds, so log number of steps and exit this experiment.

[Gaëtan Poirier]

But this stepNumber still makes the user input a number for the steps the thing can take. How would I make this randomized?

[Gaëtan Poirier]

So I'm looking to write

r=rand();
x(1)=r; %initial position – random location
a=r>=x(1); %bar two which has to be greater than the initial position
if r>0.5
x=x+1;
elseif r<0.5
x=x-1:
end

This is what I think I have to write, but I don't know the syntax...

[Image Analyst]

Where exactly does it ask the user for the number of steps? It doesn't. stepNumber is a loop iterator. maxSteps in a parameter - you can set it very high to make sure that you never run out of steps before you reach the boundary. I've done that for you. See if you can finish it yourself. I mean, you got to do something for your homework, right?

% Monte Carlo experiment for random walk
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

% maxStepLength = max length of one step, from minus this to + this.
maxStepLength = 1;
% Define Monte Carlo parameters
numExperiments = 3000;
maxSteps = 10000;
% Initialize data.
x = zeros(maxSteps, numExperiments);

% Since there is only 1 row it should only
%record the fist time it crosses the boundary
out_of_bounds=zeros(1, numExperiments);

% Define a and b.
a = 0;
B = 30; % Boundary is 30 = max allowed distance before we go to next particle
% Better is to use inputdlg() to ask the user for these numbers.

for experiment = 1 : numExperiments
  for stepNumber = 2: maxSteps
    % Find the location of the next point
    deltaX = maxStepLength * (2 * rand(1) - 1);
    x(stepNumber, experiment) = x(stepNumber-1, experiment) + deltaX;

    if x(stepNumber, experiment) < a || x(stepNumber, experiment) >= B
      % Then the particle is out of bounds.
      % Record the number of steps 
      % it took to get out of bounds.
      out_of_bounds(experiment) = stepNumber;
      break; % Quit tracking this particle.
    end  
  end
end
maxNumberOfSteps = max(out_of_bounds);

% DONE!  Now make fancy plots.
% Plot the number of steps to get out of bounds.
subplot(1,2,1);
plot(out_of_bounds, 'b-', 'LineWidth', 2);
xlim([0, numExperiments]);
caption = sprintf('Steps until out of bounds for %d experiments', numExperiments);
title(caption, 'FontSize', fontSize);
xlabel('Experiment Number', 'FontSize', fontSize);
ylabel('Number of Steps Until Out of Bounds', 'FontSize', fontSize);
grid on;
subplot(1,2,2);
histogram(out_of_bounds);
grid on;
caption = sprintf('Distribution of the Number of Steps until out of bounds');
title(caption, 'FontSize', fontSize);
xlabel('Number of Steps', 'FontSize', fontSize);
ylabel('Count of Number of Experiments', 'FontSize', fontSize);

% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Gaëtan Poirier]

Much thanks!

[Gaëtan Poirier]

So, I decided to use code that was a bit simpler and more comprehendible (for my level). I have the following, although I have to run it once with the A(N) and ans statements as comments before it returns a number:

for N=1:10000
    A(N)=n;
    x(1)=10;
    for n=1:Inf
        g(n)=rand();
        if g(n)<0.5
            x(n+1)=x(n)-1;
        else x(n+1)=x(n)+1;
        end
        if x(n+1)>29||x(n+1)<1
            break
        else continue
        end
    end
end
num=(1/N)*sum(A);

this gives me a value around 200. If anyone knows how I can run this smoothly, that would be great. The second part of the assignment is to find distribution after exactly 20 or 21 steps. Half of calculations should be with 20 and the other half should be with 21. How would I change the code to find the distribution?

[Image Analyst]

Take out my comments, use simpler/more cryptic variable names, and only allow +/- 1 step sizes, and it's essentially my code. To get the distribution back again, add back in the line that you took out that say

out_of_bounds(experiment) = stepNumber;

then, after the loop

histogram(out_of_bounds);

It's almost never a good idea to take out comments. I suggest you put comments into your code.

By the way, instead of that if/then/else statement to get x, simply get rid of g and x and simply do this:

x(n+1) = x(n+1) + sign(rand(1)-0.5)

sign(rand(1)-0.5) will give you randomly +1 or -1.

[Gaëtan Poirier]

So here it is so far:

for N=1:10000 %number of "walkers"
    x(1)=10; %starting position
out_of_bounds = stepNumber;
for stepNumber=20;  %number of steps limited to each "walker"
    g(n)=rand(); %random number
    if g(n)<0.5
        x(n+1)=x(n)-1; %if random number is less than, go left
    else x(n+1)=x(n)+1; %else go right
    end
    if x(n+1)>29||x(n+1)<1 %second bar––if "walker" reaches 30, the loop terminates
        break
    else continue % else it continues 
    end
end 
end
histogram(out_of_bounds);

My distribution is a square, completely filled so I obviously have something wrong. Any ideas?

[Image Analyst]

Look at my code again. Do you see this:

for stepNumber=20

No. That's because your for loop will execute only once, and with a value of only 20, not every integer from 1 up to 20.