Generating a discrete RW on a unit circle and calculate its long-time distribution

Hi, I am new to MATLAB and I would like to generate a random walk of a particle on a ring (circle of radius r). For each iteration, the particle should move to the "right" on the circle by 1, with a probability p=1/2 and to the "left" on the circle by -1 with a probability q=1-p=1/2. I can generate a random walk in MATLAB like in the following code, but I am not sure how to confine it on a circle. After this, I would like to calculate the long-time distribution of the particle numerically.

Thanks in advance.

 %Random Walk in 2D

 close all;
 clear all;

 nsteps=10;

 x(1)=0;
 y(1)=0;

 for i=1:nsteps

    theta=2*pi*rand()
    r=1;

    dx=r*cos(theta);
    dy=r*sin(theta);

    x(i+1)=x(i)+dx;
    y(i+1)=(x(i)/tan(theta))+dy;
end

 plot(x,y,'r-');

Answer 1

[John D'Errico]

Ok. New to MATLAB or not, I'll give you credit for having made a credible effort.

So what are you doing wrong here? You are looking at it from the wrong direction, trying to solve the problem in a bass-ackwards way. What does that mean?

This is NOT a random walk in the (x,y) plane! This is a random walk in polar angle theta, and the (x,y) coordinates are computed parameters as a function of that polar angle theta.

What you are trying to do is make it a walk in the plane, then constrain the points to lie on a circle. That is coming at the problem the wrong way around.

So just generate the angle theta as a sequence, then compute x and y. This ensures that the radius is ALWAYS 1, so the points ALWAYS lie exactly on the unit circle.

r = 1;
ntheta = 100000;
theta0 = 0;
dtheta = rand(1,ntheta-1)*2 - 1;
theta = cumsum([theta0,dtheta]);

x = r*cos(theta);
y = r*sin(theta);

plot(x,y,'.')
grid on
axis equal
axis square

I'll let you do the plots and the computations.

As far as computing the long term distribution, consider the difference between these two calls to hist:

hist(theta,100)
hist(mod(theta,2*pi),100)

Which is appropriate for your question?

Finally, could you do this using a loop? Of course. But you need to learn to use MATLAB, as it was designed to be used. To solve the problem using a loop, just create theta incrementally. Then at the end, compute x and y.

[Evangelos Mitsokapas]

Hi, thank you for taking the time to answer. It appears that you are right, I was just making the random walker to transpose on the plane, but what I really wanted is to transform the movement on the x,y-plane to the polar one. Thank you for that. The polar plot looks as expected. I believe that the first histogram is proper for my case. Two more questions: a) How is it possible to normalize this? b) Is there a way for MATLAB to somehow fit this data to an approximate analytical distribution? (I have seen this being done, but it looks quite a complicated thing to code). Again, thanks!

[Evangelos Mitsokapas]

Oh and an another question: While I perfectly understand that with dtheta = rand(1,ntheta-1)*2 - 1; you are generating a 1xN array with random numbers uniformly distributed in the interval (a,b)=(-1,1), I am not yet able to understand why do you then define theta as a cumsum. What is the reason for it? I know that you need to do that because I plotted the corresponding RW with dtheta instead of theta and I got a hyperbola instead of a circle.

[John D'Errico]

DOH. I was asleep when I wrote the rand step. If you want a +/-1 binomial random variable for a step, I forgot to add one thing.

dtheta = round(rand(1,ntheta-1))*2 - 1;

So it comes down to this:

r = 1;
ntheta = 10000000;
theta0 = 0;
dtheta = round(rand(1,ntheta-1))*2 - 1;
theta = cumsum([theta0,dtheta]);
x = r*cos(theta);
y = r*sin(theta);

See that I've used a considerably larger sample here, so the randomness can get smoothed out.

As far as the histogram goes, actually, the second histogram is what you want. In fact, the expected distribution of theta will be normally distributed over the long term. But the distribution of the points around the circle will be uniform!

hist(theta,1000)

hist(mod(theta,2*pi),100)

Why is this? Because theta+k*(2*pi) is equivalent to theta, when theta is the polar coordinate angle, and integer k. The mod computation removes multiples of 2*pi from the angle.

So the law of large numbers ensures that theta itself is normally distributed. It is just the sum of zillions of iid random variables. The result will be normally distributed.

BUT what matters is where the points lie around the circle. And that will be uniform. It will take many samples to make that happen.

[Evangelos Mitsokapas]

Thanks for coming back. I see, now rand() is a uniform random variable in the interval [-1,1] and you round the result to get the step +/-1. Smart! Thank you for this. Btw, I am not familiar with MATLAB hist(*,*) therefore, could you please explain what the axes of the histogram represent in this case? And what about a normalization? Would that be possible?

[Image Analyst]

The vertical axis is the frequency (count, number of occurrences). The x axis is the value. So it says how many times each number occurs. You can use the newer function histogram() instead of hist() if you have a recent release. It has several normalization options, like 'count' (default), 'probability', 'countdensity', 'pdf', 'cumcount', and 'cdf'.

[Evangelos Mitsokapas]

@Image Analyst: Thank you for the answer! Now it becomes more clear what has been going on behind the histogram. Would you be so kind as to provide an example of normalization of this histogram using the histogram() command? Many thanks.

Answer 2

[Image Analyst]

Your first step is right on the edge of the circle! So you need to have a different step length than the size of the outer, limiting circle. Then you need to computer y(i+1) differently (just add dy to it) and have a while loop to try again if the point would land outside the outer circle.

numberOfSteps = 2000;
maxIterations = 1000;
x(1) = 0;
y(1) = 0;
r=15; % Outer radius
stepLength = 1;

for k = 1 : numberOfSteps
  again = 1;
  loopCount = 1;
  while again && loopCount < maxIterations
    theta=2*pi*rand();
    dx = stepLength * cos(theta);
    dy = stepLength * sin(theta);
    
    x(k+1) = x(k) + dx;
    y(k+1) = y(k) + dy;
    if x(k+1)^2 + y(k+1)^2 <= r^2
      again = 0;
    end
    loopCount = loopCount + 1;
  end
end

plot(x,y,'r.-', 'MarkerSize', 18);
grid on;
xlabel('X');
ylabel('Y');
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.FontSize = 30;
ax.FontWeight = 'Bold';
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);

[Evangelos Mitsokapas]

Hi, thank you for taking the time to help. In the piece of code which you provide, the particle is able to move also within the circle, isn't it right? But what I want, is for the particle to move exactly on the unit circle, meaning on the circumference 1=sqrt(x^2+y^2). Correct me if I am understanding your code wrong, I am not yet familiar with Matlab :)

[Image Analyst]

Oh, sorry - I misunderstood.

[Evangelos Mitsokapas]

It's totally fine! Again, many thanks for your time :)