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

Question

Evangelos Mitsokapas 2017-10-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360148-generating-a-discrete-rw-on-a-unit-circle-and-calculate-its-long-time-distribution

评论： Evangelos Mitsokapas 2017-10-8

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-');

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2017-10-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360148-generating-a-discrete-rw-on-a-unit-circle-and-calculate-its-long-time-distribution#answer_284705

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

John D'Errico 2017-10-7

编辑：John D'Errico 2017-10-7

在 MATLAB Online 中打开

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.

Image Analyst 2017-10-8

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 2017-10-8

@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 2017-10-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360148-generating-a-discrete-rw-on-a-unit-circle-and-calculate-its-long-time-distribution#answer_284708

编辑：Image Analyst 2017-10-7

在 MATLAB Online 中打开

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]);